mutation_enrichment.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/peter/opt/anaconda2/lib/python2.7/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n",
      "  from pandas.core import datetools\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from __future__ import division\n",
    "import re\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import numpy.random as npr\n",
    "import glob\n",
    "import scipy.stats as st\n",
    "import statsmodels.api as sm\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary\n",
    "\n",
    "This notebook produces the results on mutation enrichment presented in Table S2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "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>fqid</th>\n",
       "      <th>position_rcrs</th>\n",
       "      <th>A</th>\n",
       "      <th>C</th>\n",
       "      <th>G</th>\n",
       "      <th>T</th>\n",
       "      <th>a</th>\n",
       "      <th>c</th>\n",
       "      <th>g</th>\n",
       "      <th>t</th>\n",
       "      <th>...</th>\n",
       "      <th>Sex</th>\n",
       "      <th>spike_in</th>\n",
       "      <th>haplogroup</th>\n",
       "      <th>fam_str</th>\n",
       "      <th>fam_cat</th>\n",
       "      <th>mot_cat</th>\n",
       "      <th>age_collection</th>\n",
       "      <th>age_birth</th>\n",
       "      <th>twin_info</th>\n",
       "      <th>heteroplasmy</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>M132-bl</td>\n",
       "      <td>14461</td>\n",
       "      <td>0</td>\n",
       "      <td>181</td>\n",
       "      <td>0</td>\n",
       "      <td>3895</td>\n",
       "      <td>0</td>\n",
       "      <td>170</td>\n",
       "      <td>0</td>\n",
       "      <td>4170</td>\n",
       "      <td>...</td>\n",
       "      <td>Female</td>\n",
       "      <td>0</td>\n",
       "      <td>J1c5d</td>\n",
       "      <td>0-0-1-2</td>\n",
       "      <td>m1c2</td>\n",
       "      <td>m1c2</td>\n",
       "      <td>16658</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>yes</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>M132-ch</td>\n",
       "      <td>14461</td>\n",
       "      <td>0</td>\n",
       "      <td>170</td>\n",
       "      <td>0</td>\n",
       "      <td>3327</td>\n",
       "      <td>0</td>\n",
       "      <td>147</td>\n",
       "      <td>0</td>\n",
       "      <td>3694</td>\n",
       "      <td>...</td>\n",
       "      <td>Female</td>\n",
       "      <td>0</td>\n",
       "      <td>J1c5d</td>\n",
       "      <td>0-0-1-2</td>\n",
       "      <td>m1c2</td>\n",
       "      <td>m1c2</td>\n",
       "      <td>16658</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>yes</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2 rows × 29 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      fqid  position_rcrs  A    C  G     T  a    c  g     t     ...       \\\n",
       "0  M132-bl          14461  0  181  0  3895  0  170  0  4170     ...        \n",
       "1  M132-ch          14461  0  170  0  3327  0  147  0  3694     ...        \n",
       "\n",
       "      Sex spike_in haplogroup  fam_str fam_cat mot_cat age_collection  \\\n",
       "0  Female        0      J1c5d  0-0-1-2    m1c2    m1c2          16658   \n",
       "1  Female        0      J1c5d  0-0-1-2    m1c2    m1c2          16658   \n",
       "\n",
       "  age_birth twin_info heteroplasmy  \n",
       "0       NaN       NaN          yes  \n",
       "1       NaN       NaN          yes  \n",
       "\n",
       "[2 rows x 29 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load in the heteroplasmy data.\n",
    "counts = pd.read_csv('table_s10_famfreqs.tsv', sep='\\t', skiprows=1)\n",
    "counts.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Load in the annotations\n",
    "annot = pd.read_csv('table_s1_heteroplasmy_data.tsv', sep='\\t', skiprows=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# sizes of the mt-genome in different categories\n",
    "d_loop_len = 1122\n",
    "trna_len = 1508\n",
    "rrna_len = 2513\n",
    "protein_s_len = 2834\n",
    "protein_n_len = 8533\n",
    "intergenic_len = 88\n",
    "all_lens = np.array([d_loop_len, trna_len, rrna_len, protein_s_len, protein_n_len, intergenic_len])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# group by family-position combination\n",
    "mut_annot = annot.groupby(['family_id', 'position_rcrs']).apply(lambda x: x.iloc[0,:])\n",
    "mut_annot['gene'].replace(np.nan, 'intergenic', inplace=True)\n",
    "tot_obs = mut_annot.shape[0]\n",
    "assert tot_obs == 346"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2377,  4332,  1282,  4520, 11088])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the mutation positions\n",
    "positions = mut_annot.index.get_level_values(1).values\n",
    "positions[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### D-loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "is_d_loop = lambda x: (16024 <= x) | (x <= 576)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "73"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_d_loop = np.array(map(is_d_loop, positions)).sum()\n",
    "num_d_loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D-loop: observed, expected, binomial 2-sided p-value\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(73, 23.43001991671193, 4.0500484611955255e-18)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from __future__ import division\n",
    "p_d_loop_expected = is_d_loop(np.arange(1,16570)).sum() / 16569\n",
    "pval = st.binom_test(num_d_loop, tot_obs, p=p_d_loop_expected, alternative='two-sided')\n",
    "print 'D-loop: observed, expected, binomial 2-sided p-value'\n",
    "is_d_loop(positions).sum(), p_d_loop_expected*346, pval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tRNA: observed, expected, binomial 2-sided p-value\n",
      "36 31.4906150039 0.3995571716659997\n"
     ]
    }
   ],
   "source": [
    "obs_trna = mut_annot['gene'].str.startswith('TR').sum()\n",
    "p_trna_exp = 1508.0/16569\n",
    "pval = st.binom_test(obs_trna, tot_obs, p_trna_exp, alternative='two-sided')\n",
    "print 'tRNA: observed, expected, binomial 2-sided p-value'\n",
    "print obs_trna, p_trna_exp*tot_obs, pval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rRNA: observed, expected, binomial 2-sided p-value\n",
      "42 52.4773975496 0.0606555594659883\n"
     ]
    }
   ],
   "source": [
    "obs_rrna = mut_annot['gene'].str.startswith('RN').sum()\n",
    "p_rrna_exp = 2513.0/16569\n",
    "pval = st.binom_test(obs_rrna, tot_obs, p_trna_exp, alternative='two-sided')\n",
    "print 'rRNA: observed, expected, binomial 2-sided p-value'\n",
    "print obs_rrna, p_rrna_exp*tot_obs, pval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "syn: observed, expected, binomial 2-sided p-value\n",
      "76 59.1806385419 0.018445595490323933\n"
     ]
    }
   ],
   "source": [
    "num_syn = (mut_annot['syn'] == 'syn').sum()\n",
    "p_syn_exp = 2834.0/16569\n",
    "pval = st.binom_test(num_syn, tot_obs, p_syn_exp)\n",
    "print 'syn: observed, expected, binomial 2-sided p-value'\n",
    "print num_syn, p_syn_exp*tot_obs, pval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nonsyn: observed, expected, binomial 2-sided p-value\n",
      "118 178.189269117 1.0142234592112752e-10\n"
     ]
    }
   ],
   "source": [
    "num_nonsyn = (mut_annot['syn'] == 'nsyn').sum()\n",
    "p_nsyn_exp = 8533/16569\n",
    "pval = st.binom_test(num_nonsyn, tot_obs, p_nsyn_exp)\n",
    "print 'nonsyn: observed, expected, binomial 2-sided p-value'\n",
    "print num_nonsyn, p_nsyn_exp*tot_obs, pval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "intergenic: observed, expected, binomial 2-sided p-value\n",
      "1 1.83764862092 1.0\n"
     ]
    }
   ],
   "source": [
    "num_intergenic = (mut_annot['gene'] == 'intergenic').sum()  # == 1\n",
    "p_inter_exp = 88.0/16569\n",
    "pval = st.binom_test(num_intergenic, tot_obs, p_inter_exp)\n",
    "print 'intergenic: observed, expected, binomial 2-sided p-value'\n",
    "print num_intergenic, p_inter_exp*tot_obs, pval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "is_coding = mut_annot['gene'].isin(['ATP6', 'ATP8', 'COX1', 'COX2', 'COX3', 'CYTB', 'ND1', 'ND2', 'ND3', 'ND4', 'ND5', 'ND6'])\n",
    "num_coding_observed = is_coding.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This parses the features for genes.\n",
    "# I downloaded features from NCBI GenBank rCRS entry:\n",
    "# https://www.ncbi.nlm.nih.gov/nuccore/251831106\n",
    "with open('features_rcrs.txt') as fin:\n",
    "    features = np.array([re.sub('[A-Za-z()]', '', re.split(' +', el.strip())[-1]).split('..') for el in fin if 'CDS' in el]).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3307,  4262],\n",
       "       [ 4470,  5511],\n",
       "       [ 5904,  7445],\n",
       "       [ 7586,  8269],\n",
       "       [ 8366,  8572],\n",
       "       [ 8527,  9207],\n",
       "       [ 9207,  9990],\n",
       "       [10059, 10404],\n",
       "       [10470, 10766],\n",
       "       [10760, 12137],\n",
       "       [12337, 14148],\n",
       "       [14149, 14673],\n",
       "       [14747, 15887]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(194, 237.68314321926488, 9.595573572388399e-07)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_coding_sites = np.diff(features, axis=1).sum()\n",
    "p_coding_exp = num_coding_sites / 16569.0\n",
    "num_coding_sites, p_coding_exp\n",
    "num_coding_observed, p_coding_exp*346, st.binom_test(num_coding_observed, 346, p_coding_exp)"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python [anaconda2]",
   "language": "python",
   "name": "Python [anaconda2]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}