004_misc_suppmat.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__purpose__ : get some numbers for tables and supplement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "##################################################################\n",
      "Current commit of pythonimports:\n",
      "commit b1d8bd7312fbf3c6afef4ad9ea2585831ec509a5\n",
      "Author: Brandon Lind <lindb@vcu.edu>\n",
      "Date:   Fri Feb 12 12:21:51 2021 -0500\n",
      "Today:\tMarch 10, 2021 - 11:29:52\n",
      "python version: 3.8.5\n",
      "##################################################################\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from pythonimports import *\n",
    "latest_commit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# determine megaSNP intersection with baseline-filtered SNPs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### load baseline-filtered SNPs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gatk\n",
      "\t JP\n",
      "\t DF\n",
      "varscan\n",
      "\t JP\n",
      "\t DF\n"
     ]
    }
   ],
   "source": [
    "# code from 001_testdata_explore.ipynb\n",
    "# get the file names\n",
    "sppfiles = {'gatk':{'JP': '/data/projects/pool_seq/non-pangenome/gatk_diploid_testdata/JP_i101-gatk/filtered_snps/JP_i101_filtered_concatenated_snps_max-missing_table_biallelic-only_translated.txt',\n",
    "                   'DF': '/data/projects/pool_seq/non-pangenome/gatk_diploid_testdata/DF_i52-gatk/filtered_snps/DF_i52_filtered_concatenated_snps_max-missing_table_biallelic-only_p000_translated.txt'},\n",
    "           'varscan':{'JP': '/data/projects/pool_seq/non-pangenome/varscan_pooled/JP_pooled/snpsANDindels/JP_pooled-varscan_all_bedfiles_SNP_translated.txt',\n",
    "                      'DF': '/data/projects/pool_seq/non-pangenome/varscan_pooled/DF_p52/snpsANDindels/DF_p52-varscan_all_bedfiles_SNP.txt'}}\n",
    "# read in the tables\n",
    "tablecols = ['CHROM','POS','REF','ALT','AF','QUAL','TYPE','FILTER','ADP','WT','HET','HOM','NC']\n",
    "snps = {}\n",
    "for method in sppfiles:\n",
    "    print(method)\n",
    "    snps[method] = {}\n",
    "    for spp in sppfiles[method]:\n",
    "        print('\\t',spp)\n",
    "        snps[method][spp] = pd.read_table(sppfiles[method][spp])\n",
    "        if method == 'varscan' and spp == 'JP':\n",
    "            # i ran all pools, isolate JP_p101 calls\n",
    "            locuscols = [col for col in snps[method][spp] if 'locus' in col]\n",
    "            snps[method][spp] = snps[method][spp][tablecols + ['JP_p101.GT','JP_p101.GQ','JP_p101.SDP',\n",
    "                                                               'JP_p101.DP','JP_p101.FREQ','JP_p101.PVAL'] +locuscols]\n",
    "            # remove AF=0 (because SNPs were called across all pools, JP_p101 can have AF=0)\n",
    "            # AF=0 wouldn't be possible if calling AFs on a single pool\n",
    "            snps[method][spp] = snps[method][spp][snps[method][spp]['JP_p101.FREQ'] != \"0%\"]\n",
    "            # remove rows with no calls\n",
    "            snps[method][spp] = snps[method][spp][snps[method][spp]['JP_p101.FREQ'] == snps[method][spp]['JP_p101.FREQ']]\n",
    "            # fill in AF column\n",
    "            snps[method][spp]['AF'] = snps[method][spp]['JP_p101.FREQ'].str.replace(\"%\", \"\").astype(float)/100\n",
    "        elif method =='varscan' and spp == 'DF':\n",
    "            df = snps[method][spp]  # view, will changes go to snps['varscan']['DF']\n",
    "            # add in a \">\" so locus names match between gatk and varscan\n",
    "            df['CHROM'] = [f\">{chrom}\" for chrom in df['CHROM'].tolist()]\n",
    "            df['locus'] = [f\">{locus}\" for locus in df['locus'].tolist()]\n",
    "# note varscan will output AF=0 (fixed for REF), I'm leaving these in so we can compare AFs of these too"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gatk JP 377080\n",
      "gatk DF 1526554\n",
      "varscan JP 3686528\n",
      "varscan DF 1601285\n"
     ]
    }
   ],
   "source": [
    "# remove noREF SNPs\n",
    "for method,spp in snps.items():\n",
    "    for sp,df in spp.items():\n",
    "        if method == 'gatk':\n",
    "            snps[method][sp] = df[~df['ALT'].str.contains('+', regex=False)]\n",
    "        print(method, sp, nrow(snps[method][sp]))\n",
    "        \n",
    "# original submission\n",
    "# gatk JP 377080\n",
    "# gatk DF 1526554\n",
    "# varscan JP 3686528\n",
    "# varscan DF 1601285"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "JP 32751\n",
      "DF 398774\n"
     ]
    }
   ],
   "source": [
    "# code from 003_testdata_validate_megaSNPs.ipynb\n",
    "# load megaSNPs\n",
    "megafiles = {'JP': '/data/projects/pool_seq/non-pangenome/varscan_mega/JP_RFmg7/snpsANDindels/JP_RFmg7-varscan_all_bedfiles_SNP_paralog_snps_translated.txt',\n",
    "             'DF': '/data/projects/pool_seq/non-pangenome/varscan_mega/DF_megaSNPs_from_download/DF_mega/snpsANDindels/02_baseline_filtered/DF_mega-varscan_all_bedfiles_SNP_paralog_snps.txt'}\n",
    "megasnps = {}\n",
    "for sp,file in megafiles.items():\n",
    "    megasnps[sp] = pd.read_table(file)\n",
    "    locuscol = 'unstitched_locus' if 'unstitched_locus' in megasnps[sp].columns else 'locus'\n",
    "    if locuscol == 'locus':\n",
    "        loci = [\">%s\" % locus for locus in megasnps[sp][locuscol].tolist()]\n",
    "    else:\n",
    "        loci = megasnps[sp][locuscol].tolist()\n",
    "    megasnps[sp][locuscol] = loci\n",
    "    megasnps[sp].index = loci\n",
    "    print(sp, nrow(megasnps[sp]))\n",
    "\n",
    "# original submission\n",
    "# JP 32751\n",
    "# DF 398774"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gatk JP 7408\n",
      "gatk DF 293\n",
      "varscan JP 25500\n",
      "varscan DF 825\n"
     ]
    }
   ],
   "source": [
    "# get intersection\n",
    "for method,spp in snps.items():\n",
    "    for sp,df in spp.items():\n",
    "        snpslocuscol = 'unstitched_locus' if 'unstitched_locus' in df.columns else 'locus'\n",
    "        megalocuscol = 'unstitched_locus' if 'unstitched_locus' in megasnps[sp].columns else 'locus'\n",
    "        inter = set(df[snpslocuscol]).intersection(megasnps[sp][megalocuscol])\n",
    "        print(method, sp, len(inter))\n",
    "# # original submission\n",
    "# gatk JP 7408\n",
    "# gatk DF 293\n",
    "# varscan JP 25500\n",
    "# varscan DF 825"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# load trimming and mapping rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evaluate_readinfo(df):\n",
    "    \"\"\"Print out some stats from the readinfo.txt file.\"\"\"\n",
    "    # q30\n",
    "    print(ColorText('\\n\\tBase Quality').bold())\n",
    "    totalbefore = np.nansum(df['total_bases-before_trimming'])\n",
    "    totalafter = np.nansum(df['total_bases-after_trimming'])\n",
    "    q30before = np.nansum(df['q30_bases-before_trimming'])\n",
    "    q30after = np.nansum(df['q30_bases-after_trimming'])\n",
    "    percbefore = '\\t\\tperc before = %s%%' % (round((q30before/totalbefore)*100,2))\n",
    "    percafter = '\\t\\tperc after = %s%%' % (round((q30after/totalafter)*100,2))\n",
    "    print(percbefore)\n",
    "    print(percafter)\n",
    "\n",
    "    # trimming extent\n",
    "    print(ColorText('\\n\\tTrimming').bold())\n",
    "    before = np.nansum(df['total_reads-before_trimming'])\n",
    "    after = np.nansum(df['total_reads-after_trimming'])\n",
    "    print('\\t\\tbefore = %s' % before)\n",
    "    perc = round(after/before, 4)\n",
    "    print('\\t\\tperc = %s%%' % (round((after/before)*100,2)))\n",
    "\n",
    "    # mapping\n",
    "    print(ColorText('\\n\\tMapping').bold())\n",
    "    try:\n",
    "        mapped = np.nansum(df['mapped_bamfile'])\n",
    "        percmapped = '\\t\\tperc mapped = %s%%' % round((mapped/after)*100, 2)\n",
    "        rmdup = np.nansum(df['dedup_bamfile'])\n",
    "        percnodups = '\\t\\tperc of mapped remaining after removing duplicates = %s%%' % round((rmdup/mapped)*100,2)\n",
    "        print(percmapped)\n",
    "        print(percnodups)\n",
    "    except:\n",
    "        pass\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### read in the files, print out data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "40\n",
      "40\n",
      "\u001b[94m\u001b[1m\n",
      "gatk DF\u001b[0m\u001b[0m\n",
      "\u001b[1m\n",
      "\tBase Quality\u001b[0m\n",
      "\t\tperc before = 88.92%\n",
      "\t\tperc after = 90.94%\n",
      "\u001b[1m\n",
      "\tTrimming\u001b[0m\n",
      "\t\tbefore = 419786724.0\n",
      "\t\tperc = 98.07%\n",
      "\u001b[1m\n",
      "\tMapping\u001b[0m\n",
      "\t\tperc mapped = 85.48%\n",
      "\t\tperc of mapped remaining after removing duplicates = 68.68%\n",
      "44\n",
      "40\n",
      "\u001b[94m\u001b[1m\n",
      "gatk JP\u001b[0m\u001b[0m\n",
      "\u001b[1m\n",
      "\tBase Quality\u001b[0m\n",
      "\t\tperc before = 87.28%\n",
      "\t\tperc after = 89.92%\n",
      "\u001b[1m\n",
      "\tTrimming\u001b[0m\n",
      "\t\tbefore = 392241534.0\n",
      "\t\tperc = 97.4%\n",
      "\u001b[1m\n",
      "\tMapping\u001b[0m\n",
      "\t\tperc mapped = 46.14%\n",
      "\t\tperc of mapped remaining after removing duplicates = 79.4%\n",
      "42\n",
      "2\n",
      "\u001b[94m\u001b[1m\n",
      "varscan DF\u001b[0m\u001b[0m\n",
      "\u001b[1m\n",
      "\tBase Quality\u001b[0m\n",
      "\t\tperc before = 88.76%\n",
      "\t\tperc after = 90.44%\n",
      "\u001b[1m\n",
      "\tTrimming\u001b[0m\n",
      "\t\tbefore = 109897366.0\n",
      "\t\tperc = 99.05%\n",
      "\u001b[1m\n",
      "\tMapping\u001b[0m\n",
      "\t\tperc mapped = 86.31%\n",
      "\t\tperc of mapped remaining after removing duplicates = 41.55%\n",
      "118\n",
      "2\n",
      "\u001b[94m\u001b[1m\n",
      "varscan JP\u001b[0m\u001b[0m\n",
      "\u001b[1m\n",
      "\tBase Quality\u001b[0m\n",
      "\t\tperc before = 86.55%\n",
      "\t\tperc after = 89.04%\n",
      "\u001b[1m\n",
      "\tTrimming\u001b[0m\n",
      "\t\tbefore = 152988120.0\n",
      "\t\tperc = 98.06%\n",
      "\u001b[1m\n",
      "\tMapping\u001b[0m\n",
      "7\n",
      "7\n",
      "\u001b[94m\u001b[1m\n",
      "varscan mega DF\u001b[0m\u001b[0m\n",
      "\u001b[1m\n",
      "\tBase Quality\u001b[0m\n",
      "\t\tperc before = 95.63%\n",
      "\t\tperc after = 97.55%\n",
      "\u001b[1m\n",
      "\tTrimming\u001b[0m\n",
      "\t\tbefore = 203296026.0\n",
      "\t\tperc = 97.16%\n",
      "\u001b[1m\n",
      "\tMapping\u001b[0m\n",
      "\t\tperc mapped = 83.52%\n",
      "\t\tperc of mapped remaining after removing duplicates = 99.65%\n",
      "6\n",
      "2\n",
      "\u001b[94m\u001b[1m\n",
      "varscan mega JP\u001b[0m\u001b[0m\n",
      "\u001b[1m\n",
      "\tBase Quality\u001b[0m\n",
      "\t\tperc before = 74.82%\n",
      "\t\tperc after = 78.69%\n",
      "\u001b[1m\n",
      "\tTrimming\u001b[0m\n",
      "\t\tbefore = 201317556.0\n",
      "\t\tperc = 96.49%\n",
      "\u001b[1m\n",
      "\tMapping\u001b[0m\n",
      "\t\tperc mapped = 12.72%\n",
      "\t\tperc of mapped remaining after removing duplicates = 36.7%\n"
     ]
    }
   ],
   "source": [
    "# in some instances, the table I'm loading has other info I don't need here, use convert to reduce df\n",
    "convert = {'gatk JP': 'JP_101',\n",
    "           'gatk DF': 'DF_52',\n",
    "           'varscan JP': 'JP_p101',\n",
    "           'varscan DF': 'DF_p52',\n",
    "           'varscan mega JP': 'JP_RFmg7'}\n",
    "files = ['/data/projects/pool_seq/non-pangenome/gatk_diploid_testdata/DF_i52-gatk/readinfo.txt',\n",
    "         '/data/projects/pool_seq/non-pangenome/gatk_diploid_testdata/JP_i101-gatk/readinfo.txt',\n",
    "         '/data/projects/pool_seq/non-pangenome/varscan_pooled/DF_p52/readinfo.txt',\n",
    "         '/data/projects/pool_seq/non-pangenome/varscan_pooled/JP_pooled/readinfo.txt',\n",
    "         '/data/projects/pool_seq/non-pangenome/varscan_mega/DF_megaSNPs_from_download/DF_mega/readinfo.txt',\n",
    "         '/data/projects/pool_seq/non-pangenome/varscan_mega/JP_RFmg7/readinfo.txt']\n",
    "# read in the files, reduce, and print sequencing statistics\n",
    "dfs = {}\n",
    "for f in files:\n",
    "    df = pd.read_table(f)\n",
    "    print(nrow(df))\n",
    "    method = 'gatk' if 'gatk' in f else 'varscan'\n",
    "    if 'mega' in f:\n",
    "        method = method + ' mega'\n",
    "    sp = 'JP' if 'JP_' in f else 'DF'\n",
    "    msp = method + ' ' + sp\n",
    "    if msp in convert:\n",
    "        df = df[df['samp'].str.contains(convert[msp], regex=False)]\n",
    "    dfs[msp] = df\n",
    "    print(nrow(df))\n",
    "    print(ColorText(f'\\n{method} {sp}').bold().blue())\n",
    "    evaluate_readinfo(df)\n",
    "#     display(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}