1KG_HGDP_array_imputation_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4a32ff5c-98f2-478e-bbfc-42b151fc4e8e",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17fa91c0-34ee-4301-bdaa-a046f1293ecf",
   "metadata": {},
   "source": [
    "## 1. Help functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7ed167cd-f9e5-4b64-aaca-322d9ef62e1e",
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec4bed38-36a4-437d-9608-65c30e55a358",
   "metadata": {},
   "source": [
    "## 2. Input files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "051f9bff-7cba-49bb-9d57-00edc75d373d",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "SAMPLE_ANNOTATIONS_FILE = 'sample_annotations.txt'\n",
    "ARRAY_BENCHMARK_DIR = '1KG_HGDP_array_imputation_analysis'\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c1b2fbe-d705-4b4e-9903-d74f2b2f5f79",
   "metadata": {},
   "source": [
    "## 3. Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e782394a-a3ce-42b5-87fc-a36a5ec3dba5",
   "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",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a67f8cf-7cfc-4ba9-a7b6-0c3b01f7a517",
   "metadata": {},
   "source": [
    "### 3.1 Load 1KG+HGDP array imputed variants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3f688ef8-cf3b-4686-9df7-e4487d51a022",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_array_imputed_variants = []\n",
    "for row_index, row in df_samples[['HGID', 'Family']].drop_duplicates().iterrows():\n",
    "    filepath = f'{ARRAY_BENCHMARK_DIR}/Genomewide/{row.HGID}_GRCh37_1_22_v4.2.1_GSAv3_positions.impuded.dose.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_array_imputed_variants.append(df)\n",
    "df_array_imputed_variants = pd.concat(df_array_imputed_variants).reset_index(drop = True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdf18fd4-3233-40ed-b1c5-7ca1a2d8e1eb",
   "metadata": {},
   "source": [
    "## 4. Figures and tables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5515344-6c62-4064-b7cc-4188a03c2b29",
   "metadata": {},
   "source": [
    "### TABLE: Precision and recall rates of variants imputed using local reference panel genome-wide."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f7e2b67d-18fa-44a4-954c-7f0e65a99ffa",
   "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>Imputed SNVs TP</th>\n",
       "      <th>Imputed SNVs FN</th>\n",
       "      <th>Imputed SNVs FP</th>\n",
       "      <th>Imputed SNVs Precision</th>\n",
       "      <th>Imputed SNVs Recall</th>\n",
       "      <th>Imputed InDels TP</th>\n",
       "      <th>Imputed InDels FN</th>\n",
       "      <th>Imputed InDels FP</th>\n",
       "      <th>Imputed InDels Precision</th>\n",
       "      <th>Imputed InDels Recall</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>HG002</td>\n",
       "      <td>2,942,195</td>\n",
       "      <td>410,484</td>\n",
       "      <td>316,718</td>\n",
       "      <td>0.9028</td>\n",
       "      <td>0.8776</td>\n",
       "      <td>326,360</td>\n",
       "      <td>196,029</td>\n",
       "      <td>47,144</td>\n",
       "      <td>0.8738</td>\n",
       "      <td>0.6247</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>HG003</td>\n",
       "      <td>2,908,634</td>\n",
       "      <td>405,509</td>\n",
       "      <td>317,430</td>\n",
       "      <td>0.9016</td>\n",
       "      <td>0.8776</td>\n",
       "      <td>317,001</td>\n",
       "      <td>183,792</td>\n",
       "      <td>46,494</td>\n",
       "      <td>0.8721</td>\n",
       "      <td>0.6330</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>HG004</td>\n",
       "      <td>2,934,901</td>\n",
       "      <td>405,570</td>\n",
       "      <td>313,069</td>\n",
       "      <td>0.9036</td>\n",
       "      <td>0.8786</td>\n",
       "      <td>320,526</td>\n",
       "      <td>187,748</td>\n",
       "      <td>45,098</td>\n",
       "      <td>0.8767</td>\n",
       "      <td>0.6306</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Sample Imputed SNVs TP Imputed SNVs FN Imputed SNVs FP  \\\n",
       "0  HG002       2,942,195         410,484         316,718   \n",
       "1  HG003       2,908,634         405,509         317,430   \n",
       "2  HG004       2,934,901         405,570         313,069   \n",
       "\n",
       "  Imputed SNVs Precision Imputed SNVs Recall Imputed InDels TP  \\\n",
       "0                 0.9028              0.8776           326,360   \n",
       "1                 0.9016              0.8776           317,001   \n",
       "2                 0.9036              0.8786           320,526   \n",
       "\n",
       "  Imputed InDels FN Imputed InDels FP Imputed InDels Precision  \\\n",
       "0           196,029            47,144                   0.8738   \n",
       "1           183,792            46,494                   0.8721   \n",
       "2           187,748            45,098                   0.8767   \n",
       "\n",
       "  Imputed InDels Recall  \n",
       "0                0.6247  \n",
       "1                0.6330  \n",
       "2                0.6306  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = {\n",
    "    'Sample': [],\n",
    "    'Imputed SNVs TP': [],\n",
    "    'Imputed SNVs FN': [],\n",
    "    'Imputed SNVs FP': [],\n",
    "    'Imputed SNVs Precision': [],\n",
    "    'Imputed SNVs Recall': [],\n",
    "    'Imputed InDels TP': [],\n",
    "    'Imputed InDels FN': [],\n",
    "    'Imputed InDels FP': [],\n",
    "    'Imputed InDels Precision': [],\n",
    "    'Imputed InDels Recall': [],\n",
    "}\n",
    "\n",
    "for hgid, naid in hgid2naid.items():\n",
    "    data['Sample'].append(hgid)\n",
    "    \n",
    "    df_temp = df_array_imputed_variants[df_array_imputed_variants.HGID == hgid]\n",
    "    df_temp = df_temp[(df_temp.Type == 'SNP') & (df_temp.Filter == 'ALL')][['TP', 'FN', 'FP', 'Recall', 'Precision']]\n",
    "    data['Imputed SNVs TP'].append(f'{df_temp.TP.values[0]:,}')\n",
    "    data['Imputed SNVs FP'].append(f'{df_temp.FP.values[0]:,}')\n",
    "    data['Imputed SNVs FN'].append(f'{df_temp.FN.values[0]:,}')\n",
    "    data['Imputed SNVs Precision'].append(f'{df_temp.Precision.values[0]:.4f}')\n",
    "    data['Imputed SNVs Recall'].append(f'{df_temp.Recall.values[0]:.4f}')\n",
    "    \n",
    "    df_temp = df_array_imputed_variants[df_array_imputed_variants.HGID == hgid]\n",
    "    df_temp = df_temp[(df_temp.Type == 'INDEL') & (df_temp.Filter == 'ALL')][['TP', 'FN', 'FP', 'Recall', 'Precision']]\n",
    "    data['Imputed InDels TP'].append(f'{df_temp.TP.values[0]:,}')\n",
    "    data['Imputed InDels FP'].append(f'{df_temp.FP.values[0]:,}')\n",
    "    data['Imputed InDels FN'].append(f'{df_temp.FN.values[0]:,}')\n",
    "    data['Imputed InDels Precision'].append(f'{df_temp.Precision.values[0]:.4f}')\n",
    "    data['Imputed InDels Recall'].append(f'{df_temp.Recall.values[0]:.4f}')\n",
    "    \n",
    "pd.DataFrame(data)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b810d763-6940-4031-a1e7-6213297d5e1a",
   "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
}