WEGS_vs_30XWGS_comparison.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6ad7a48e-643f-44d4-8cf1-d10478e724ed",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ignoring pip: markers 'python_version < \"3\"' don't match your environment\n",
      "Looking in links: /cvmfs/soft.computecanada.ca/custom/python/wheelhouse/gentoo/avx2, /cvmfs/soft.computecanada.ca/custom/python/wheelhouse/gentoo/generic, /cvmfs/soft.computecanada.ca/custom/python/wheelhouse/generic\n",
      "Requirement already satisfied: pysam in /localscratch/dtaliun.15230538.0/lib/python3.9/site-packages (0.20.0+computecanada)\n",
      "Ignoring pip: markers 'python_version < \"3\"' don't match your environment\n",
      "Looking in links: /cvmfs/soft.computecanada.ca/custom/python/wheelhouse/gentoo/avx2, /cvmfs/soft.computecanada.ca/custom/python/wheelhouse/gentoo/generic, /cvmfs/soft.computecanada.ca/custom/python/wheelhouse/generic\n",
      "Requirement already satisfied: intervaltree in /localscratch/dtaliun.15230538.0/lib/python3.9/site-packages (3.1.0+computecanada)\n",
      "Requirement already satisfied: sortedcontainers<3.0,>=2.0 in /localscratch/dtaliun.15230538.0/lib/python3.9/site-packages (from intervaltree) (2.4.0+computecanada)\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "!{sys.executable} -m pip install pysam --no-index\n",
    "!{sys.executable} -m pip install intervaltree --no-index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c9779cd4-0df8-4d10-9219-d3a2ef4a7629",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from intervaltree import Interval, IntervalTree\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "pd.set_option('display.max_columns', None)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "176357e9-42e4-4044-b07b-729824ea1595",
   "metadata": {},
   "source": [
    "## 1. Help functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f02dcc7e-aae2-4116-b8c6-f369e3eeaf8b",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "HAPPY_COLUMNS = ['Type', 'Filter', 'TRUTH.TP', 'QUERY.FP', 'TRUTH.FN', 'METRIC.Recall', 'METRIC.Precision', 'METRIC.F1_Score']\n",
    "\n",
    "def load_targets_bed(filename):\n",
    "    targets = dict()\n",
    "    with open(filename, 'rt') as ifile:\n",
    "        for line in ifile:\n",
    "            chrom, start, stop = line.rstrip().split()\n",
    "            targets.setdefault(chrom, IntervalTree()).addi(int(start) + 1, int(stop) + 1) # because BED is 0-based and IntervalTree stores half-open intervals.\n",
    "    return targets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e24375c3-eb61-4cfa-9794-c62eba66a737",
   "metadata": {},
   "source": [
    "## 2. Input files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "746aae87-a4cd-4472-8851-b544b606a4be",
   "metadata": {},
   "outputs": [],
   "source": [
    "SAMPLE_ANNOTATIONS_FILE = 'sample_annotations.txt'\n",
    "WEGS_TARGETS_BED_FILE = 'WES_targets/SureSelectHumanAllExonV7.b37.Target.bed'\n",
    "NOPLEXWES_BENCHMARK_DIR = 'Variant_calling_analysis/WES/1plex'\n",
    "WEGS_BENCHMARK_DIR = 'Variant_calling_analysis/WEGS'\n",
    "WGS30X_BENCHMARK_DIR = 'WEGS_vs_30XWGS_comparison'\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f90cb8f-017c-4218-8dad-9626c5d2dff8",
   "metadata": {},
   "source": [
    "## 3. Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "12504829-12e7-417d-9319-1da089cc2398",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Load sample annotations\n",
    "df_samples = pd.read_csv(SAMPLE_ANNOTATIONS_FILE, header = 0, sep = '\\t')\n",
    "\n",
    "hgid2naid = {\n",
    "    'HG002': 'NA24385',\n",
    "    'HG003': 'NA24149',\n",
    "    'HG004': 'NA24143'\n",
    "}\n",
    "\n",
    "hgid2rel = {\n",
    "    'HG002': 'Son',\n",
    "    'HG003': 'Father',\n",
    "    'HG004': 'Mother'\n",
    "}\n",
    "\n",
    "WEGS_TARGETS = load_targets_bed(WEGS_TARGETS_BED_FILE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60a265df-360b-4f69-ad51-67602abff072",
   "metadata": {},
   "source": [
    "### 3.1 Load variant calls inside targets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "13992abf-99a7-411d-9835-6f32e0f86dff",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_noplexwes_target_variants = []\n",
    "for row_index, row in df_samples.iterrows():\n",
    "    if row.Plexing != 1:\n",
    "        continue\n",
    "    filepath = f'{NOPLEXWES_BENCHMARK_DIR}/{row.Filename}.norm.hard_filter.happy_benchmark.summary.csv'\n",
    "    df = pd.read_csv(filepath, usecols = HAPPY_COLUMNS).rename(columns = {\n",
    "            'METRIC.Recall': 'Recall', \n",
    "            'METRIC.Precision': 'Precision',\n",
    "            'METRIC.F1_Score': 'F1_Score',\n",
    "            'TRUTH.TP': 'TP',\n",
    "            'QUERY.FP': 'FP',\n",
    "            'TRUTH.FN': 'FN'\n",
    "        })\n",
    "    df['Total_called'] = df['TP'] + df['FP']\n",
    "    df['ID'] = row.Filename\n",
    "    df['HGID'] = row.HGID\n",
    "    df['Batch'] = row.Batch\n",
    "    df['Plexing'] = row.Plexing\n",
    "    df_noplexwes_target_variants.append(df)\n",
    "df_noplexwes_target_variants = pd.concat(df_noplexwes_target_variants).reset_index(drop = True)\n",
    "\n",
    "\n",
    "df_wegs_target_variants = []\n",
    "for row_index, row in df_samples.iterrows():\n",
    "    if row.Plexing == 1:\n",
    "        continue\n",
    "    for wgs_depth in [2, 5]:\n",
    "        filepath = f'{WEGS_BENCHMARK_DIR}/{row.Plexing}plex_{wgs_depth}X/{row.Filename}.norm.hard_filter.happy_benchmark.summary.csv'\n",
    "        df = pd.read_csv(filepath, usecols = HAPPY_COLUMNS).rename(columns = {\n",
    "            'METRIC.Recall': 'Recall', \n",
    "            'METRIC.Precision': 'Precision',\n",
    "            'METRIC.F1_Score': 'F1_Score',\n",
    "            'TRUTH.TP': 'TP',\n",
    "            'QUERY.FP': 'FP',\n",
    "            'TRUTH.FN': 'FN'\n",
    "        })\n",
    "        df['Total_called'] = df['TP'] + df['FP']\n",
    "        df['ID'] = row.Filename\n",
    "        df['HGID'] = row.HGID\n",
    "        df['Batch'] = row.Batch\n",
    "        df['Plexing'] = row.Plexing\n",
    "        df['WGS_depth'] = wgs_depth\n",
    "        df_wegs_target_variants.append(df)\n",
    "df_wegs_target_variants = pd.concat(df_wegs_target_variants).reset_index(drop = True)\n",
    "\n",
    "df_wgs30x_target_variants = []\n",
    "for row_index, row in df_samples[['HGID', 'Family']].drop_duplicates().iterrows():\n",
    "    for i in range(1, 6):\n",
    "        filepath = f'{WGS30X_BENCHMARK_DIR}/Targets/{row.HGID}_{hgid2naid[row.HGID]}_{row.Family}_S{i}.norm.hard_filter.happy_benchmark.summary.csv'\n",
    "        df = pd.read_csv(filepath, usecols = HAPPY_COLUMNS).rename(columns = {\n",
    "            'METRIC.Recall': 'Recall', \n",
    "            'METRIC.Precision': 'Precision',\n",
    "            'METRIC.F1_Score': 'F1_Score',\n",
    "            'TRUTH.TP': 'TP',\n",
    "            'QUERY.FP': 'FP',\n",
    "            'TRUTH.FN': 'FN'\n",
    "        })\n",
    "        df['Total_called'] = df['TP'] + df['FP']\n",
    "        df['HGID'] = row.HGID\n",
    "        df['Downsample'] = i\n",
    "        df_wgs30x_target_variants.append(df)\n",
    "df_wgs30x_target_variants = pd.concat(df_wgs30x_target_variants).reset_index(drop = True)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a77277b5-8be9-43c7-84f6-0e3a6bef3e7f",
   "metadata": {},
   "source": [
    "### 3.2 Load variant calls genome-wide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "060aa0d8-8818-4aaf-bd89-76a6ca8d3b81",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_wegs_variants = []\n",
    "for row_index, row in df_samples.iterrows():\n",
    "    if row.Plexing == 1:\n",
    "        continue\n",
    "    for wgs_depth in [2, 5]:\n",
    "        filepath = f'{WEGS_BENCHMARK_DIR}/{row.Plexing}plex_{wgs_depth}X_genomewide/{row.Filename}.norm.hard_filter.happy_benchmark.summary.csv'\n",
    "        df = pd.read_csv(filepath, usecols = HAPPY_COLUMNS).rename(columns = {\n",
    "            'METRIC.Recall': 'Recall', \n",
    "            'METRIC.Precision': 'Precision',\n",
    "            'METRIC.F1_Score': 'F1_Score',\n",
    "            'TRUTH.TP': 'TP',\n",
    "            'QUERY.FP': 'FP',\n",
    "            'TRUTH.FN': 'FN'\n",
    "        })\n",
    "        df['Total_called'] = df['TP'] + df['FP']\n",
    "        df['ID'] = row.Filename\n",
    "        df['HGID'] = row.HGID\n",
    "        df['Batch'] = row.Batch\n",
    "        df['Plexing'] = row.Plexing\n",
    "        df['WGS_depth'] = wgs_depth\n",
    "        df_wegs_variants.append(df)\n",
    "df_wegs_variants = pd.concat(df_wegs_variants).reset_index(drop = True)\n",
    "\n",
    "df_wgs30x_variants = []\n",
    "for row_index, row in df_samples[['HGID', 'Family']].drop_duplicates().iterrows():\n",
    "    for i in range(1, 6):\n",
    "        filepath = f'{WGS30X_BENCHMARK_DIR}/Genomewide/{row.HGID}_{hgid2naid[row.HGID]}_{row.Family}_S{i}.norm.hard_filter.happy_benchmark.summary.csv'\n",
    "        df = pd.read_csv(filepath, usecols = HAPPY_COLUMNS).rename(columns = {\n",
    "            'METRIC.Recall': 'Recall', \n",
    "            'METRIC.Precision': 'Precision',\n",
    "            'METRIC.F1_Score': 'F1_Score',\n",
    "            'TRUTH.TP': 'TP',\n",
    "            'QUERY.FP': 'FP',\n",
    "            'TRUTH.FN': 'FN'\n",
    "        })\n",
    "        df['Total_called'] = df['TP'] + df['FP']\n",
    "        df['HGID'] = row.HGID\n",
    "        df['Downsample'] = i\n",
    "        df_wgs30x_variants.append(df)\n",
    "df_wgs30x_variants = pd.concat(df_wgs30x_variants).reset_index(drop = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b5d8355-0334-4c2a-9ef0-6a686e0cbbd2",
   "metadata": {},
   "source": [
    "## 4. Figures and tables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f3bb49a-3eec-4357-abc6-1672c2f3e07c",
   "metadata": {},
   "source": [
    "### TABLE: Average variant recall and precision rates in 30X WGS, WES, and WEGS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5ec0d420-b128-4045-923e-d2f4e55db876",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Label</th>\n",
       "      <th>SNV TP</th>\n",
       "      <th>SNV FP</th>\n",
       "      <th>SNV FN</th>\n",
       "      <th>SNV Precision</th>\n",
       "      <th>SNV Recall</th>\n",
       "      <th>InDel TP</th>\n",
       "      <th>InDel FP</th>\n",
       "      <th>InDel FN</th>\n",
       "      <th>InDel Precision</th>\n",
       "      <th>InDel Recall</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>WGS 30X</td>\n",
       "      <td>22,338 (19)</td>\n",
       "      <td>260 (5)</td>\n",
       "      <td>287 (5)</td>\n",
       "      <td>0.9885 (0.0002)</td>\n",
       "      <td>0.9873 (0.0002)</td>\n",
       "      <td>661 (5)</td>\n",
       "      <td>12 (1)</td>\n",
       "      <td>11 (1)</td>\n",
       "      <td>0.9823 (0.0017)</td>\n",
       "      <td>0.9841 (0.0011)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>WES 100X</td>\n",
       "      <td>22,241 (30)</td>\n",
       "      <td>406 (7)</td>\n",
       "      <td>384 (8)</td>\n",
       "      <td>0.9821 (0.0003)</td>\n",
       "      <td>0.9830 (0.0004)</td>\n",
       "      <td>631 (6)</td>\n",
       "      <td>53 (2)</td>\n",
       "      <td>41 (2)</td>\n",
       "      <td>0.9232 (0.0028)</td>\n",
       "      <td>0.9390 (0.0029)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>WEGS 4P2X</td>\n",
       "      <td>22,268 (24)</td>\n",
       "      <td>426 (6)</td>\n",
       "      <td>357 (5)</td>\n",
       "      <td>0.9812 (0.0002)</td>\n",
       "      <td>0.9842 (0.0002)</td>\n",
       "      <td>638 (5)</td>\n",
       "      <td>51 (2)</td>\n",
       "      <td>34 (2)</td>\n",
       "      <td>0.9269 (0.0024)</td>\n",
       "      <td>0.9493 (0.0028)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>WEGS 4P5X</td>\n",
       "      <td>22,291 (25)</td>\n",
       "      <td>427 (6)</td>\n",
       "      <td>334 (3)</td>\n",
       "      <td>0.9812 (0.0002)</td>\n",
       "      <td>0.9852 (0.0001)</td>\n",
       "      <td>642 (5)</td>\n",
       "      <td>49 (2)</td>\n",
       "      <td>30 (1)</td>\n",
       "      <td>0.9300 (0.0024)</td>\n",
       "      <td>0.9552 (0.0019)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>WEGS 8P2X</td>\n",
       "      <td>22,247 (21)</td>\n",
       "      <td>421 (6)</td>\n",
       "      <td>376 (4)</td>\n",
       "      <td>0.9814 (0.0003)</td>\n",
       "      <td>0.9834 (0.0002)</td>\n",
       "      <td>633 (4)</td>\n",
       "      <td>52 (2)</td>\n",
       "      <td>40 (2)</td>\n",
       "      <td>0.9240 (0.0024)</td>\n",
       "      <td>0.9407 (0.0026)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>WEGS 8P5X</td>\n",
       "      <td>22,277 (21)</td>\n",
       "      <td>418 (6)</td>\n",
       "      <td>346 (3)</td>\n",
       "      <td>0.9816 (0.0002)</td>\n",
       "      <td>0.9847 (0.0001)</td>\n",
       "      <td>638 (4)</td>\n",
       "      <td>50 (2)</td>\n",
       "      <td>34 (2)</td>\n",
       "      <td>0.9277 (0.0027)</td>\n",
       "      <td>0.9490 (0.0020)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Label       SNV TP   SNV FP   SNV FN    SNV Precision       SNV Recall  \\\n",
       "0    WGS 30X  22,338 (19)  260 (5)  287 (5)  0.9885 (0.0002)  0.9873 (0.0002)   \n",
       "1   WES 100X  22,241 (30)  406 (7)  384 (8)  0.9821 (0.0003)  0.9830 (0.0004)   \n",
       "2  WEGS 4P2X  22,268 (24)  426 (6)  357 (5)  0.9812 (0.0002)  0.9842 (0.0002)   \n",
       "3  WEGS 4P5X  22,291 (25)  427 (6)  334 (3)  0.9812 (0.0002)  0.9852 (0.0001)   \n",
       "4  WEGS 8P2X  22,247 (21)  421 (6)  376 (4)  0.9814 (0.0003)  0.9834 (0.0002)   \n",
       "5  WEGS 8P5X  22,277 (21)  418 (6)  346 (3)  0.9816 (0.0002)  0.9847 (0.0001)   \n",
       "\n",
       "  InDel TP InDel FP InDel FN  InDel Precision     InDel Recall  \n",
       "0  661 (5)   12 (1)   11 (1)  0.9823 (0.0017)  0.9841 (0.0011)  \n",
       "1  631 (6)   53 (2)   41 (2)  0.9232 (0.0028)  0.9390 (0.0029)  \n",
       "2  638 (5)   51 (2)   34 (2)  0.9269 (0.0024)  0.9493 (0.0028)  \n",
       "3  642 (5)   49 (2)   30 (1)  0.9300 (0.0024)  0.9552 (0.0019)  \n",
       "4  633 (4)   52 (2)   40 (2)  0.9240 (0.0024)  0.9407 (0.0026)  \n",
       "5  638 (4)   50 (2)   34 (2)  0.9277 (0.0027)  0.9490 (0.0020)  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dfs = [\n",
    "    ( 'WGS 30X', df_wgs30x_target_variants), \n",
    "    ( 'WES 100X', df_noplexwes_target_variants),\n",
    "    ( 'WEGS 4P2X', df_wegs_target_variants[(df_wegs_target_variants.Plexing == 4) & (df_wegs_target_variants.WGS_depth == 2)]),\n",
    "    ( 'WEGS 4P5X', df_wegs_target_variants[(df_wegs_target_variants.Plexing == 4) & (df_wegs_target_variants.WGS_depth == 5)]),\n",
    "    ( 'WEGS 8P2X', df_wegs_target_variants[(df_wegs_target_variants.Plexing == 8) & (df_wegs_target_variants.WGS_depth == 2)]),\n",
    "    ( 'WEGS 8P5X', df_wegs_target_variants[(df_wegs_target_variants.Plexing == 8) & (df_wegs_target_variants.WGS_depth == 5)])\n",
    "]\n",
    "\n",
    "data = {\n",
    "    'Label': [],\n",
    "    'SNV TP': [],\n",
    "    'SNV FP': [],\n",
    "    'SNV FN': [],\n",
    "    'SNV Precision': [],\n",
    "    'SNV Recall': [],\n",
    "    'InDel TP': [],\n",
    "    'InDel FP': [],\n",
    "    'InDel FN': [],\n",
    "    'InDel Precision': [],\n",
    "    'InDel Recall': []\n",
    "}\n",
    "\n",
    "for label, df in dfs:\n",
    "    data['Label'].append(label)\n",
    "    \n",
    "    df_temp = df[(df.Type == 'SNP') & (df.Filter == 'ALL')]\n",
    "    data['SNV TP'].append(f'{np.mean(df_temp.TP):,.0f} ({stats.sem(df_temp.TP):,.0f})')\n",
    "    data['SNV FP'].append(f'{np.mean(df_temp.FP):,.0f} ({stats.sem(df_temp.FP):,.0f})')\n",
    "    data['SNV FN'].append(f'{np.mean(df_temp.FN):,.0f} ({stats.sem(df_temp.FN):,.0f})')\n",
    "    data['SNV Precision'].append(f'{np.mean(df_temp.Precision):.4f} ({stats.sem(df_temp.Precision):,.4f})')\n",
    "    data['SNV Recall'].append(f'{np.mean(df_temp.Recall):.4f} ({stats.sem(df_temp.Recall):,.4f})')\n",
    "\n",
    "    df_temp = df[(df.Type == 'INDEL') & (df.Filter == 'ALL')]\n",
    "    data['InDel TP'].append(f'{np.mean(df_temp.TP):,.0f} ({stats.sem(df_temp.TP):,.0f})')\n",
    "    data['InDel FP'].append(f'{np.mean(df_temp.FP):,.0f} ({stats.sem(df_temp.FP):,.0f})')\n",
    "    data['InDel FN'].append(f'{np.mean(df_temp.FN):,.0f} ({stats.sem(df_temp.FN):,.0f})')\n",
    "    data['InDel Precision'].append(f'{np.mean(df_temp.Precision):.4f} ({stats.sem(df_temp.Precision):,.4f})')\n",
    "    data['InDel Recall'].append(f'{np.mean(df_temp.Recall):.4f} ({stats.sem(df_temp.Recall):,.4f})')\n",
    "\n",
    "pd.DataFrame(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc0ba5a2-5bf1-4741-bd6b-a7dd80782548",
   "metadata": {},
   "source": [
    "### EXTRA TABLE: Average SNV recall and precision rates in 30X WGS, WES, and WEGS stratified by individual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "386d1a05-90ce-44fb-bd08-8d7e6a291fe9",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Sample</th>\n",
       "      <th>WGS 30X TP</th>\n",
       "      <th>WGS 30X FN</th>\n",
       "      <th>WGS 30X FP</th>\n",
       "      <th>WGS 30X Precision</th>\n",
       "      <th>WGS 30X Recall</th>\n",
       "      <th>WES 100X TP</th>\n",
       "      <th>WES 100X FN</th>\n",
       "      <th>WES 100X FP</th>\n",
       "      <th>WES 100X Precision</th>\n",
       "      <th>WES 100X Recall</th>\n",
       "      <th>WEGS 4P2X TP</th>\n",
       "      <th>WEGS 4P2X FN</th>\n",
       "      <th>WEGS 4P2X FP</th>\n",
       "      <th>WEGS 4P2X Precision</th>\n",
       "      <th>WEGS 4P2X Recall</th>\n",
       "      <th>WEGS 8P5X TP</th>\n",
       "      <th>WEGS 8P5X FN</th>\n",
       "      <th>WEGS 8P5X FP</th>\n",
       "      <th>WEGS 8P5X Precision</th>\n",
       "      <th>WEGS 8P5X Recall</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>HG002</td>\n",
       "      <td>22,433 (4)</td>\n",
       "      <td>310 (4)</td>\n",
       "      <td>284 (5)</td>\n",
       "      <td>0.9875 (0.0002)</td>\n",
       "      <td>0.9864 (0.0002)</td>\n",
       "      <td>22,356 (3)</td>\n",
       "      <td>387 (3)</td>\n",
       "      <td>408 (3)</td>\n",
       "      <td>0.9821 (0.0001)</td>\n",
       "      <td>0.9830 (0.0001)</td>\n",
       "      <td>22,379 (7)</td>\n",
       "      <td>364 (7)</td>\n",
       "      <td>446 (6)</td>\n",
       "      <td>0.9805 (0.0002)</td>\n",
       "      <td>0.9840 (0.0003)</td>\n",
       "      <td>22,390 (5)</td>\n",
       "      <td>353 (5)</td>\n",
       "      <td>449 (7)</td>\n",
       "      <td>0.9803 (0.0003)</td>\n",
       "      <td>0.9845 (0.0002)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>HG003</td>\n",
       "      <td>22,312 (4)</td>\n",
       "      <td>277 (4)</td>\n",
       "      <td>245 (3)</td>\n",
       "      <td>0.9892 (0.0001)</td>\n",
       "      <td>0.9877 (0.0002)</td>\n",
       "      <td>22,202 (17)</td>\n",
       "      <td>387 (17)</td>\n",
       "      <td>405 (6)</td>\n",
       "      <td>0.9821 (0.0003)</td>\n",
       "      <td>0.9829 (0.0008)</td>\n",
       "      <td>22,231 (9)</td>\n",
       "      <td>358 (9)</td>\n",
       "      <td>425 (4)</td>\n",
       "      <td>0.9813 (0.0001)</td>\n",
       "      <td>0.9841 (0.0004)</td>\n",
       "      <td>22,249 (6)</td>\n",
       "      <td>340 (6)</td>\n",
       "      <td>406 (5)</td>\n",
       "      <td>0.9821 (0.0002)</td>\n",
       "      <td>0.9850 (0.0003)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>HG004</td>\n",
       "      <td>22,269 (5)</td>\n",
       "      <td>274 (5)</td>\n",
       "      <td>251 (6)</td>\n",
       "      <td>0.9888 (0.0003)</td>\n",
       "      <td>0.9878 (0.0002)</td>\n",
       "      <td>22,167 (22)</td>\n",
       "      <td>376 (22)</td>\n",
       "      <td>406 (24)</td>\n",
       "      <td>0.9820 (0.0010)</td>\n",
       "      <td>0.9833 (0.0010)</td>\n",
       "      <td>22,195 (9)</td>\n",
       "      <td>348 (9)</td>\n",
       "      <td>406 (7)</td>\n",
       "      <td>0.9820 (0.0003)</td>\n",
       "      <td>0.9846 (0.0004)</td>\n",
       "      <td>22,197 (5)</td>\n",
       "      <td>346 (5)</td>\n",
       "      <td>402 (3)</td>\n",
       "      <td>0.9822 (0.0001)</td>\n",
       "      <td>0.9847 (0.0002)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Sample  WGS 30X TP WGS 30X FN WGS 30X FP WGS 30X Precision   WGS 30X Recall  \\\n",
       "0  HG002  22,433 (4)    310 (4)    284 (5)   0.9875 (0.0002)  0.9864 (0.0002)   \n",
       "1  HG003  22,312 (4)    277 (4)    245 (3)   0.9892 (0.0001)  0.9877 (0.0002)   \n",
       "2  HG004  22,269 (5)    274 (5)    251 (6)   0.9888 (0.0003)  0.9878 (0.0002)   \n",
       "\n",
       "   WES 100X TP WES 100X FN WES 100X FP WES 100X Precision  WES 100X Recall  \\\n",
       "0   22,356 (3)     387 (3)     408 (3)    0.9821 (0.0001)  0.9830 (0.0001)   \n",
       "1  22,202 (17)    387 (17)     405 (6)    0.9821 (0.0003)  0.9829 (0.0008)   \n",
       "2  22,167 (22)    376 (22)    406 (24)    0.9820 (0.0010)  0.9833 (0.0010)   \n",
       "\n",
       "  WEGS 4P2X TP WEGS 4P2X FN WEGS 4P2X FP WEGS 4P2X Precision WEGS 4P2X Recall  \\\n",
       "0   22,379 (7)      364 (7)      446 (6)     0.9805 (0.0002)  0.9840 (0.0003)   \n",
       "1   22,231 (9)      358 (9)      425 (4)     0.9813 (0.0001)  0.9841 (0.0004)   \n",
       "2   22,195 (9)      348 (9)      406 (7)     0.9820 (0.0003)  0.9846 (0.0004)   \n",
       "\n",
       "  WEGS 8P5X TP WEGS 8P5X FN WEGS 8P5X FP WEGS 8P5X Precision WEGS 8P5X Recall  \n",
       "0   22,390 (5)      353 (5)      449 (7)     0.9803 (0.0003)  0.9845 (0.0002)  \n",
       "1   22,249 (6)      340 (6)      406 (5)     0.9821 (0.0002)  0.9850 (0.0003)  \n",
       "2   22,197 (5)      346 (5)      402 (3)     0.9822 (0.0001)  0.9847 (0.0002)  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = {\n",
    "    'Sample': [],\n",
    "    'WGS 30X TP': [],\n",
    "    'WGS 30X FN': [],\n",
    "    'WGS 30X FP': [],\n",
    "    'WGS 30X Precision': [],\n",
    "    'WGS 30X Recall': [],\n",
    "    'WES 100X TP': [],\n",
    "    'WES 100X FN': [],\n",
    "    'WES 100X FP': [],\n",
    "    'WES 100X Precision': [],\n",
    "    'WES 100X Recall': [],\n",
    "    'WEGS 4P2X TP': [],\n",
    "    'WEGS 4P2X FN': [],\n",
    "    'WEGS 4P2X FP': [],\n",
    "    'WEGS 4P2X Precision': [],\n",
    "    'WEGS 4P2X Recall': [],\n",
    "    'WEGS 8P5X TP': [],\n",
    "    'WEGS 8P5X FN': [],\n",
    "    'WEGS 8P5X FP': [],\n",
    "    'WEGS 8P5X Precision': [],\n",
    "    'WEGS 8P5X Recall': [],\n",
    "}\n",
    "\n",
    "for hgid, naid in hgid2naid.items():\n",
    "    data['Sample'].append(hgid)\n",
    "    \n",
    "    df_temp = df_wgs30x_target_variants[df_wgs30x_target_variants.HGID == hgid]\n",
    "    df_temp = df_temp[(df_temp.Type == 'SNP') & (df_temp.Filter == 'ALL')][['TP', 'FN', 'FP', 'Recall', 'Precision']]\n",
    "    data['WGS 30X TP'].append(f'{np.mean(df_temp.TP):,.0f} ({stats.sem(df_temp.TP):,.0f})')\n",
    "    data['WGS 30X FP'].append(f'{np.mean(df_temp.FP):,.0f} ({stats.sem(df_temp.FP):,.0f})')\n",
    "    data['WGS 30X FN'].append(f'{np.mean(df_temp.FN):,.0f} ({stats.sem(df_temp.FN):,.0f})')\n",
    "    data['WGS 30X Precision'].append(f'{np.mean(df_temp.Precision):.4f} ({stats.sem(df_temp.Precision):,.4f})')\n",
    "    data['WGS 30X Recall'].append(f'{np.mean(df_temp.Recall):.4f} ({stats.sem(df_temp.Recall):,.4f})')\n",
    "\n",
    "    df_temp = df_noplexwes_target_variants[df_noplexwes_target_variants.HGID == hgid]\n",
    "    df_temp = df_temp[(df_temp.Type == 'SNP') & (df_temp.Filter == 'ALL')][['TP', 'FN', 'FP', 'Recall', 'Precision']]\n",
    "    data['WES 100X TP'].append(f'{np.mean(df_temp.TP):,.0f} ({stats.sem(df_temp.TP):,.0f})')\n",
    "    data['WES 100X FP'].append(f'{np.mean(df_temp.FP):,.0f} ({stats.sem(df_temp.FP):,.0f})')\n",
    "    data['WES 100X FN'].append(f'{np.mean(df_temp.FN):,.0f} ({stats.sem(df_temp.FN):,.0f})')\n",
    "    data['WES 100X Precision'].append(f'{np.mean(df_temp.Precision):.4f} ({stats.sem(df_temp.Precision):,.4f})')\n",
    "    data['WES 100X Recall'].append(f'{np.mean(df_temp.Recall):.4f} ({stats.sem(df_temp.Recall):,.4f})')\n",
    "\n",
    "    df_temp = df_wegs_target_variants[df_wegs_target_variants.HGID == hgid]\n",
    "    df_temp = df_temp[(df_temp.Type == 'SNP') & (df_temp.Filter == 'ALL') & (df_temp.Plexing == 4) & (df_temp.WGS_depth == 2)]\n",
    "    data['WEGS 4P2X TP'].append(f'{np.mean(df_temp.TP):,.0f} ({stats.sem(df_temp.TP):,.0f})')\n",
    "    data['WEGS 4P2X FP'].append(f'{np.mean(df_temp.FP):,.0f} ({stats.sem(df_temp.FP):,.0f})')\n",
    "    data['WEGS 4P2X FN'].append(f'{np.mean(df_temp.FN):,.0f} ({stats.sem(df_temp.FN):,.0f})')\n",
    "    data['WEGS 4P2X Precision'].append(f'{np.mean(df_temp.Precision):.4f} ({stats.sem(df_temp.Precision):,.4f})')\n",
    "    data['WEGS 4P2X Recall'].append(f'{np.mean(df_temp.Recall):.4f} ({stats.sem(df_temp.Recall):,.4f})')\n",
    "    \n",
    "    df_temp = df_wegs_target_variants[df_wegs_target_variants.HGID == hgid]\n",
    "    df_temp = df_temp[(df_temp.Type == 'SNP') & (df_temp.Filter == 'ALL') & (df_temp.Plexing == 8) & (df_temp.WGS_depth == 5)]\n",
    "    data['WEGS 8P5X TP'].append(f'{np.mean(df_temp.TP):,.0f} ({stats.sem(df_temp.TP):,.0f})')\n",
    "    data['WEGS 8P5X FP'].append(f'{np.mean(df_temp.FP):,.0f} ({stats.sem(df_temp.FP):,.0f})')\n",
    "    data['WEGS 8P5X FN'].append(f'{np.mean(df_temp.FN):,.0f} ({stats.sem(df_temp.FN):,.0f})')\n",
    "    data['WEGS 8P5X Precision'].append(f'{np.mean(df_temp.Precision):.4f} ({stats.sem(df_temp.Precision):,.4f})')\n",
    "    data['WEGS 8P5X Recall'].append(f'{np.mean(df_temp.Recall):.4f} ({stats.sem(df_temp.Recall):,.4f})')\n",
    "    \n",
    "\n",
    "pd.DataFrame(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9499df8f-832b-40fd-b66a-a0daa08f499c",
   "metadata": {},
   "source": [
    "### EXTRA TABLE: Average InDel recall and precision rates in 30X WGS, WES, and WEGS stratified by individual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "cb872cd7-1d90-430d-b821-cc17b7386197",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Sample</th>\n",
       "      <th>WGS 30X TP</th>\n",
       "      <th>WGS 30X FN</th>\n",
       "      <th>WGS 30X FP</th>\n",
       "      <th>WGS 30X Precision</th>\n",
       "      <th>WGS 30X Recall</th>\n",
       "      <th>WES 100X TP</th>\n",
       "      <th>WES 100X FN</th>\n",
       "      <th>WES 100X FP</th>\n",
       "      <th>WES 100X Precision</th>\n",
       "      <th>WES 100X Recall</th>\n",
       "      <th>WEGS 4P2X TP</th>\n",
       "      <th>WEGS 4P2X FN</th>\n",
       "      <th>WEGS 4P2X FP</th>\n",
       "      <th>WEGS 4P2X Precision</th>\n",
       "      <th>WEGS 4P2X Recall</th>\n",
       "      <th>WEGS 8P5X TP</th>\n",
       "      <th>WEGS 8P5X FN</th>\n",
       "      <th>WEGS 8P5X FP</th>\n",
       "      <th>WEGS 8P5X Precision</th>\n",
       "      <th>WEGS 8P5X Recall</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>HG002</td>\n",
       "      <td>673 (1)</td>\n",
       "      <td>9 (1)</td>\n",
       "      <td>16 (1)</td>\n",
       "      <td>0.9765 (0.0009)</td>\n",
       "      <td>0.9865 (0.0009)</td>\n",
       "      <td>641 (1)</td>\n",
       "      <td>41 (1)</td>\n",
       "      <td>60 (1)</td>\n",
       "      <td>0.9147 (0.0014)</td>\n",
       "      <td>0.9399 (0.0015)</td>\n",
       "      <td>648 (3)</td>\n",
       "      <td>34 (3)</td>\n",
       "      <td>57 (1)</td>\n",
       "      <td>0.9197 (0.0017)</td>\n",
       "      <td>0.9498 (0.0050)</td>\n",
       "      <td>647 (2)</td>\n",
       "      <td>35 (2)</td>\n",
       "      <td>58 (3)</td>\n",
       "      <td>0.9177 (0.0038)</td>\n",
       "      <td>0.9481 (0.0028)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>HG003</td>\n",
       "      <td>676 (1)</td>\n",
       "      <td>14 (1)</td>\n",
       "      <td>13 (2)</td>\n",
       "      <td>0.9809 (0.0028)</td>\n",
       "      <td>0.9797 (0.0013)</td>\n",
       "      <td>642 (3)</td>\n",
       "      <td>48 (3)</td>\n",
       "      <td>50 (3)</td>\n",
       "      <td>0.9280 (0.0036)</td>\n",
       "      <td>0.9309 (0.0040)</td>\n",
       "      <td>649 (2)</td>\n",
       "      <td>41 (2)</td>\n",
       "      <td>52 (2)</td>\n",
       "      <td>0.9254 (0.0019)</td>\n",
       "      <td>0.9409 (0.0027)</td>\n",
       "      <td>651 (1)</td>\n",
       "      <td>39 (1)</td>\n",
       "      <td>52 (1)</td>\n",
       "      <td>0.9265 (0.0018)</td>\n",
       "      <td>0.9432 (0.0019)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>HG004</td>\n",
       "      <td>634 (1)</td>\n",
       "      <td>9 (1)</td>\n",
       "      <td>7 (1)</td>\n",
       "      <td>0.9895 (0.0009)</td>\n",
       "      <td>0.9860 (0.0017)</td>\n",
       "      <td>608 (3)</td>\n",
       "      <td>35 (3)</td>\n",
       "      <td>48 (3)</td>\n",
       "      <td>0.9271 (0.0048)</td>\n",
       "      <td>0.9461 (0.0049)</td>\n",
       "      <td>616 (2)</td>\n",
       "      <td>28 (2)</td>\n",
       "      <td>43 (2)</td>\n",
       "      <td>0.9357 (0.0035)</td>\n",
       "      <td>0.9572 (0.0026)</td>\n",
       "      <td>615 (2)</td>\n",
       "      <td>28 (2)</td>\n",
       "      <td>40 (2)</td>\n",
       "      <td>0.9393 (0.0033)</td>\n",
       "      <td>0.9568 (0.0031)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Sample WGS 30X TP WGS 30X FN WGS 30X FP WGS 30X Precision   WGS 30X Recall  \\\n",
       "0  HG002    673 (1)      9 (1)     16 (1)   0.9765 (0.0009)  0.9865 (0.0009)   \n",
       "1  HG003    676 (1)     14 (1)     13 (2)   0.9809 (0.0028)  0.9797 (0.0013)   \n",
       "2  HG004    634 (1)      9 (1)      7 (1)   0.9895 (0.0009)  0.9860 (0.0017)   \n",
       "\n",
       "  WES 100X TP WES 100X FN WES 100X FP WES 100X Precision  WES 100X Recall  \\\n",
       "0     641 (1)      41 (1)      60 (1)    0.9147 (0.0014)  0.9399 (0.0015)   \n",
       "1     642 (3)      48 (3)      50 (3)    0.9280 (0.0036)  0.9309 (0.0040)   \n",
       "2     608 (3)      35 (3)      48 (3)    0.9271 (0.0048)  0.9461 (0.0049)   \n",
       "\n",
       "  WEGS 4P2X TP WEGS 4P2X FN WEGS 4P2X FP WEGS 4P2X Precision WEGS 4P2X Recall  \\\n",
       "0      648 (3)       34 (3)       57 (1)     0.9197 (0.0017)  0.9498 (0.0050)   \n",
       "1      649 (2)       41 (2)       52 (2)     0.9254 (0.0019)  0.9409 (0.0027)   \n",
       "2      616 (2)       28 (2)       43 (2)     0.9357 (0.0035)  0.9572 (0.0026)   \n",
       "\n",
       "  WEGS 8P5X TP WEGS 8P5X FN WEGS 8P5X FP WEGS 8P5X Precision WEGS 8P5X Recall  \n",
       "0      647 (2)       35 (2)       58 (3)     0.9177 (0.0038)  0.9481 (0.0028)  \n",
       "1      651 (1)       39 (1)       52 (1)     0.9265 (0.0018)  0.9432 (0.0019)  \n",
       "2      615 (2)       28 (2)       40 (2)     0.9393 (0.0033)  0.9568 (0.0031)  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = {\n",
    "    'Sample': [],\n",
    "    'WGS 30X TP': [],\n",
    "    'WGS 30X FN': [],\n",
    "    'WGS 30X FP': [],\n",
    "    'WGS 30X Precision': [],\n",
    "    'WGS 30X Recall': [],\n",
    "    'WES 100X TP': [],\n",
    "    'WES 100X FN': [],\n",
    "    'WES 100X FP': [],\n",
    "    'WES 100X Precision': [],\n",
    "    'WES 100X Recall': [],\n",
    "    'WEGS 4P2X TP': [],\n",
    "    'WEGS 4P2X FN': [],\n",
    "    'WEGS 4P2X FP': [],\n",
    "    'WEGS 4P2X Precision': [],\n",
    "    'WEGS 4P2X Recall': [],\n",
    "    'WEGS 8P5X TP': [],\n",
    "    'WEGS 8P5X FN': [],\n",
    "    'WEGS 8P5X FP': [],\n",
    "    'WEGS 8P5X Precision': [],\n",
    "    'WEGS 8P5X Recall': [],\n",
    "}\n",
    "\n",
    "for hgid, naid in hgid2naid.items():\n",
    "    data['Sample'].append(hgid)\n",
    "\n",
    "    df_temp = df_wgs30x_target_variants[df_wgs30x_target_variants.HGID == hgid]\n",
    "    df_temp = df_temp[(df_temp.Type == 'INDEL') & (df_temp.Filter == 'ALL')][['TP', 'FN', 'FP', 'Recall', 'Precision']]\n",
    "    data['WGS 30X TP'].append(f'{np.mean(df_temp.TP):,.0f} ({stats.sem(df_temp.TP):,.0f})')\n",
    "    data['WGS 30X FP'].append(f'{np.mean(df_temp.FP):,.0f} ({stats.sem(df_temp.FP):,.0f})')\n",
    "    data['WGS 30X FN'].append(f'{np.mean(df_temp.FN):,.0f} ({stats.sem(df_temp.FN):,.0f})')\n",
    "    data['WGS 30X Precision'].append(f'{np.mean(df_temp.Precision):.4f} ({stats.sem(df_temp.Precision):,.4f})')\n",
    "    data['WGS 30X Recall'].append(f'{np.mean(df_temp.Recall):.4f} ({stats.sem(df_temp.Recall):,.4f})')\n",
    "\n",
    "    df_temp = df_noplexwes_target_variants[df_noplexwes_target_variants.HGID == hgid]\n",
    "    df_temp = df_temp[(df_temp.Type == 'INDEL') & (df_temp.Filter == 'ALL')][['TP', 'FN', 'FP', 'Recall', 'Precision']]\n",
    "    data['WES 100X TP'].append(f'{np.mean(df_temp.TP):,.0f} ({stats.sem(df_temp.TP):,.0f})')\n",
    "    data['WES 100X FP'].append(f'{np.mean(df_temp.FP):,.0f} ({stats.sem(df_temp.FP):,.0f})')\n",
    "    data['WES 100X FN'].append(f'{np.mean(df_temp.FN):,.0f} ({stats.sem(df_temp.FN):,.0f})')\n",
    "    data['WES 100X Precision'].append(f'{np.mean(df_temp.Precision):.4f} ({stats.sem(df_temp.Precision):,.4f})')\n",
    "    data['WES 100X Recall'].append(f'{np.mean(df_temp.Recall):.4f} ({stats.sem(df_temp.Recall):,.4f})')\n",
    "\n",
    "    df_temp = df_wegs_target_variants[df_wegs_target_variants.HGID == hgid]\n",
    "    df_temp = df_temp[(df_temp.Type == 'INDEL') & (df_temp.Filter == 'ALL') & (df_temp.Plexing == 4) & (df_temp.WGS_depth == 2)]\n",
    "    data['WEGS 4P2X TP'].append(f'{np.mean(df_temp.TP):,.0f} ({stats.sem(df_temp.TP):,.0f})')\n",
    "    data['WEGS 4P2X FP'].append(f'{np.mean(df_temp.FP):,.0f} ({stats.sem(df_temp.FP):,.0f})')\n",
    "    data['WEGS 4P2X FN'].append(f'{np.mean(df_temp.FN):,.0f} ({stats.sem(df_temp.FN):,.0f})')\n",
    "    data['WEGS 4P2X Precision'].append(f'{np.mean(df_temp.Precision):.4f} ({stats.sem(df_temp.Precision):,.4f})')\n",
    "    data['WEGS 4P2X Recall'].append(f'{np.mean(df_temp.Recall):.4f} ({stats.sem(df_temp.Recall):,.4f})')\n",
    "    \n",
    "    df_temp = df_wegs_target_variants[df_wegs_target_variants.HGID == hgid]\n",
    "    df_temp = df_temp[(df_temp.Type == 'INDEL') & (df_temp.Filter == 'ALL') & (df_temp.Plexing == 8) & (df_temp.WGS_depth == 5)]\n",
    "    data['WEGS 8P5X TP'].append(f'{np.mean(df_temp.TP):,.0f} ({stats.sem(df_temp.TP):,.0f})')\n",
    "    data['WEGS 8P5X FP'].append(f'{np.mean(df_temp.FP):,.0f} ({stats.sem(df_temp.FP):,.0f})')\n",
    "    data['WEGS 8P5X FN'].append(f'{np.mean(df_temp.FN):,.0f} ({stats.sem(df_temp.FN):,.0f})')\n",
    "    data['WEGS 8P5X Precision'].append(f'{np.mean(df_temp.Precision):.4f} ({stats.sem(df_temp.Precision):,.4f})')\n",
    "    data['WEGS 8P5X Recall'].append(f'{np.mean(df_temp.Recall):.4f} ({stats.sem(df_temp.Recall):,.4f})')\n",
    "    \n",
    "pd.DataFrame(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cdb5c45-e57c-4c9a-926a-8bba66efe948",
   "metadata": {},
   "source": [
    "### TABLE: Average genome-wide variant recall and precision rates in 30X WGS and WEGS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "77e5df3b-e6ec-4684-8034-ee70492bc814",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Label</th>\n",
       "      <th>SNV TP</th>\n",
       "      <th>SNV FP</th>\n",
       "      <th>SNV FN</th>\n",
       "      <th>SNV Precision</th>\n",
       "      <th>SNV Recall</th>\n",
       "      <th>InDel TP</th>\n",
       "      <th>InDel FP</th>\n",
       "      <th>InDel FN</th>\n",
       "      <th>InDel Precision</th>\n",
       "      <th>InDel Recall</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>WGS 30X</td>\n",
       "      <td>3,309,667 (4,295)</td>\n",
       "      <td>15,268 (174)</td>\n",
       "      <td>26,097 (82)</td>\n",
       "      <td>0.9954 (0.0001)</td>\n",
       "      <td>0.9922 (0.0000)</td>\n",
       "      <td>500,728 (2,382)</td>\n",
       "      <td>5,605 (35)</td>\n",
       "      <td>9,757 (111)</td>\n",
       "      <td>0.9889 (0.0001)</td>\n",
       "      <td>0.9809 (0.0002)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>WEGS 4P2X</td>\n",
       "      <td>1,333,840 (44,603)</td>\n",
       "      <td>303,452 (4,673)</td>\n",
       "      <td>2,001,925 (48,123)</td>\n",
       "      <td>0.8137 (0.0033)</td>\n",
       "      <td>0.4000 (0.0138)</td>\n",
       "      <td>132,918 (5,359)</td>\n",
       "      <td>56,441 (1,663)</td>\n",
       "      <td>377,567 (6,801)</td>\n",
       "      <td>0.7010 (0.0025)</td>\n",
       "      <td>0.2607 (0.0112)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>WEGS 4P5X</td>\n",
       "      <td>1,909,331 (58,921)</td>\n",
       "      <td>240,392 (3,054)</td>\n",
       "      <td>1,426,434 (62,566)</td>\n",
       "      <td>0.8869 (0.0043)</td>\n",
       "      <td>0.5726 (0.0183)</td>\n",
       "      <td>201,593 (7,908)</td>\n",
       "      <td>68,206 (1,221)</td>\n",
       "      <td>308,892 (9,249)</td>\n",
       "      <td>0.7457 (0.0041)</td>\n",
       "      <td>0.3954 (0.0166)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>WEGS 8P2X</td>\n",
       "      <td>1,326,086 (39,124)</td>\n",
       "      <td>307,039 (4,029)</td>\n",
       "      <td>2,008,327 (42,308)</td>\n",
       "      <td>0.8109 (0.0029)</td>\n",
       "      <td>0.3979 (0.0121)</td>\n",
       "      <td>131,783 (4,687)</td>\n",
       "      <td>56,193 (1,460)</td>\n",
       "      <td>378,097 (6,013)</td>\n",
       "      <td>0.7001 (0.0022)</td>\n",
       "      <td>0.2588 (0.0098)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>WEGS 8P5X</td>\n",
       "      <td>1,914,743 (51,710)</td>\n",
       "      <td>240,678 (2,681)</td>\n",
       "      <td>1,419,670 (54,999)</td>\n",
       "      <td>0.8870 (0.0038)</td>\n",
       "      <td>0.5745 (0.0161)</td>\n",
       "      <td>201,995 (6,927)</td>\n",
       "      <td>68,102 (1,067)</td>\n",
       "      <td>307,885 (8,174)</td>\n",
       "      <td>0.7463 (0.0036)</td>\n",
       "      <td>0.3967 (0.0146)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Label              SNV TP           SNV FP              SNV FN  \\\n",
       "0    WGS 30X   3,309,667 (4,295)     15,268 (174)         26,097 (82)   \n",
       "1  WEGS 4P2X  1,333,840 (44,603)  303,452 (4,673)  2,001,925 (48,123)   \n",
       "2  WEGS 4P5X  1,909,331 (58,921)  240,392 (3,054)  1,426,434 (62,566)   \n",
       "3  WEGS 8P2X  1,326,086 (39,124)  307,039 (4,029)  2,008,327 (42,308)   \n",
       "4  WEGS 8P5X  1,914,743 (51,710)  240,678 (2,681)  1,419,670 (54,999)   \n",
       "\n",
       "     SNV Precision       SNV Recall         InDel TP        InDel FP  \\\n",
       "0  0.9954 (0.0001)  0.9922 (0.0000)  500,728 (2,382)      5,605 (35)   \n",
       "1  0.8137 (0.0033)  0.4000 (0.0138)  132,918 (5,359)  56,441 (1,663)   \n",
       "2  0.8869 (0.0043)  0.5726 (0.0183)  201,593 (7,908)  68,206 (1,221)   \n",
       "3  0.8109 (0.0029)  0.3979 (0.0121)  131,783 (4,687)  56,193 (1,460)   \n",
       "4  0.8870 (0.0038)  0.5745 (0.0161)  201,995 (6,927)  68,102 (1,067)   \n",
       "\n",
       "          InDel FN  InDel Precision     InDel Recall  \n",
       "0      9,757 (111)  0.9889 (0.0001)  0.9809 (0.0002)  \n",
       "1  377,567 (6,801)  0.7010 (0.0025)  0.2607 (0.0112)  \n",
       "2  308,892 (9,249)  0.7457 (0.0041)  0.3954 (0.0166)  \n",
       "3  378,097 (6,013)  0.7001 (0.0022)  0.2588 (0.0098)  \n",
       "4  307,885 (8,174)  0.7463 (0.0036)  0.3967 (0.0146)  "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dfs = [\n",
    "    ( 'WGS 30X', df_wgs30x_variants), \n",
    "    ( 'WEGS 4P2X', df_wegs_variants[(df_wegs_variants.Plexing == 4) & (df_wegs_variants.WGS_depth == 2)]),\n",
    "    ( 'WEGS 4P5X', df_wegs_variants[(df_wegs_variants.Plexing == 4) & (df_wegs_variants.WGS_depth == 5)]),\n",
    "    ( 'WEGS 8P2X', df_wegs_variants[(df_wegs_variants.Plexing == 8) & (df_wegs_variants.WGS_depth == 2)]),\n",
    "    ( 'WEGS 8P5X', df_wegs_variants[(df_wegs_variants.Plexing == 8) & (df_wegs_variants.WGS_depth == 5)])\n",
    "]\n",
    "\n",
    "data = {\n",
    "    'Label': [],\n",
    "    'SNV TP': [],\n",
    "    'SNV FP': [],\n",
    "    'SNV FN': [],\n",
    "    'SNV Precision': [],\n",
    "    'SNV Recall': [],\n",
    "    'InDel TP': [],\n",
    "    'InDel FP': [],\n",
    "    'InDel FN': [],\n",
    "    'InDel Precision': [],\n",
    "    'InDel Recall': []\n",
    "}\n",
    "\n",
    "for label, df in dfs:\n",
    "    data['Label'].append(label)\n",
    "    \n",
    "    df_temp = df[(df.Type == 'SNP') & (df.Filter == 'ALL')]\n",
    "    data['SNV TP'].append(f'{np.mean(df_temp.TP):,.0f} ({stats.sem(df_temp.TP):,.0f})')\n",
    "    data['SNV FP'].append(f'{np.mean(df_temp.FP):,.0f} ({stats.sem(df_temp.FP):,.0f})')\n",
    "    data['SNV FN'].append(f'{np.mean(df_temp.FN):,.0f} ({stats.sem(df_temp.FN):,.0f})')\n",
    "    data['SNV Precision'].append(f'{np.mean(df_temp.Precision):.4f} ({stats.sem(df_temp.Precision):,.4f})')\n",
    "    data['SNV Recall'].append(f'{np.mean(df_temp.Recall):.4f} ({stats.sem(df_temp.Recall):,.4f})')\n",
    "\n",
    "    df_temp = df[(df.Type == 'INDEL') & (df.Filter == 'ALL')]\n",
    "    data['InDel TP'].append(f'{np.mean(df_temp.TP):,.0f} ({stats.sem(df_temp.TP):,.0f})')\n",
    "    data['InDel FP'].append(f'{np.mean(df_temp.FP):,.0f} ({stats.sem(df_temp.FP):,.0f})')\n",
    "    data['InDel FN'].append(f'{np.mean(df_temp.FN):,.0f} ({stats.sem(df_temp.FN):,.0f})')\n",
    "    data['InDel Precision'].append(f'{np.mean(df_temp.Precision):.4f} ({stats.sem(df_temp.Precision):,.4f})')\n",
    "    data['InDel Recall'].append(f'{np.mean(df_temp.Recall):.4f} ({stats.sem(df_temp.Recall):,.4f})')\n",
    "\n",
    "pd.DataFrame(data)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5fe2c21-18f0-4929-aae5-31e2e6fce8f7",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9",
   "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}