weeks/week09/Module 1 Exercise 1 - Population Divergence.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Exercise 1: Read in population allele frequencies generated by Plink\n",
    "# Use PLINK to obtain allele frequencies for each of the four HapMap Populations.  \n",
    "# Files with family and individual IDs for the populations:\n",
    "# hapmap-CEU.txt, hapmap-YRI.txt, hapmap-CHB.txt, and hapmap-ASW.txt.  \n",
    "# You can use the following PLINK command to obtain allele frequencies for only the CEU individuals.\n",
    "# ./plink --bfile  plink_hapmap_CEU_YRI_ASW_CHB –-keep hapmap-CEU.txt --freq --out CEU\n",
    "CEU.data=read.table(\"~/Dropbox/Teaching/2015/BMS262/Module1/Excercises/CEU.frq\",header=TRUE)\n",
    "YRI.data=read.table(\"~/Dropbox/Teaching/2015/BMS262/Module1/Excercises/YRI.frq\",header=TRUE)\n",
    "ASW.data=read.table(\"~/Dropbox/Teaching/2015/BMS262/Module1/Excercises/ASW.frq\",header=TRUE)\n",
    "CHB.data=read.table(\"~/Dropbox/Teaching/2015/BMS262/Module1/Excercises/CHB.frq\",header=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>CHR</th><th scope=col>SNP</th><th scope=col>A1</th><th scope=col>A2</th><th scope=col>MAF</th><th scope=col>NCHROBS</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>1</td><td>rs9442372</td><td>A</td><td>G</td><td>0.4107</td><td>224</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>1</td><td>rs2887286</td><td>C</td><td>T</td><td>0.1607</td><td>224</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>1</td><td>rs3813199</td><td>A</td><td>G</td><td>0.09375</td><td>224</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>1</td><td>rs6685064</td><td>T</td><td>C</td><td>0.07143</td><td>224</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>1</td><td>rs9439462</td><td>T</td><td>C</td><td>0.04505</td><td>222</td></tr>\n",
       "\t<tr><th scope=row>6</th><td>1</td><td>rs3820011</td><td>A</td><td>C</td><td>0.3063</td><td>222</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llllll}\n",
       "  & CHR & SNP & A1 & A2 & MAF & NCHROBS\\\\\n",
       "\\hline\n",
       "\t1 & 1 & rs9442372 & A & G & 0.4107 & 224\\\\\n",
       "\t2 & 1 & rs2887286 & C & T & 0.1607 & 224\\\\\n",
       "\t3 & 1 & rs3813199 & A & G & 0.09375 & 224\\\\\n",
       "\t4 & 1 & rs6685064 & T & C & 0.07143 & 224\\\\\n",
       "\t5 & 1 & rs9439462 & T & C & 0.04505 & 222\\\\\n",
       "\t6 & 1 & rs3820011 & A & C & 0.3063 & 222\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "  CHR       SNP A1 A2     MAF NCHROBS\n",
       "1   1 rs9442372  A  G 0.41070     224\n",
       "2   1 rs2887286  C  T 0.16070     224\n",
       "3   1 rs3813199  A  G 0.09375     224\n",
       "4   1 rs6685064  T  C 0.07143     224\n",
       "5   1 rs9439462  T  C 0.04505     222\n",
       "6   1 rs3820011  A  C 0.30630     222"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>CHR</th><th scope=col>SNP</th><th scope=col>A1</th><th scope=col>A2</th><th scope=col>MAF</th><th scope=col>NCHROBS</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>1</td><td>rs9442372</td><td>A</td><td>G</td><td>0.4898</td><td>294</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>1</td><td>rs2887286</td><td>T</td><td>C</td><td>0.4864</td><td>294</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>1</td><td>rs3813199</td><td>A</td><td>G</td><td>0.1837</td><td>294</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>1</td><td>rs6685064</td><td>C</td><td>T</td><td>0.4048</td><td>294</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>1</td><td>rs9439462</td><td>C</td><td>T</td><td>0.2483</td><td>294</td></tr>\n",
       "\t<tr><th scope=row>6</th><td>1</td><td>rs3820011</td><td>A</td><td>C</td><td>0.1497</td><td>294</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llllll}\n",
       "  & CHR & SNP & A1 & A2 & MAF & NCHROBS\\\\\n",
       "\\hline\n",
       "\t1 & 1 & rs9442372 & A & G & 0.4898 & 294\\\\\n",
       "\t2 & 1 & rs2887286 & T & C & 0.4864 & 294\\\\\n",
       "\t3 & 1 & rs3813199 & A & G & 0.1837 & 294\\\\\n",
       "\t4 & 1 & rs6685064 & C & T & 0.4048 & 294\\\\\n",
       "\t5 & 1 & rs9439462 & C & T & 0.2483 & 294\\\\\n",
       "\t6 & 1 & rs3820011 & A & C & 0.1497 & 294\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "  CHR       SNP A1 A2    MAF NCHROBS\n",
       "1   1 rs9442372  A  G 0.4898     294\n",
       "2   1 rs2887286  T  C 0.4864     294\n",
       "3   1 rs3813199  A  G 0.1837     294\n",
       "4   1 rs6685064  C  T 0.4048     294\n",
       "5   1 rs9439462  C  T 0.2483     294\n",
       "6   1 rs3820011  A  C 0.1497     294"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "head(CEU.data)\n",
    "head(YRI.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>CHR</th><th scope=col>SNP</th><th scope=col>A1</th><th scope=col>A2</th><th scope=col>MAF</th><th scope=col>NCHROBS</th><th scope=col>NewFreq</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>1</td><td>rs9442372</td><td>A</td><td>G</td><td>0.4898</td><td>294</td><td>0.4898</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>1</td><td>rs2887286</td><td>T</td><td>C</td><td>0.4864</td><td>294</td><td>0.5136</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>1</td><td>rs3813199</td><td>A</td><td>G</td><td>0.1837</td><td>294</td><td>0.1837</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>1</td><td>rs6685064</td><td>C</td><td>T</td><td>0.4048</td><td>294</td><td>0.5952</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>1</td><td>rs9439462</td><td>C</td><td>T</td><td>0.2483</td><td>294</td><td>0.7517</td></tr>\n",
       "\t<tr><th scope=row>6</th><td>1</td><td>rs3820011</td><td>A</td><td>C</td><td>0.1497</td><td>294</td><td>0.1497</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllllll}\n",
       "  & CHR & SNP & A1 & A2 & MAF & NCHROBS & NewFreq\\\\\n",
       "\\hline\n",
       "\t1 & 1 & rs9442372 & A & G & 0.4898 & 294 & 0.4898\\\\\n",
       "\t2 & 1 & rs2887286 & T & C & 0.4864 & 294 & 0.5136\\\\\n",
       "\t3 & 1 & rs3813199 & A & G & 0.1837 & 294 & 0.1837\\\\\n",
       "\t4 & 1 & rs6685064 & C & T & 0.4048 & 294 & 0.5952\\\\\n",
       "\t5 & 1 & rs9439462 & C & T & 0.2483 & 294 & 0.7517\\\\\n",
       "\t6 & 1 & rs3820011 & A & C & 0.1497 & 294 & 0.1497\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "  CHR       SNP A1 A2    MAF NCHROBS NewFreq\n",
       "1   1 rs9442372  A  G 0.4898     294  0.4898\n",
       "2   1 rs2887286  T  C 0.4864     294  0.5136\n",
       "3   1 rs3813199  A  G 0.1837     294  0.1837\n",
       "4   1 rs6685064  C  T 0.4048     294  0.5952\n",
       "5   1 rs9439462  C  T 0.2483     294  0.7517\n",
       "6   1 rs3820011  A  C 0.1497     294  0.1497"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#### Use the same reference allele ####\n",
    "CEU.frq=CEU.data$MAF\n",
    "SWITCH1=CEU.data$A1!=YRI.data$A1\n",
    "NewFreq=YRI.data$MAF\n",
    "NewFreq[SWITCH1]=1-YRI.data$MAF[SWITCH1]\n",
    "YRI.data$NewFreq=NewFreq\n",
    "YRI.frq=NewFreq\n",
    "head(YRI.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>0.4107</li>\n",
       "\t<li>0.1607</li>\n",
       "\t<li>0.09375</li>\n",
       "\t<li>0.07143</li>\n",
       "\t<li>0.04505</li>\n",
       "\t<li>0.3063</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0.4107\n",
       "\\item 0.1607\n",
       "\\item 0.09375\n",
       "\\item 0.07143\n",
       "\\item 0.04505\n",
       "\\item 0.3063\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0.4107\n",
       "2. 0.1607\n",
       "3. 0.09375\n",
       "4. 0.07143\n",
       "5. 0.04505\n",
       "6. 0.3063\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 0.41070 0.16070 0.09375 0.07143 0.04505 0.30630"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>0.4898</li>\n",
       "\t<li>0.5136</li>\n",
       "\t<li>0.1837</li>\n",
       "\t<li>0.5952</li>\n",
       "\t<li>0.7517</li>\n",
       "\t<li>0.1497</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0.4898\n",
       "\\item 0.5136\n",
       "\\item 0.1837\n",
       "\\item 0.5952\n",
       "\\item 0.7517\n",
       "\\item 0.1497\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0.4898\n",
       "2. 0.5136\n",
       "3. 0.1837\n",
       "4. 0.5952\n",
       "5. 0.7517\n",
       "6. 0.1497\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 0.4898 0.5136 0.1837 0.5952 0.7517 0.1497"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>0.45025</li>\n",
       "\t<li>0.33715</li>\n",
       "\t<li>0.138725</li>\n",
       "\t<li>0.333315</li>\n",
       "\t<li>0.398375</li>\n",
       "\t<li>0.228</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0.45025\n",
       "\\item 0.33715\n",
       "\\item 0.138725\n",
       "\\item 0.333315\n",
       "\\item 0.398375\n",
       "\\item 0.228\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0.45025\n",
       "2. 0.33715\n",
       "3. 0.138725\n",
       "4. 0.333315\n",
       "5. 0.398375\n",
       "6. 0.228\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 0.450250 0.337150 0.138725 0.333315 0.398375 0.228000"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Exercise 2: Obtain Fst measures for CEU versus YRI Hapmap Populations\n",
    "FREQpop1=CEU.frq\n",
    "FREQpop2=YRI.frq\n",
    "FREQ=(FREQpop1+FREQpop2)/2\n",
    "\n",
    "INCLUDE=(FREQ>0 & FREQ<1)\n",
    "FREQpop1=FREQpop1[INCLUDE]\n",
    "FREQpop2=FREQpop2[INCLUDE]\n",
    "FREQ=FREQ[INCLUDE]\n",
    "\n",
    "head(FREQpop1)\n",
    "head(FREQpop2)\n",
    "head(FREQ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.144598943242045"
      ],
      "text/latex": [
       "0.144598943242045"
      ],
      "text/markdown": [
       "0.144598943242045"
      ],
      "text/plain": [
       "[1] 0.1445989"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "0.147010883351759"
      ],
      "text/latex": [
       "0.147010883351759"
      ],
      "text/markdown": [
       "0.147010883351759"
      ],
      "text/plain": [
       "[1] 0.1470109"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Two ways to estimate FST\n",
    "# FST Estimate 1\n",
    "r=2   # r is the number of populations\n",
    "s2=((FREQpop1-FREQ)^2+(FREQpop2-FREQ)^2)/((r-1))\n",
    "FSTVEC=s2/(FREQ*(1-FREQ))\n",
    "FST1=mean(FSTVEC)\n",
    "\n",
    "# FST Estimate 2: Weir & Cockerham's ESTIMATOR (take into account sample size)\n",
    "n1=max(CEU.data$NCHROBS)/2\n",
    "n2=max(YRI.data$NCHROBS)/2\n",
    "nc=sum(n1+n2)-sum(n1^(2)+n2^(2))/sum(n1+n2)\n",
    "MSA= (1/(r-1))* (n1*(FREQpop1-FREQ)^2+n2*(FREQpop2-FREQ)^2)\n",
    "MSW=(1/(n1-1+n2-1))*(n1* FREQpop1*(1-FREQpop1)+n2*FREQpop2*(1-FREQpop2))\n",
    "numer=sum(MSA-MSW)\n",
    "denom=(sum(MSA+(nc-1)*MSW))\n",
    "FST2=numer/denom\n",
    "\n",
    "FST1\n",
    "FST2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.0953635907786836"
      ],
      "text/latex": [
       "0.0953635907786836"
      ],
      "text/markdown": [
       "0.0953635907786836"
      ],
      "text/plain": [
       "[1] 0.09536359"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "0.102535359565031"
      ],
      "text/latex": [
       "0.102535359565031"
      ],
      "text/markdown": [
       "0.102535359565031"
      ],
      "text/plain": [
       "[1] 0.1025354"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Exercise 3: Obtain Fst values for other HapMap population samples\n",
    "\n",
    "SWITCH2=CEU.data$A1!=ASW.data$A1\n",
    "NewFreq=ASW.data$MAF\n",
    "NewFreq[SWITCH2]=1-ASW.data$MAF[SWITCH2]\n",
    "ASW.data$NewFreq=NewFreq\n",
    "ASW.frq=NewFreq\n",
    "\n",
    "FREQpop1=CEU.frq\n",
    "FREQpop2=ASW.frq\n",
    "FREQ=(FREQpop1+FREQpop2)/2\n",
    "\n",
    "INCLUDE=(FREQ>0 & FREQ<1)\n",
    "FREQpop1=FREQpop1[INCLUDE]\n",
    "FREQpop2=FREQpop2[INCLUDE]\n",
    "FREQ=FREQ[INCLUDE]\n",
    "\n",
    "r=2   # r is the number of populations\n",
    "s2=((FREQpop1-FREQ)^2+(FREQpop2-FREQ)^2)/((r-1))\n",
    "FSTVEC=s2/(FREQ*(1-FREQ))\n",
    "FST1=mean(FSTVEC)\n",
    "\n",
    "# FST Estimate 2: Weir & Cockerham's ESTIMATOR (take into account sample size)\n",
    "n1=max(CEU.data$NCHROBS)/2\n",
    "n2=max(ASW.data$NCHROBS)/2\n",
    "nc=sum(n1+n2)-sum(n1^(2)+n2^(2))/sum(n1+n2)\n",
    "MSA= (1/(r-1))* (n1*(FREQpop1-FREQ)^2+n2*(FREQpop2-FREQ)^2)\n",
    "MSW=(1/(n1-1+n2-1))*(n1* FREQpop1*(1-FREQpop1)+n2*FREQpop2*(1-FREQpop2))\n",
    "numer=sum(MSA-MSW)\n",
    "denom=(sum(MSA+(nc-1)*MSW))\n",
    "FST2=numer/denom\n",
    "\n",
    "FST1\n",
    "FST2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.107286676728355"
      ],
      "text/latex": [
       "0.107286676728355"
      ],
      "text/markdown": [
       "0.107286676728355"
      ],
      "text/plain": [
       "[1] 0.1072867"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "0.107080415538111"
      ],
      "text/latex": [
       "0.107080415538111"
      ],
      "text/markdown": [
       "0.107080415538111"
      ],
      "text/plain": [
       "[1] 0.1070804"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SWITCH3=CEU.data$A1!=CHB.data$A1\n",
    "NewFreq=CHB.data$MAF\n",
    "NewFreq[SWITCH3]=1-CHB.data$MAF[SWITCH3]\n",
    "CHB.data$NewFreq=NewFreq\n",
    "CHB.frq=NewFreq\n",
    "\n",
    "FREQpop1=CEU.frq\n",
    "FREQpop2=CHB.frq\n",
    "FREQ=(FREQpop1+FREQpop2)/2\n",
    "\n",
    "INCLUDE=(FREQ>0 & FREQ<1)\n",
    "FREQpop1=FREQpop1[INCLUDE]\n",
    "FREQpop2=FREQpop2[INCLUDE]\n",
    "FREQ=FREQ[INCLUDE]\n",
    "\n",
    "r=2   # r is the number of populations\n",
    "s2=((FREQpop1-FREQ)^2+(FREQpop2-FREQ)^2)/((r-1))\n",
    "FSTVEC=s2/(FREQ*(1-FREQ))\n",
    "FST1=mean(FSTVEC)\n",
    "\n",
    "# FST Estimate 2: Weir & Cockerham's ESTIMATOR (take into account sample size)\n",
    "n1=max(CEU.data$NCHROBS)/2\n",
    "n2=max(CHB.data$NCHROBS)/2\n",
    "nc=sum(n1+n2)-sum(n1^(2)+n2^(2))/sum(n1+n2)\n",
    "MSA= (1/(r-1))* (n1*(FREQpop1-FREQ)^2+n2*(FREQpop2-FREQ)^2)\n",
    "MSW=(1/(n1-1+n2-1))*(n1* FREQpop1*(1-FREQpop1)+n2*FREQpop2*(1-FREQpop2))\n",
    "numer=sum(MSA-MSW)\n",
    "denom=(sum(MSA+(nc-1)*MSW))\n",
    "FST2=numer/denom\n",
    "\n",
    "FST1\n",
    "FST2\n",
    "\n",
    "# Which of the HapMap population pairs are the most divergent?  \n",
    "# Which populations are the most similar based on Fst? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.2.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}