Prostate Cancer Risk Screening Using Germline Variant Data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prostate Cancer Risk Screening Using Germline Variant Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Collecting scikit-allel\n",
      "  Downloading scikit_allel-1.3.5-cp38-cp38-macosx_10_9_x86_64.whl (1.4 MB)\n",
      "\u001b[K     |████████████████████████████████| 1.4 MB 6.4 MB/s eta 0:00:01\n",
      "\u001b[?25hRequirement already satisfied: numpy in /Users/crystalshin/opt/anaconda3/lib/python3.8/site-packages (from scikit-allel) (1.19.2)\n",
      "Requirement already satisfied: dask[array] in /Users/crystalshin/opt/anaconda3/lib/python3.8/site-packages (from scikit-allel) (2.30.0)\n",
      "Requirement already satisfied: pyyaml in /Users/crystalshin/opt/anaconda3/lib/python3.8/site-packages (from dask[array]->scikit-allel) (5.3.1)\n",
      "Requirement already satisfied: toolz>=0.8.2 in /Users/crystalshin/opt/anaconda3/lib/python3.8/site-packages (from dask[array]->scikit-allel) (0.11.1)\n",
      "Installing collected packages: scikit-allel\n",
      "Successfully installed scikit-allel-1.3.5\n",
      "\u001b[33mWARNING: You are using pip version 21.0.1; however, version 22.3.1 is available.\n",
      "You should consider upgrading via the '/Users/crystalshin/opt/anaconda3/bin/python -m pip install --upgrade pip' command.\u001b[0m\n",
      "Collecting zarr\n",
      "  Downloading zarr-2.13.3-py3-none-any.whl (191 kB)\n",
      "\u001b[K     |████████████████████████████████| 191 kB 10.3 MB/s eta 0:00:01\n",
      "\u001b[?25hCollecting fasteners\n",
      "  Downloading fasteners-0.18-py3-none-any.whl (18 kB)\n",
      "Collecting numcodecs>=0.10.0\n",
      "  Downloading numcodecs-0.11.0-cp38-cp38-macosx_10_9_x86_64.whl (1.4 MB)\n",
      "\u001b[K     |████████████████████████████████| 1.4 MB 9.7 MB/s eta 0:00:01\n",
      "\u001b[?25hRequirement already satisfied: numpy>=1.7 in /Users/crystalshin/opt/anaconda3/lib/python3.8/site-packages (from zarr) (1.19.2)\n",
      "Collecting asciitree\n",
      "  Downloading asciitree-0.3.3.tar.gz (4.0 kB)\n",
      "Requirement already satisfied: entrypoints in /Users/crystalshin/opt/anaconda3/lib/python3.8/site-packages (from numcodecs>=0.10.0->zarr) (0.3)\n",
      "Building wheels for collected packages: asciitree\n",
      "  Building wheel for asciitree (setup.py) ... \u001b[?25ldone\n",
      "\u001b[?25h  Created wheel for asciitree: filename=asciitree-0.3.3-py3-none-any.whl size=5035 sha256=cad9ecd67baa7eaf1946cd4a7b127b564a18cbcc6a8817d04858f7eebdd8876d\n",
      "  Stored in directory: /Users/crystalshin/Library/Caches/pip/wheels/a3/d7/75/19cd0d2a893cad4bb0b2b16dd572ad2916d19c0d5ee9612511\n",
      "Successfully built asciitree\n",
      "Installing collected packages: numcodecs, fasteners, asciitree, zarr\n",
      "Successfully installed asciitree-0.3.3 fasteners-0.18 numcodecs-0.11.0 zarr-2.13.3\n",
      "\u001b[33mWARNING: You are using pip version 21.0.1; however, version 22.3.1 is available.\n",
      "You should consider upgrading via the '/Users/crystalshin/opt/anaconda3/bin/python -m pip install --upgrade pip' command.\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "!pip install scikit-allel\n",
    "!pip install zarr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "cwd = os.getcwd()\n",
    "DATA = cwd + '/data'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Identifying Pathogenic Variants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can identify potential disease-causing variants using germline variant data. 'Clinvar' is a clinical variants database for pathogenic disease variant annotation. In this analysis, I will identify pathogenic germline variants that are likely to result in prostate cancer in 1000 genome project subjects. 1000 genome project contains variants from whole genome sequencing of 2,504 volunteers of 26 sub-population. I have downloaded the Clinvar database as a VCF file. VCF (variant call format) is the summary of all variants compared to the reference genome, and each row contains the annotation for each variant. We are using variant mapped to human reference genome version GRCh37 since the 1000 genome database is also mapped to the same version. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[read_vcf] 65536 rows in 0.66s; chunk in 0.66s (98744 rows/s); 1\u0000:150469297\n",
      "[read_vcf] 131072 rows in 1.31s; chunk in 0.65s (100707 rows/s); 2\u0000:26467440\n",
      "[read_vcf] 196608 rows in 1.97s; chunk in 0.66s (99585 rows/s); 2\u0000:165946783\n",
      "[read_vcf] 262144 rows in 2.66s; chunk in 0.69s (95524 rows/s); 2\u0000:233113958\n",
      "[read_vcf] 327680 rows in 3.35s; chunk in 0.69s (94696 rows/s); 3\u0000:128627919\n",
      "[read_vcf] 393216 rows in 4.03s; chunk in 0.68s (96557 rows/s); 4\u0000:184615101\n",
      "[read_vcf] 458752 rows in 4.66s; chunk in 0.63s (103692 rows/s); 5\u0000:148377020\n",
      "[read_vcf] 524288 rows in 5.33s; chunk in 0.67s (97701 rows/s); 6\u0000:129774241\n",
      "[read_vcf] 589824 rows in 5.97s; chunk in 0.64s (102314 rows/s); 7\u0000:124532342\n",
      "[read_vcf] 655360 rows in 6.66s; chunk in 0.68s (95891 rows/s); 8\u0000:145739033\n",
      "[read_vcf] 720896 rows in 7.28s; chunk in 0.62s (105669 rows/s); 9\u0000:139418044\n",
      "[read_vcf] 786432 rows in 7.88s; chunk in 0.60s (108906 rows/s); 11\u0000:3877646\n",
      "[read_vcf] 851968 rows in 8.52s; chunk in 0.65s (101525 rows/s); 11\u0000:108218000\n",
      "[read_vcf] 917504 rows in 9.16s; chunk in 0.64s (102565 rows/s); 12\u0000:112892622\n",
      "[read_vcf] 983040 rows in 9.81s; chunk in 0.65s (101412 rows/s); 14\u0000:64457839\n",
      "[read_vcf] 1048576 rows in 10.45s; chunk in 0.64s (101987 rows/s); 15\u0000:80465343\n",
      "[read_vcf] 1114112 rows in 11.12s; chunk in 0.66s (98663 rows/s); 16\u0000:56899006\n",
      "[read_vcf] 1179648 rows in 11.83s; chunk in 0.71s (91703 rows/s); 17\u0000:33446172\n",
      "[read_vcf] 1245184 rows in 12.52s; chunk in 0.69s (95209 rows/s); 18\u0000:19070202\n",
      "[read_vcf] 1310720 rows in 13.17s; chunk in 0.65s (100337 rows/s); 19\u0000:39061277\n",
      "[read_vcf] 1376256 rows in 13.80s; chunk in 0.62s (105143 rows/s); 21\u0000:46929440\n",
      "[read_vcf] 1441792 rows in 14.48s; chunk in 0.68s (96189 rows/s); X\u0000:70329122\n",
      "[read_vcf] 1475077 rows in 14.83s; chunk in 0.35s (94499 rows/s)\n",
      "[read_vcf] all done (99469 rows/s)\n"
     ]
    }
   ],
   "source": [
    "import allel\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "clinvar_file = DATA + \"/clinvar_20220910_grch37.vcf.gz\"\n",
    "clinvar = allel.read_vcf(str(clinvar_file), fields=\"*\", log=sys.stdout)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we would like to find all variants that are annotated as pathogenic. For more thorough search in real life, we would want to search all chromosomes for pathogenic variants, but for simplicity of this analysis, we will restrict the boundary to chromosome X. \n",
    "\n",
    "Pathogenicity can be found in column variants/CLNSIG, and the disease name could be found in variants/CLNDN. I used regular expression to find all entries containing the word \"Pathogenic\" in the column `variants/CLNDN`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "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>CHROM</th>\n",
       "      <th>POS</th>\n",
       "      <th>REF</th>\n",
       "      <th>GENEINFO</th>\n",
       "      <th>CLNDN</th>\n",
       "      <th>CLNSIG</th>\n",
       "      <th>AF_TGP</th>\n",
       "      <th>ALT</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>X</td>\n",
       "      <td>66765148</td>\n",
       "      <td>T</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate|Kennedy_disease|An...</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>TTGC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>X</td>\n",
       "      <td>66765196</td>\n",
       "      <td>C</td>\n",
       "      <td>AR:367|LOC109504725:109504725</td>\n",
       "      <td>Malignant_tumor_of_prostate</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>X</td>\n",
       "      <td>66765463</td>\n",
       "      <td>G</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate|Androgen_resistanc...</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>X</td>\n",
       "      <td>66765802</td>\n",
       "      <td>C</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate|not_provided</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>X</td>\n",
       "      <td>66766013</td>\n",
       "      <td>C</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>X</td>\n",
       "      <td>66766051</td>\n",
       "      <td>G</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>X</td>\n",
       "      <td>66766163</td>\n",
       "      <td>C</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>X</td>\n",
       "      <td>66766183</td>\n",
       "      <td>T</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>X</td>\n",
       "      <td>66766196</td>\n",
       "      <td>C</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>X</td>\n",
       "      <td>66863125</td>\n",
       "      <td>G</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>X</td>\n",
       "      <td>66863132</td>\n",
       "      <td>G</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>X</td>\n",
       "      <td>66905875</td>\n",
       "      <td>A</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Non-obstructive_azoospermia|Malignant_tumor_of...</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>X</td>\n",
       "      <td>66931301</td>\n",
       "      <td>G</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Prostate_cancer</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>X</td>\n",
       "      <td>66937326</td>\n",
       "      <td>G</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate|Prostate_cancer_su...</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>0.00053</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>X</td>\n",
       "      <td>66937330</td>\n",
       "      <td>C</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>X</td>\n",
       "      <td>66937337</td>\n",
       "      <td>G</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate|Prostate_cancer</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>0.00053</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>X</td>\n",
       "      <td>66942818</td>\n",
       "      <td>G</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate|Hypospadias_1</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>X</td>\n",
       "      <td>66943543</td>\n",
       "      <td>C</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Prostate_cancer</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>T</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>X</td>\n",
       "      <td>66943552</td>\n",
       "      <td>A</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Prostate_cancer</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>X</td>\n",
       "      <td>66943553</td>\n",
       "      <td>C</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Prostate_cancer</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>X</td>\n",
       "      <td>66943628</td>\n",
       "      <td>A</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Prostate_cancer</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "      <td>G</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   CHROM       POS REF                       GENEINFO  \\\n",
       "0      X  66765148   T                         AR:367   \n",
       "1      X  66765196   C  AR:367|LOC109504725:109504725   \n",
       "2      X  66765463   G                         AR:367   \n",
       "3      X  66765802   C                         AR:367   \n",
       "4      X  66766013   C                         AR:367   \n",
       "5      X  66766051   G                         AR:367   \n",
       "6      X  66766163   C                         AR:367   \n",
       "7      X  66766183   T                         AR:367   \n",
       "8      X  66766196   C                         AR:367   \n",
       "9      X  66863125   G                         AR:367   \n",
       "10     X  66863132   G                         AR:367   \n",
       "11     X  66905875   A                         AR:367   \n",
       "12     X  66931301   G                         AR:367   \n",
       "13     X  66937326   G                         AR:367   \n",
       "14     X  66937330   C                         AR:367   \n",
       "15     X  66937337   G                         AR:367   \n",
       "16     X  66942818   G                         AR:367   \n",
       "17     X  66943543   C                         AR:367   \n",
       "18     X  66943552   A                         AR:367   \n",
       "19     X  66943553   C                         AR:367   \n",
       "20     X  66943628   A                         AR:367   \n",
       "\n",
       "                                                CLNDN      CLNSIG   AF_TGP  \\\n",
       "0   Malignant_tumor_of_prostate|Kennedy_disease|An...  Pathogenic      NaN   \n",
       "1                         Malignant_tumor_of_prostate  Pathogenic      NaN   \n",
       "2   Malignant_tumor_of_prostate|Androgen_resistanc...  Pathogenic      NaN   \n",
       "3            Malignant_tumor_of_prostate|not_provided  Pathogenic      NaN   \n",
       "4                         Malignant_tumor_of_prostate  Pathogenic      NaN   \n",
       "5                         Malignant_tumor_of_prostate  Pathogenic      NaN   \n",
       "6                         Malignant_tumor_of_prostate  Pathogenic      NaN   \n",
       "7                         Malignant_tumor_of_prostate  Pathogenic      NaN   \n",
       "8                         Malignant_tumor_of_prostate  Pathogenic      NaN   \n",
       "9                         Malignant_tumor_of_prostate  Pathogenic      NaN   \n",
       "10                        Malignant_tumor_of_prostate  Pathogenic      NaN   \n",
       "11  Non-obstructive_azoospermia|Malignant_tumor_of...  Pathogenic      NaN   \n",
       "12                                    Prostate_cancer  Pathogenic      NaN   \n",
       "13  Malignant_tumor_of_prostate|Prostate_cancer_su...  Pathogenic  0.00053   \n",
       "14                        Malignant_tumor_of_prostate  Pathogenic      NaN   \n",
       "15        Malignant_tumor_of_prostate|Prostate_cancer  Pathogenic  0.00053   \n",
       "16          Malignant_tumor_of_prostate|Hypospadias_1  Pathogenic      NaN   \n",
       "17                                    Prostate_cancer  Pathogenic      NaN   \n",
       "18                                    Prostate_cancer  Pathogenic      NaN   \n",
       "19                                    Prostate_cancer  Pathogenic      NaN   \n",
       "20                                    Prostate_cancer  Pathogenic      NaN   \n",
       "\n",
       "     ALT  \n",
       "0   TTGC  \n",
       "1      T  \n",
       "2      A  \n",
       "3      T  \n",
       "4      T  \n",
       "5      C  \n",
       "6      G  \n",
       "7      C  \n",
       "8      T  \n",
       "9      T  \n",
       "10     C  \n",
       "11     G  \n",
       "12     A  \n",
       "13     T  \n",
       "14     G  \n",
       "15     A  \n",
       "16     A  \n",
       "17     T  \n",
       "18     G  \n",
       "19     G  \n",
       "20     G  "
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import re\n",
    "import pandas as pd\n",
    "\n",
    "chrom_x_idx = clinvar[\"variants/CHROM\"] == 'X'\n",
    "pathogenic_idx = np.array([x == \"Pathogenic\" for x in clinvar[\"variants/CLNSIG\"]])\n",
    "prostate_idx = np.array([\"prostate\" in x.lower() for x in clinvar[\"variants/CLNDN\"]])\n",
    "\n",
    "idx_set = chrom_x_idx & pathogenic_idx & prca_idx\n",
    "\n",
    "pathogenic_df = pd.DataFrame(\n",
    "    dict([(k, clinvar[\"variants/\" + k][idx_set]) \n",
    "    for k in [\"CHROM\", \"POS\", \"REF\", \"GENEINFO\", \"CLNDN\", \"CLNSIG\", \"AF_TGP\"]])\n",
    ")\n",
    "pathogenic_df[\"ALT\"] = clinvar[\"variants/ALT\"][idx_set, 0]\n",
    "\n",
    "pathogenic_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There were 21 pathogenic variants to prostate cancer in chromosome X. These variants came from the AR:367 gene.\n",
    "\n",
    "The AR gene makes the androgen receptor. Androgen receptor binds to androgenic hormones including testosterone. Its functions include being a DNA binding transcription factor that regulates gene expression and also development and maintenance of the male sexual phenotype. According to literature, the androgen receptor is linked with pathogenesis and progression of prostate cancer. An overexpression of AR leads to progression of prostate cancer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Identifying subjects carrying pathogenic variants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, I will be screening the 1000 genome project to identify subjects carrying prostate cancer pathogenic variants. First, I will shortlist the variants that both appear in the variants identified in the Clinvar database and 1000 genome project. I will be converting vcf file to a zarr file to make handling of large data more manageable.\n",
    "\n",
    "I will also find the maximum number of alternate alleles across all positions, since this information could be useful when I try to calculate alternate allele frequencies of different alleles. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[read_vcf] 65536 rows in 3.60s; chunk in 3.60s (18210 rows/s); X\u0000:1734128\n",
      "[read_vcf] 131072 rows in 6.91s; chunk in 3.31s (19777 rows/s); X\u0000:3915040\n",
      "[read_vcf] 196608 rows in 9.86s; chunk in 2.95s (22219 rows/s); X\u0000:6474138\n",
      "[read_vcf] 262144 rows in 12.98s; chunk in 3.12s (21032 rows/s); X\u0000:9519142\n",
      "[read_vcf] 327680 rows in 16.04s; chunk in 3.06s (21388 rows/s); X\u0000:12551848\n",
      "[read_vcf] 393216 rows in 19.11s; chunk in 3.07s (21359 rows/s); X\u0000:15550284\n",
      "[read_vcf] 458752 rows in 22.03s; chunk in 2.92s (22436 rows/s); X\u0000:18776831\n",
      "[read_vcf] 524288 rows in 24.96s; chunk in 2.93s (22385 rows/s); X\u0000:21881750\n",
      "[read_vcf] 589824 rows in 28.30s; chunk in 3.34s (19637 rows/s); X\u0000:24574486\n",
      "[read_vcf] 655360 rows in 31.27s; chunk in 2.98s (22010 rows/s); X\u0000:27659891\n",
      "[read_vcf] 720896 rows in 34.11s; chunk in 2.84s (23101 rows/s); X\u0000:30671296\n",
      "[read_vcf] 786432 rows in 37.00s; chunk in 2.89s (22652 rows/s); X\u0000:33194562\n",
      "[read_vcf] 851968 rows in 40.10s; chunk in 3.10s (21131 rows/s); X\u0000:36031556\n",
      "[read_vcf] 917504 rows in 43.36s; chunk in 3.25s (20140 rows/s); X\u0000:39125534\n",
      "[read_vcf] 983040 rows in 46.72s; chunk in 3.36s (19497 rows/s); X\u0000:42186473\n",
      "[read_vcf] 1048576 rows in 49.65s; chunk in 2.93s (22392 rows/s); X\u0000:45309329\n",
      "[read_vcf] 1114112 rows in 52.77s; chunk in 3.13s (20964 rows/s); X\u0000:48359153\n",
      "[read_vcf] 1179648 rows in 56.20s; chunk in 3.42s (19137 rows/s); X\u0000:51881411\n",
      "[read_vcf] 1245184 rows in 59.09s; chunk in 2.89s (22673 rows/s); X\u0000:55705421\n",
      "[read_vcf] 1310720 rows in 61.96s; chunk in 2.87s (22809 rows/s); X\u0000:57996652\n",
      "[read_vcf] 1376256 rows in 64.88s; chunk in 2.91s (22485 rows/s); X\u0000:63082608\n",
      "[read_vcf] 1441792 rows in 67.70s; chunk in 2.83s (23187 rows/s); X\u0000:65358749\n",
      "[read_vcf] 1507328 rows in 70.73s; chunk in 3.02s (21668 rows/s); X\u0000:67634861\n",
      "[read_vcf] 1572864 rows in 73.58s; chunk in 2.86s (22937 rows/s); X\u0000:70712278\n",
      "[read_vcf] 1638400 rows in 76.43s; chunk in 2.84s (23038 rows/s); X\u0000:74091783\n",
      "[read_vcf] 1703936 rows in 79.28s; chunk in 2.85s (22982 rows/s); X\u0000:76280043\n",
      "[read_vcf] 1769472 rows in 82.14s; chunk in 2.86s (22904 rows/s); X\u0000:78926545\n",
      "[read_vcf] 1835008 rows in 85.45s; chunk in 3.31s (19788 rows/s); X\u0000:81422228\n",
      "[read_vcf] 1900544 rows in 88.39s; chunk in 2.94s (22277 rows/s); X\u0000:83965308\n",
      "[read_vcf] 1966080 rows in 91.27s; chunk in 2.87s (22829 rows/s); X\u0000:86441870\n",
      "[read_vcf] 2031616 rows in 94.21s; chunk in 2.94s (22290 rows/s); X\u0000:89012763\n",
      "[read_vcf] 2097152 rows in 97.23s; chunk in 3.02s (21676 rows/s); X\u0000:91967747\n",
      "[read_vcf] 2162688 rows in 100.11s; chunk in 2.88s (22771 rows/s); X\u0000:94617360\n",
      "[read_vcf] 2228224 rows in 103.02s; chunk in 2.91s (22520 rows/s); X\u0000:97514551\n",
      "[read_vcf] 2293760 rows in 105.91s; chunk in 2.89s (22655 rows/s); X\u0000:100546234\n",
      "[read_vcf] 2359296 rows in 108.84s; chunk in 2.93s (22377 rows/s); X\u0000:104064469\n",
      "[read_vcf] 2424832 rows in 111.68s; chunk in 2.84s (23036 rows/s); X\u0000:107375123\n",
      "[read_vcf] 2490368 rows in 114.54s; chunk in 2.85s (22969 rows/s); X\u0000:110344724\n",
      "[read_vcf] 2555904 rows in 117.39s; chunk in 2.86s (22942 rows/s); X\u0000:113267805\n",
      "[read_vcf] 2621440 rows in 120.28s; chunk in 2.88s (22721 rows/s); X\u0000:116181346\n",
      "[read_vcf] 2686976 rows in 123.14s; chunk in 2.87s (22852 rows/s); X\u0000:118882692\n",
      "[read_vcf] 2752512 rows in 125.99s; chunk in 2.85s (23025 rows/s); X\u0000:121911362\n",
      "[read_vcf] 2818048 rows in 128.90s; chunk in 2.91s (22494 rows/s); X\u0000:124928450\n",
      "[read_vcf] 2883584 rows in 131.81s; chunk in 2.91s (22535 rows/s); X\u0000:127808992\n",
      "[read_vcf] 2949120 rows in 134.65s; chunk in 2.84s (23058 rows/s); X\u0000:130956155\n",
      "[read_vcf] 3014656 rows in 137.51s; chunk in 2.85s (22986 rows/s); X\u0000:134318747\n",
      "[read_vcf] 3080192 rows in 140.35s; chunk in 2.85s (23018 rows/s); X\u0000:137603907\n",
      "[read_vcf] 3145728 rows in 143.20s; chunk in 2.85s (23011 rows/s); X\u0000:140749939\n",
      "[read_vcf] 3211264 rows in 146.06s; chunk in 2.86s (22915 rows/s); X\u0000:143421964\n",
      "[read_vcf] 3276800 rows in 149.06s; chunk in 3.00s (21851 rows/s); X\u0000:146052307\n",
      "[read_vcf] 3342336 rows in 151.93s; chunk in 2.87s (22795 rows/s); X\u0000:149173075\n",
      "[read_vcf] 3407872 rows in 155.09s; chunk in 3.15s (20783 rows/s); X\u0000:152165183\n",
      "[read_vcf] 3468093 rows in 157.95s; chunk in 2.86s (21077 rows/s)\n",
      "[read_vcf] all done (21957 rows/s)\n"
     ]
    }
   ],
   "source": [
    "import zarr\n",
    "\n",
    "if \"vcf\" in dir():\n",
    "    del vcf\n",
    "\n",
    "vcf_filename = DATA + \"/ALL.chrX.phase3_shapeit2_mvncall_integrated_v1b.20130502.genotypes.vcf.gz\"\n",
    "\n",
    "vcf = allel.read_vcf(str(vcf_filename), [\"variants/ID\", \"variants/numalt\"],\n",
    "                     log=sys.stdout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numalt = vcf['variants/numalt']\n",
    "max_alt = np.max(numalt)\n",
    "max_alt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[vcf_to_zarr] 65536 rows in 18.37s; chunk in 18.37s (3567 rows/s); X\u0000:1734128\n",
      "[vcf_to_zarr] 131072 rows in 35.40s; chunk in 17.03s (3847 rows/s); X\u0000:3915040\n",
      "[vcf_to_zarr] 196608 rows in 50.55s; chunk in 15.15s (4326 rows/s); X\u0000:6474138\n",
      "[vcf_to_zarr] 262144 rows in 63.79s; chunk in 13.24s (4949 rows/s); X\u0000:9519142\n",
      "[vcf_to_zarr] 327680 rows in 76.67s; chunk in 12.88s (5086 rows/s); X\u0000:12551848\n",
      "[vcf_to_zarr] 393216 rows in 89.81s; chunk in 13.13s (4990 rows/s); X\u0000:15550284\n",
      "[vcf_to_zarr] 458752 rows in 102.19s; chunk in 12.39s (5290 rows/s); X\u0000:18776831\n",
      "[vcf_to_zarr] 524288 rows in 114.18s; chunk in 11.99s (5466 rows/s); X\u0000:21881750\n",
      "[vcf_to_zarr] 589824 rows in 126.26s; chunk in 12.08s (5426 rows/s); X\u0000:24574486\n",
      "[vcf_to_zarr] 655360 rows in 138.35s; chunk in 12.09s (5418 rows/s); X\u0000:27659891\n",
      "[vcf_to_zarr] 720896 rows in 150.55s; chunk in 12.20s (5370 rows/s); X\u0000:30671296\n",
      "[vcf_to_zarr] 786432 rows in 163.66s; chunk in 13.11s (4999 rows/s); X\u0000:33194562\n",
      "[vcf_to_zarr] 851968 rows in 177.75s; chunk in 14.09s (4652 rows/s); X\u0000:36031556\n",
      "[vcf_to_zarr] 917504 rows in 191.38s; chunk in 13.63s (4807 rows/s); X\u0000:39125534\n",
      "[vcf_to_zarr] 983040 rows in 204.88s; chunk in 13.49s (4856 rows/s); X\u0000:42186473\n",
      "[vcf_to_zarr] 1048576 rows in 219.50s; chunk in 14.62s (4481 rows/s); X\u0000:45309329\n",
      "[vcf_to_zarr] 1114112 rows in 232.82s; chunk in 13.33s (4917 rows/s); X\u0000:48359153\n",
      "[vcf_to_zarr] 1179648 rows in 245.89s; chunk in 13.07s (5014 rows/s); X\u0000:51881411\n",
      "[vcf_to_zarr] 1245184 rows in 259.15s; chunk in 13.26s (4943 rows/s); X\u0000:55705421\n",
      "[vcf_to_zarr] 1310720 rows in 273.22s; chunk in 14.06s (4659 rows/s); X\u0000:57996652\n",
      "[vcf_to_zarr] 1376256 rows in 286.04s; chunk in 12.83s (5108 rows/s); X\u0000:63082608\n",
      "[vcf_to_zarr] 1441792 rows in 298.62s; chunk in 12.58s (5210 rows/s); X\u0000:65358749\n",
      "[vcf_to_zarr] 1507328 rows in 311.25s; chunk in 12.62s (5191 rows/s); X\u0000:67634861\n",
      "[vcf_to_zarr] 1572864 rows in 323.59s; chunk in 12.34s (5309 rows/s); X\u0000:70712278\n",
      "[vcf_to_zarr] 1638400 rows in 336.21s; chunk in 12.62s (5191 rows/s); X\u0000:74091783\n",
      "[vcf_to_zarr] 1703936 rows in 348.67s; chunk in 12.46s (5261 rows/s); X\u0000:76280043\n",
      "[vcf_to_zarr] 1769472 rows in 361.05s; chunk in 12.38s (5293 rows/s); X\u0000:78926545\n",
      "[vcf_to_zarr] 1835008 rows in 373.81s; chunk in 12.76s (5136 rows/s); X\u0000:81422228\n",
      "[vcf_to_zarr] 1900544 rows in 387.15s; chunk in 13.34s (4912 rows/s); X\u0000:83965308\n",
      "[vcf_to_zarr] 1966080 rows in 399.99s; chunk in 12.84s (5103 rows/s); X\u0000:86441870\n",
      "[vcf_to_zarr] 2031616 rows in 412.45s; chunk in 12.46s (5259 rows/s); X\u0000:89012763\n",
      "[vcf_to_zarr] 2097152 rows in 425.30s; chunk in 12.85s (5100 rows/s); X\u0000:91967747\n",
      "[vcf_to_zarr] 2162688 rows in 437.85s; chunk in 12.55s (5220 rows/s); X\u0000:94617360\n",
      "[vcf_to_zarr] 2228224 rows in 450.25s; chunk in 12.40s (5284 rows/s); X\u0000:97514551\n",
      "[vcf_to_zarr] 2293760 rows in 463.18s; chunk in 12.93s (5068 rows/s); X\u0000:100546234\n",
      "[vcf_to_zarr] 2359296 rows in 475.96s; chunk in 12.77s (5130 rows/s); X\u0000:104064469\n",
      "[vcf_to_zarr] 2424832 rows in 488.93s; chunk in 12.97s (5052 rows/s); X\u0000:107375123\n",
      "[vcf_to_zarr] 2490368 rows in 501.72s; chunk in 12.79s (5124 rows/s); X\u0000:110344724\n",
      "[vcf_to_zarr] 2555904 rows in 514.26s; chunk in 12.54s (5224 rows/s); X\u0000:113267805\n",
      "[vcf_to_zarr] 2621440 rows in 526.71s; chunk in 12.44s (5266 rows/s); X\u0000:116181346\n",
      "[vcf_to_zarr] 2686976 rows in 539.17s; chunk in 12.46s (5258 rows/s); X\u0000:118882692\n",
      "[vcf_to_zarr] 2752512 rows in 551.53s; chunk in 12.37s (5300 rows/s); X\u0000:121911362\n",
      "[vcf_to_zarr] 2818048 rows in 563.95s; chunk in 12.41s (5279 rows/s); X\u0000:124928450\n",
      "[vcf_to_zarr] 2883584 rows in 576.40s; chunk in 12.45s (5263 rows/s); X\u0000:127808992\n",
      "[vcf_to_zarr] 2949120 rows in 588.82s; chunk in 12.42s (5275 rows/s); X\u0000:130956155\n",
      "[vcf_to_zarr] 3014656 rows in 601.21s; chunk in 12.39s (5288 rows/s); X\u0000:134318747\n",
      "[vcf_to_zarr] 3080192 rows in 613.58s; chunk in 12.37s (5299 rows/s); X\u0000:137603907\n",
      "[vcf_to_zarr] 3145728 rows in 625.91s; chunk in 12.33s (5313 rows/s); X\u0000:140749939\n",
      "[vcf_to_zarr] 3211264 rows in 639.38s; chunk in 13.47s (4866 rows/s); X\u0000:143421964\n",
      "[vcf_to_zarr] 3276800 rows in 653.17s; chunk in 13.78s (4754 rows/s); X\u0000:146052307\n",
      "[vcf_to_zarr] 3342336 rows in 667.36s; chunk in 14.20s (4616 rows/s); X\u0000:149173075\n",
      "[vcf_to_zarr] 3407872 rows in 682.42s; chunk in 15.06s (4352 rows/s); X\u0000:152165183\n",
      "[vcf_to_zarr] 3468093 rows in 695.06s; chunk in 12.64s (4765 rows/s)\n",
      "[vcf_to_zarr] all done (4985 rows/s)\n"
     ]
    }
   ],
   "source": [
    "zarr_file = vcf_filename.replace(\".vcf.gz\", \".zarr\")\n",
    "if not os.path.exists(zarr_file):\n",
    "    allel.vcf_to_zarr(str(vcf_filename),\n",
    "                      str(zarr_file),\n",
    "                      fields='*', alt_number=max_alt,\n",
    "                      log=sys.stdout, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the VCF file has been converted to Zarr file, I can now create a dataframe extracting the rows that are pathogenic variants of prostate cancer from 1000 genome project Zarr file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "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>CHROM</th>\n",
       "      <th>POS</th>\n",
       "      <th>REF</th>\n",
       "      <th>ALT</th>\n",
       "      <th>GENEINFO</th>\n",
       "      <th>CLNDN</th>\n",
       "      <th>CLNSIG</th>\n",
       "      <th>AF_TGP</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>X</td>\n",
       "      <td>66765802</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate|not_provided</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>X</td>\n",
       "      <td>66766163</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>X</td>\n",
       "      <td>66937326</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate|Prostate_cancer_su...</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>0.00053</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>X</td>\n",
       "      <td>66937337</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>AR:367</td>\n",
       "      <td>Malignant_tumor_of_prostate|Prostate_cancer</td>\n",
       "      <td>Pathogenic</td>\n",
       "      <td>0.00053</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  CHROM       POS REF ALT GENEINFO  \\\n",
       "0     X  66765802   C   T   AR:367   \n",
       "1     X  66766163   C   G   AR:367   \n",
       "2     X  66937326   G   T   AR:367   \n",
       "3     X  66937337   G   A   AR:367   \n",
       "\n",
       "                                               CLNDN      CLNSIG   AF_TGP  \n",
       "0           Malignant_tumor_of_prostate|not_provided  Pathogenic      NaN  \n",
       "1                        Malignant_tumor_of_prostate  Pathogenic      NaN  \n",
       "2  Malignant_tumor_of_prostate|Prostate_cancer_su...  Pathogenic  0.00053  \n",
       "3        Malignant_tumor_of_prostate|Prostate_cancer  Pathogenic  0.00053  "
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vcfzarr = zarr.open_group(str(zarr_file), mode=\"r\")\n",
    "\n",
    "pos_idx = np.where(np.isin(vcfzarr[\"variants/POS\"], pathogenic_df.POS))[0]\n",
    "\n",
    "vdf =  pd.DataFrame(\n",
    "    dict([(k, vcfzarr[\"variants/\" + k][pos_idx]) \n",
    "    for k in [\"CHROM\", \"POS\", \"REF\"]])\n",
    ")\n",
    "vdf[\"ALT\"] = vcfzarr[\"variants/ALT\"][pos_idx, 0]\n",
    "\n",
    "vdf = vdf.merge(pathogenic_df)\n",
    "\n",
    "vdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a total of 4 variants observed in 1000 genome project. Now that we have an idea of which variants appear in 1000 genome project, we can extract subjects with non-homozygous reference genotype to identify subjects carrying pathogenic variants of prostate cancer. In the genotype array, we would be essentially looking for subjects with genotypes that are not 0/0 for female or 0/. for male."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "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>HG00321</th>\n",
       "      <th>HG00371</th>\n",
       "      <th>HG01850</th>\n",
       "      <th>HG02031</th>\n",
       "      <th>HG02069</th>\n",
       "      <th>HG02787</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>X:66765802</th>\n",
       "      <td>b'0/.'</td>\n",
       "      <td>b'0/.'</td>\n",
       "      <td>b'0/0'</td>\n",
       "      <td>b'0/0'</td>\n",
       "      <td>b'0/0'</td>\n",
       "      <td>b'0/1'</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>X:66766163</th>\n",
       "      <td>b'0/.'</td>\n",
       "      <td>b'0/.'</td>\n",
       "      <td>b'0/0'</td>\n",
       "      <td>b'1/0'</td>\n",
       "      <td>b'0/0'</td>\n",
       "      <td>b'0/0'</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>X:66937326</th>\n",
       "      <td>b'1/.'</td>\n",
       "      <td>b'1/.'</td>\n",
       "      <td>b'0/0'</td>\n",
       "      <td>b'0/0'</td>\n",
       "      <td>b'0/0'</td>\n",
       "      <td>b'0/0'</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>X:66937337</th>\n",
       "      <td>b'0/.'</td>\n",
       "      <td>b'0/.'</td>\n",
       "      <td>b'0/1'</td>\n",
       "      <td>b'0/0'</td>\n",
       "      <td>b'0/1'</td>\n",
       "      <td>b'0/0'</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           HG00321 HG00371 HG01850 HG02031 HG02069 HG02787\n",
       "X:66765802  b'0/.'  b'0/.'  b'0/0'  b'0/0'  b'0/0'  b'0/1'\n",
       "X:66766163  b'0/.'  b'0/.'  b'0/0'  b'1/0'  b'0/0'  b'0/0'\n",
       "X:66937326  b'1/.'  b'1/.'  b'0/0'  b'0/0'  b'0/0'  b'0/0'\n",
       "X:66937337  b'0/.'  b'0/.'  b'0/1'  b'0/0'  b'0/1'  b'0/0'"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# reset idx\n",
    "pos_idx = np.where(np.isin(vcfzarr[\"variants/POS\"], vdf.POS))[0]\n",
    "\n",
    "gt = allel.GenotypeArray(np.array([vcfzarr[\"calldata/GT\"][i] for i in pos_idx]))\n",
    "\n",
    "# convert genotype info to pandas dataframe\n",
    "gt_df = pd.DataFrame(\n",
    "    gt.to_gt(),\n",
    "    columns=np.array(vcfzarr[\"samples\"]),\n",
    "    index=\"X:\" + pd.Series(vcfzarr[\"variants/POS\"]).loc[pos_idx].astype(str)\n",
    ")\n",
    "\n",
    "bool_df = (gt_df != b'0/.') & (gt_df != b'0/0')\n",
    "bool_df.columns.values[bool_df.any(axis=0)]\n",
    "\n",
    "gt_df[bool_df.columns.values[bool_df.any(axis=0)]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 6 subjects that have variants of interest. Now we can create a new dataframe summarizing genotype information from the provided csv file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = gt_df.loc[:, bool_df.any(axis=0)].melt(ignore_index=False, var_name=\"subject_id\", value_name=\"genotype\")\n",
    "\n",
    "# remove reference entries\n",
    "df = df[(df.genotype != b\"0/0\") & (df.genotype != b\"0/.\")].reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge population table to get super population code\n",
    "meta = pd.read_csv(DATA + \"/1kgn_phase3_population.csv\")\n",
    "out = df.merge(meta[[\"Sample name\", \"Superpopulation code\", \"Sex\"]], left_on=\"subject_id\", right_on=\"Sample name\")\n",
    "out = out.rename(columns={\n",
    "    \"index\": \"variant_of_interest\",\n",
    "    \"Sex\": \"gender\",\n",
    "    \"Superpopulation code\": \"superpopulation_code\",\n",
    "}).drop(\"Sample name\", axis=1)[[\"subject_id\", \"gender\", \"superpopulation_code\", \"variant_of_interest\", \"genotype\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "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>subject_id</th>\n",
       "      <th>gender</th>\n",
       "      <th>superpopulation_code</th>\n",
       "      <th>variant_of_interest</th>\n",
       "      <th>genotype</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>HG00321</td>\n",
       "      <td>male</td>\n",
       "      <td>EUR</td>\n",
       "      <td>X:66937326</td>\n",
       "      <td>b'1/.'</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>HG00371</td>\n",
       "      <td>male</td>\n",
       "      <td>EUR</td>\n",
       "      <td>X:66937326</td>\n",
       "      <td>b'1/.'</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>HG01850</td>\n",
       "      <td>female</td>\n",
       "      <td>EAS</td>\n",
       "      <td>X:66937337</td>\n",
       "      <td>b'0/1'</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>HG02031</td>\n",
       "      <td>female</td>\n",
       "      <td>EAS</td>\n",
       "      <td>X:66766163</td>\n",
       "      <td>b'1/0'</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>HG02069</td>\n",
       "      <td>female</td>\n",
       "      <td>EAS</td>\n",
       "      <td>X:66937337</td>\n",
       "      <td>b'0/1'</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>HG02787</td>\n",
       "      <td>female</td>\n",
       "      <td>SAS</td>\n",
       "      <td>X:66765802</td>\n",
       "      <td>b'0/1'</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  subject_id  gender superpopulation_code variant_of_interest genotype\n",
       "0    HG00321    male                  EUR          X:66937326   b'1/.'\n",
       "1    HG00371    male                  EUR          X:66937326   b'1/.'\n",
       "2    HG01850  female                  EAS          X:66937337   b'0/1'\n",
       "3    HG02031  female                  EAS          X:66766163   b'1/0'\n",
       "4    HG02069  female                  EAS          X:66937337   b'0/1'\n",
       "5    HG02787  female                  SAS          X:66765802   b'0/1'"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the table above, we see that two male subjects had pathogenic variants of prostate cancer. The two subjects are from European subpopulation."
   ]
  },
  {
   "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"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}