PRS-Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PRS Tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load in base data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "# from: https://www.med.unc.edu/pgc/download-results/\n",
    "filepath = '../files/daner_natgen_pgc_eas' # Dataset: scz2019asi, PMID: 31740837\n",
    "base_df = pd.read_table(filepath, sep=' ') # load in file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filter base data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10,694,924 SNPs in complete dataframe\n",
      "-9,947,346 SNPs are out of the (None, 0.05) p-value bounds (747,578 SNPs remain)\n",
      "-111,352 of the remaining SNPs are strand-ambiguous (636,226 SNPs remain)\n",
      "-35,113 of the remaining SNPs are indels (601,113 SNPs remain)\n",
      "-90,390 of the remaining SNPs have imputation scores lower than 0.8\n",
      "\n",
      "510,723 SNPs remain\n"
     ]
    }
   ],
   "source": [
    "import basefilter\n",
    "\n",
    "BASE_DATA = basefilter.Filter(base_df, \\\n",
    "                pval_range=(None, 0.05), \\\n",
    "                min_imputation_info_score=0.8, \\\n",
    "                    # https://choishingwan.github.io/PRS-Tutorial/base/#standard-gwas-qc\n",
    "                    # https://www.nature.com/articles/s41596-020-0353-1\n",
    "                remove_strand_ambig_snps=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preview base data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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>CHR</th>\n",
       "      <th>BP</th>\n",
       "      <th>A1</th>\n",
       "      <th>A2</th>\n",
       "      <th>FRQ_A_22778</th>\n",
       "      <th>FRQ_U_35362</th>\n",
       "      <th>INFO</th>\n",
       "      <th>OR</th>\n",
       "      <th>SE</th>\n",
       "      <th>P</th>\n",
       "      <th>ngt</th>\n",
       "      <th>HetISqt</th>\n",
       "      <th>HetDf</th>\n",
       "      <th>HetPVa</th>\n",
       "      <th>Nca</th>\n",
       "      <th>Nco</th>\n",
       "      <th>Neff</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SNP</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>rs190028338</th>\n",
       "      <td>1</td>\n",
       "      <td>100034461</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>0.9880</td>\n",
       "      <td>0.9810</td>\n",
       "      <td>0.954</td>\n",
       "      <td>1.62434</td>\n",
       "      <td>0.2165</td>\n",
       "      <td>0.025050</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>1477.0</td>\n",
       "      <td>1680.0</td>\n",
       "      <td>1563.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rs144918958</th>\n",
       "      <td>1</td>\n",
       "      <td>100072586</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>0.9740</td>\n",
       "      <td>0.9890</td>\n",
       "      <td>0.892</td>\n",
       "      <td>0.34249</td>\n",
       "      <td>0.3817</td>\n",
       "      <td>0.004997</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>547.0</td>\n",
       "      <td>540.0</td>\n",
       "      <td>543.48</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rs12136177</th>\n",
       "      <td>1</td>\n",
       "      <td>100149952</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>0.0654</td>\n",
       "      <td>0.0695</td>\n",
       "      <td>0.984</td>\n",
       "      <td>0.93128</td>\n",
       "      <td>0.0336</td>\n",
       "      <td>0.033930</td>\n",
       "      <td>6</td>\n",
       "      <td>45.3</td>\n",
       "      <td>3.0</td>\n",
       "      <td>0.1398</td>\n",
       "      <td>14583.0</td>\n",
       "      <td>17261.0</td>\n",
       "      <td>14916.61</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rs6664416</th>\n",
       "      <td>1</td>\n",
       "      <td>100162963</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>0.0657</td>\n",
       "      <td>0.0698</td>\n",
       "      <td>0.989</td>\n",
       "      <td>0.93473</td>\n",
       "      <td>0.0334</td>\n",
       "      <td>0.043730</td>\n",
       "      <td>3</td>\n",
       "      <td>41.8</td>\n",
       "      <td>3.0</td>\n",
       "      <td>0.1608</td>\n",
       "      <td>14583.0</td>\n",
       "      <td>17261.0</td>\n",
       "      <td>14916.61</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rs6686520</th>\n",
       "      <td>1</td>\n",
       "      <td>100167505</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>0.0656</td>\n",
       "      <td>0.0698</td>\n",
       "      <td>0.986</td>\n",
       "      <td>0.93379</td>\n",
       "      <td>0.0335</td>\n",
       "      <td>0.040620</td>\n",
       "      <td>1</td>\n",
       "      <td>38.2</td>\n",
       "      <td>3.0</td>\n",
       "      <td>0.1831</td>\n",
       "      <td>14583.0</td>\n",
       "      <td>17261.0</td>\n",
       "      <td>14916.61</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             CHR         BP A1 A2  FRQ_A_22778  FRQ_U_35362   INFO       OR  \\\n",
       "SNP                                                                           \n",
       "rs190028338    1  100034461  A  C       0.9880       0.9810  0.954  1.62434   \n",
       "rs144918958    1  100072586  T  C       0.9740       0.9890  0.892  0.34249   \n",
       "rs12136177     1  100149952  T  C       0.0654       0.0695  0.984  0.93128   \n",
       "rs6664416      1  100162963  A  C       0.0657       0.0698  0.989  0.93473   \n",
       "rs6686520      1  100167505  A  G       0.0656       0.0698  0.986  0.93379   \n",
       "\n",
       "                 SE         P  ngt  HetISqt  HetDf  HetPVa      Nca      Nco  \\\n",
       "SNP                                                                            \n",
       "rs190028338  0.2165  0.025050    0      0.0    0.0  1.0000   1477.0   1680.0   \n",
       "rs144918958  0.3817  0.004997    0      0.0    0.0  1.0000    547.0    540.0   \n",
       "rs12136177   0.0336  0.033930    6     45.3    3.0  0.1398  14583.0  17261.0   \n",
       "rs6664416    0.0334  0.043730    3     41.8    3.0  0.1608  14583.0  17261.0   \n",
       "rs6686520    0.0335  0.040620    1     38.2    3.0  0.1831  14583.0  17261.0   \n",
       "\n",
       "                 Neff  \n",
       "SNP                    \n",
       "rs190028338   1563.00  \n",
       "rs144918958    543.48  \n",
       "rs12136177   14916.61  \n",
       "rs6664416    14916.61  \n",
       "rs6686520    14916.61  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BASE_DATA.FRAME.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up ScoreEngine with base data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scoring\n",
    "\n",
    "engine = scoring.ScoreEngine(BASE_DATA.DICT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate risk score at one SNP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"./img/prs_eq.png\" width=\"800\" height=\"150\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(url= \"./img/prs_eq.png\", width=800, height=150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.4468121715467305"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "engine.score_snp('rs189548203', 'AA')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rs189548203 (chr14:70775453)\n",
      "Discov alleles:\tTG (T=risk)\n",
      "Target alleles:\tTT (reversed? True)\n",
      "Risk allele count: 2\n",
      "Odds ratio: 0.8 => ln odds ratio: -0.223\n",
      "Risk score: 2 * -0.223 = -0.4468\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-0.4468121715467305"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "engine.score_snp('rs189548203', 'AA', print=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate PRS on random fictional sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "alleles = ['A','G','T','C']\n",
    "\n",
    "def random_index():\n",
    "    \"\"\" Generate a random integer from 0 to 3. \"\"\"\n",
    "    return random.randint(0, len(alleles)-1)\n",
    "\n",
    "def rand_alleles():\n",
    "    \"\"\" Put together two random alleles from the list of possibilities. \"\"\"\n",
    "    return alleles[random_index()] + alleles[random_index()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generator of tuples composed of SNPs from base data and random alleles\n",
    "count = len(BASE_DATA.DICT.keys())\n",
    "pairings = zip(BASE_DATA.DICT.keys(), [rand_alleles() for _ in range(count)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate the polygenic risk score of all SNP & allele pairings in `pairings` generator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 510723/510723 [00:02<00:00, 184101.24it/s]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "620.5582032666241"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "engine.score_polygenic_risk(pairings, generator_count=count)"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
  },
  "kernelspec": {
   "display_name": "Python 3.7.9 64-bit",
   "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.7.9"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}