1000GenomesProject/AnalysisCmedi.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from browser import Pops, SPops\n",
    "\n",
    "import numpy as np\n",
    "from tools1kg import path_1000G\n",
    "from tqdm.notebook import trange, tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "import itertools\n",
    "import seaborn as sns\n",
    "from scipy.ndimage import gaussian_filter1d\n",
    "plt.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "Pops = {i:dict(zip(['pop', 'description'], Pops[i].split('--'))) for i in range(26)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('/home/enes/pyslurm/') # https://github.com/enesdilber/pyslurm/\n",
    "sys.path.append('/home/enes/bim/') # https://github.com/enesdilber/bim/\n",
    "from pyslurm import Slurm\n",
    "from utils import intersect_with_weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "Train = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "BIM = '/home/enes/bim/BIM.py'\n",
    "locvars = ['Chromosome', 'start', 'end']\n",
    "#sys.path.append(bim_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "slurm = Slurm()\n",
    "srun = slurm.batch('module load python3.8-anaconda/2020.07', \n",
    "                   '#mem-per-cpu=6000',\n",
    "                   '#job-name=\"1kg\"',\n",
    "                   '#time=0-10:00:00',\n",
    "                   '#cpus-per-task=1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get statistics\n",
    "Below calculates the several genome scan statistic including btree and bsfs. You should give the `Settings` file. It's a json file and will be red as a python dictionary. Its md5 hash will make sure the results are reproducible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "SFS = np.load(\"popSFS.npy\", allow_pickle=True).item()\n",
    "N = {i: np.array(SFS[i]).shape[1]+1 for i in range(26)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "if Train:\n",
    "    gsettings = ['--eta=etas.json', '--r2t=0.0025', '--r2s=0.0025', '--wsz=100000', '--ssz=25000', '--log_pdf=logfr']\n",
    "    jobs = {}\n",
    "    for chrno in range(1, 23):\n",
    "        for popid in range(26):     \n",
    "            \n",
    "            Ni = str(N[popid])\n",
    "            popid = str(popid)\n",
    "            chrno = str(chrno)\n",
    "            job = ' '.join(['python', BIM, path_1000G(chrno), Ni, '--pop='+popid, '--out=raw/none'+popid+'_'+chrno+'Cmedi.csv']+gsettings)\n",
    "            jobs[chrno, popid] = srun.run(job)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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>JOBID</th>\n",
       "      <th>PARTITION</th>\n",
       "      <th>NAME</th>\n",
       "      <th>USER</th>\n",
       "      <th>ACCOUNT</th>\n",
       "      <th>ST</th>\n",
       "      <th>TIME</th>\n",
       "      <th>NODES</th>\n",
       "      <th>NODELIST(REASON)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>19564801</td>\n",
       "      <td>standard</td>\n",
       "      <td>sys/dash</td>\n",
       "      <td>enes</td>\n",
       "      <td>stats_de</td>\n",
       "      <td>R</td>\n",
       "      <td>1-09:57:02</td>\n",
       "      <td>1</td>\n",
       "      <td>gl3122</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      JOBID PARTITION      NAME  USER   ACCOUNT ST        TIME  NODES  \\\n",
       "0  19564801  standard  sys/dash  enes  stats_de  R  1-09:57:02      1   \n",
       "\n",
       "  NODELIST(REASON)  \n",
       "0           gl3122  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "slurm.queue()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "if Train:\n",
    "    gsettings = ['--eta=etas.json', '--r2t=0.0025', '--r2s=0.0025', '--wsz=10000', '--ssz=5000', '--log_pdf=logfr']\n",
    "    jobs = {}\n",
    "    for chrno in range(1, 23):\n",
    "        for popid in range(26):     \n",
    "            \n",
    "            Ni = str(N[popid])\n",
    "            popid = str(popid)\n",
    "            chrno = str(chrno)\n",
    "            job = ' '.join(['python', BIM, path_1000G(chrno), Ni, '--pop='+popid, '--out=raw/'+popid+'_'+chrno+'Cmedi.csv']+gsettings)\n",
    "            jobs[chrno, popid] = srun.run(job)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merge dataframes\n",
    "Above code will create 26x22 csv files for each population chromosome combination, I will merge them for each population and calculate Neutrality tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "56e4a2371f1d4410af48be06d16865dd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=26.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "if Train:\n",
    "    def read_and_load(pop, chrno):\n",
    "        df2 = pd.read_csv('raw/'+str(pop)+'_'+str(chrno)+'Cmedi.csv', comment = '#')\n",
    "        df2.insert(0, \"Chromosome\", chrno)\n",
    "        return df2\n",
    "    for pop in trange(26):\n",
    "        pd.concat([read_and_load(pop, chrno) for chrno in range(1, 23)]).\\\n",
    "        to_csv('merged/s'+str(pop)+'Cmedi.csv', index = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.read_csv()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Take avarage statistic along the chromosome for each window\n",
    "I will calculate a base value for each windows by calculating mean and median for each statistic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f7bb5025d02e4c02b3b79af3dda62350",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=9.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "59dfc9651d8d47b1b1e6d57b56962725",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=26.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "35b5180026f3421bb61c9be4fafc5067",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=26.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1bbee3b509b9461e9fa972da7f3ca937",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=26.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c3634ef26603414db2596f6ede548186",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=26.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c2a88455e38d41cabbc4d2e1433896d0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=26.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f74c163490434804b5f7adfe9360e82c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=26.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "66041d0a526047ff98615e6508f135fb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=26.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e4cba1dc413b4f6a8fc8bfa3655d2d8e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=26.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1fe74073062c4fc1b4ed2bdedeff934a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=26.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "if Train:\n",
    "    statnames = ['SS', 'FayH', 'FulD', 'ZngE', 'TajD', 'bsfs', 'FerL', 'Colless', 'btree']\n",
    "    statmedis = {}\n",
    "    statmeans = {}\n",
    "    for stat in tqdm(statnames):\n",
    "        Stat = np.zeros((576194, 26)) # number of rows and population\n",
    "        for i in trange(26):\n",
    "            the_path = 'merged/s'+str(i)+'Cmedi.csv'\n",
    "            df = pd.read_csv(the_path)\n",
    "            Stat[:,i] = df[stat]\n",
    "        statmedis[stat] = np.median(Stat, 1)\n",
    "        statmeans[stat] = np.mean(Stat, 1)\n",
    "\n",
    "    dxmedi = df[['Chromosome', 'start', 'end']].copy()\n",
    "    dxmean = df[['Chromosome', 'start', 'end']].copy()\n",
    "    for stat in statnames:\n",
    "        dxmedi['medi_'+stat] = statmedis[stat]\n",
    "        dxmean['mean_'+stat] = statmeans[stat]\n",
    "    dxmedi = dxmedi[dxmedi['medi_SS']>5]\n",
    "    dxmean = dxmean[dxmean['mean_SS']>5]\n",
    "    dxmean.to_csv('meanstats.csv', index = False)\n",
    "    dxmedi.to_csv('medianstats.csv', index = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Change point detection\n",
    "I will calculate the change points for each statitic and for each population. I will assume a fixed number of windows for each statistic and chromosome which is $n_{cp}=$ `chromosome length/1e7`. This number will give >300 regions along the chromosome, which is roughly the number of beneficial mutations along the genome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "if Train:\n",
    "    avglen = 1e7\n",
    "    jobs = {}\n",
    "    for pop in range(26):\n",
    "        for chrno in range(1, 23):\n",
    "            job = ' '.join(['python segmenter.py', str(pop), str(chrno),str(avglen), 'raw/', 'medi'])\n",
    "            jobs[pop, chrno] = srun.run(job)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After the jobs are done, I merge them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b7d35850c6b84a46af8d9a0e21fcfda2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=26.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "if Train:\n",
    "    \n",
    "    for pop in trange(26):\n",
    "        dfs = []\n",
    "        for chrno in range(1, 23):\n",
    "            path = 'raw/s'+str(pop)+'_'+str(chrno)+'medi.csv'\n",
    "            dfs.append(pd.read_csv(path))\n",
    "\n",
    "        pd.concat(dfs).to_csv('merged/s'+str(pop)+'Cmedi.csv', index = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate Lag for each population\n",
    "Because of LD all statistics have autocorrelation. In order to calculate avaraged statistic for each region we need the variance of the avaraged statistic. To calculate this variance, we need to take account the autocovariance matrix. Here we calculate it for all statistics and for all populations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "if Train:\n",
    "\n",
    "    stats = ['SSCmedi', 'TajDCmedi', 'CollessCmedi', 'bsfsCmedi', 'btreeCmedi', \n",
    "             'FulDCmedi', 'ZngECmedi', 'FerLCmedi', 'FayHCmedi']\n",
    "\n",
    "    jobs = []\n",
    "    for stat in stats:\n",
    "        job = 'python calc_lag.py ' + stat + ' 2000 merged Cmedi'\n",
    "        jobs.append(srun.run(job))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Empirical p-values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bd6e08afc4a24d8390e57f2f50c55e80",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=26.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "if Train:\n",
    "    \n",
    "    stats = ['SSCmedi', 'TajDCmedi', 'CollessCmedi', 'bsfsCmedi', 'btreeCmedi', \n",
    "             'FulDCmedi', 'ZngECmedi', 'FerLCmedi', 'FayHCmedi']\n",
    "\n",
    "    ann = pd.read_csv('gene_result.txt', sep=\"\\t\")\n",
    "    ann = ann.rename(columns = {'start_position_on_the_genomic_accession': 'start',\n",
    "                                'end_position_on_the_genomic_accession': 'end'})[['chromosome', 'start', 'end', 'Symbol']]\n",
    "    ann = ann[~((ann['chromosome'] == 'X') | ((ann['chromosome'] == 'Y')))]\n",
    "\n",
    "    ann['chromosome'] = ann['chromosome'].astype('int')\n",
    "    ann = ann.sort_values(['chromosome', 'start', 'end']).reset_index(drop = True)\n",
    "\n",
    "    for pop in trange(26):\n",
    "        pop = str(pop)\n",
    "\n",
    "        df = pd.read_csv('merged/s'+pop+'Cmedi.csv')\n",
    "        df['start'] = df['Chromosome']+(df['start']/1e9)\n",
    "        df['end'] = df['Chromosome']+(df['end']/1e9)\n",
    "\n",
    "        bs = df['start'].to_numpy()\n",
    "        be = df['end'].to_numpy()\n",
    "\n",
    "        cps = (ann['chromosome']+(ann['start']/1e9)).to_numpy()\n",
    "        cpe = (ann['chromosome']+(ann['end']/1e9)).to_numpy()\n",
    "\n",
    "        dx = pd.DataFrame(intersect_with_weights(bs, be, cps, cpe, df[stats].to_numpy()))\n",
    "        for i, stat in enumerate(stats):\n",
    "            dx['ep_'+stat] = [x[i] for x in dx['val']]\n",
    "            dx['ep_'+stat] = dx['ep_'+stat].rank(pct=True)\n",
    "\n",
    "        dx.insert(0,\"Chromosome\", dx['start'].astype('int'))\n",
    "        dx['start'] = np.round(1e9*(dx['start']-dx['Chromosome'])).astype('int')\n",
    "        dx['end'] = np.round(1e9*(dx['end']-dx['Chromosome'])).astype('int')\n",
    "\n",
    "        df = ann.rename(columns={'chromosome':'Chromosome'})\n",
    "        df = df.merge(dx, on = locvars)[locvars+['Symbol']+['ep_'+stat for stat in stats]]\n",
    "        df.to_csv('gene_scores/pop/'+pop+'Cmedi.csv', index = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stats by the Region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [],
   "source": [
    "genes = [{'Gene':'EDAR', 'pop':0},\n",
    "         {'Gene':'CD5', 'pop':0},\n",
    "         {'Gene':'LCT', 'pop':5},\n",
    "         {'Gene':'SLC45A2', 'pop':5},\n",
    "         {'Gene':'HERC2', 'pop':5},\n",
    "         {'Gene':'SLC24A5', 'pop':5},\n",
    "         {'Gene':'CD36', 'pop':10},\n",
    "         {'Gene':'APOL1', 'pop':10}]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "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>Gene</th>\n",
       "      <th>pop</th>\n",
       "      <th>Chromosome</th>\n",
       "      <th>start</th>\n",
       "      <th>end</th>\n",
       "      <th>btreeCmedi</th>\n",
       "      <th>bsfsCmedi</th>\n",
       "      <th>TajDCmedi</th>\n",
       "      <th>CollessCmedi</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>LCT</td>\n",
       "      <td>CEU</td>\n",
       "      <td>2</td>\n",
       "      <td>135787850</td>\n",
       "      <td>135837195</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.454545</td>\n",
       "      <td>0.636364</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>SLC45A2</td>\n",
       "      <td>CEU</td>\n",
       "      <td>5</td>\n",
       "      <td>33944623</td>\n",
       "      <td>33984693</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.666667</td>\n",
       "      <td>0.444444</td>\n",
       "      <td>0.888889</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>HERC2</td>\n",
       "      <td>CEU</td>\n",
       "      <td>15</td>\n",
       "      <td>28111040</td>\n",
       "      <td>28322179</td>\n",
       "      <td>0.069767</td>\n",
       "      <td>0.023256</td>\n",
       "      <td>0.069767</td>\n",
       "      <td>0.116279</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>SLC24A5</td>\n",
       "      <td>CEU</td>\n",
       "      <td>15</td>\n",
       "      <td>48120990</td>\n",
       "      <td>48142672</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>EDAR</td>\n",
       "      <td>CHB</td>\n",
       "      <td>2</td>\n",
       "      <td>108894471</td>\n",
       "      <td>108989256</td>\n",
       "      <td>0.583333</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.166667</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>CD5</td>\n",
       "      <td>CHB</td>\n",
       "      <td>11</td>\n",
       "      <td>61093963</td>\n",
       "      <td>61127852</td>\n",
       "      <td>0.125000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>CD36</td>\n",
       "      <td>YRI</td>\n",
       "      <td>7</td>\n",
       "      <td>80602207</td>\n",
       "      <td>80679277</td>\n",
       "      <td>0.062500</td>\n",
       "      <td>0.062500</td>\n",
       "      <td>0.187500</td>\n",
       "      <td>0.125000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>APOL1</td>\n",
       "      <td>YRI</td>\n",
       "      <td>22</td>\n",
       "      <td>36253071</td>\n",
       "      <td>36267531</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Gene  pop  Chromosome      start        end  btreeCmedi  bsfsCmedi  \\\n",
       "2      LCT  CEU           2  135787850  135837195    1.000000   0.454545   \n",
       "3  SLC45A2  CEU           5   33944623   33984693    1.000000   0.666667   \n",
       "4    HERC2  CEU          15   28111040   28322179    0.069767   0.023256   \n",
       "5  SLC24A5  CEU          15   48120990   48142672    0.000000   0.000000   \n",
       "0     EDAR  CHB           2  108894471  108989256    0.583333   0.000000   \n",
       "1      CD5  CHB          11   61093963   61127852    0.125000   0.000000   \n",
       "6     CD36  YRI           7   80602207   80679277    0.062500   0.062500   \n",
       "7    APOL1  YRI          22   36253071   36267531    0.000000   0.000000   \n",
       "\n",
       "   TajDCmedi  CollessCmedi  \n",
       "2   0.636364      1.000000  \n",
       "3   0.444444      0.888889  \n",
       "4   0.069767      0.116279  \n",
       "5   0.000000      0.000000  \n",
       "0   0.166667      1.000000  \n",
       "1   0.000000      0.000000  \n",
       "6   0.187500      0.125000  \n",
       "7   0.000000      0.000000  "
      ]
     },
     "execution_count": 127,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats = ['btreeCmedi', 'bsfsCmedi', 'TajDCmedi', 'CollessCmedi']\n",
    "signif = 0.05\n",
    "pop = ''\n",
    "res = []\n",
    "for gene in genes:\n",
    "    \n",
    "    pop = gene['pop']\n",
    "    cur = gene.copy()\n",
    "    cur['pop'] = Pops[pop]['pop']\n",
    "    \n",
    "    df = pd.read_csv('gene_scores/pop/'+str(gene['pop'])+'Cmedi.csv')    \n",
    "    df = df[df['Symbol'] == gene['Gene']]\n",
    "    df = df[locvars+['ep_'+stat for stat in stats]]\n",
    "    \n",
    "    \n",
    "    cur['Chromosome'] =  df['Chromosome'].iloc[0]\n",
    "    cur['start'] = df['start'].iloc[0]\n",
    "    cur['end'] = df['end'].iloc[-1]\n",
    "    \n",
    "    for stat in stats:\n",
    "        if stat[:7] == 'Colless':\n",
    "            cur[stat] = (df['ep_'+stat]>(1-signif)).mean()        \n",
    "        else:\n",
    "            cur[stat] = (df['ep_'+stat]<signif).mean()        \n",
    "    \n",
    "    res.append(cur)\n",
    "\n",
    "res = pd.DataFrame(res)\n",
    "res.sort_values(['pop', 'Chromosome', 'start'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Replication of 1000 genomes selection paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [],
   "source": [
    "st2 = [{'Gene':'EDAR', 'Chromosome':2, 'start':109_510_927, 'end':109_605_828, 'pop':0},\n",
    "       {'Gene':'CD5', 'Chromosome':11,'start':60_869_867, 'end':60_896_324, 'pop':0},\n",
    "       {'Gene':'LCT', 'Chromosome':2,'start':136_545_410, 'end':136_594_750, 'pop':5},\n",
    "       {'Gene':'SLC45A2', 'Chromosome':5,'start':33_944_721, 'end':33_984_835, 'pop':5},\n",
    "       {'Gene':'HERC2', 'Chromosome':15,'start':28_356_186, 'end':28_567_298, 'pop':5},\n",
    "       {'Gene':'SLC24A5', 'Chromosome':15,'start':48_413_169, 'end':48_434_869, 'pop':5},\n",
    "       {'Gene':'CD36', 'Chromosome':7,'start':79_998_891, 'end':80_308_593, 'pop':10},\n",
    "       {'Gene':'APOL1', 'Chromosome':22,'start':36_649_117, 'end':36_663_576, 'pop':10}]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "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>Gene</th>\n",
       "      <th>Chromosome</th>\n",
       "      <th>start</th>\n",
       "      <th>end</th>\n",
       "      <th>pop</th>\n",
       "      <th>btreeCmedi</th>\n",
       "      <th>bsfsCmedi</th>\n",
       "      <th>TajDCmedi</th>\n",
       "      <th>CollessCmedi</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>LCT</td>\n",
       "      <td>2</td>\n",
       "      <td>136545410</td>\n",
       "      <td>136594750</td>\n",
       "      <td>CEU</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.545455</td>\n",
       "      <td>0.454545</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>SLC45A2</td>\n",
       "      <td>5</td>\n",
       "      <td>33944721</td>\n",
       "      <td>33984835</td>\n",
       "      <td>CEU</td>\n",
       "      <td>0.888889</td>\n",
       "      <td>0.666667</td>\n",
       "      <td>0.444444</td>\n",
       "      <td>0.888889</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>HERC2</td>\n",
       "      <td>15</td>\n",
       "      <td>28356186</td>\n",
       "      <td>28567298</td>\n",
       "      <td>CEU</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.279070</td>\n",
       "      <td>0.720930</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>SLC24A5</td>\n",
       "      <td>15</td>\n",
       "      <td>48413169</td>\n",
       "      <td>48434869</td>\n",
       "      <td>CEU</td>\n",
       "      <td>0.800000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.600000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>EDAR</td>\n",
       "      <td>2</td>\n",
       "      <td>109510927</td>\n",
       "      <td>109605828</td>\n",
       "      <td>CHB</td>\n",
       "      <td>0.473684</td>\n",
       "      <td>0.473684</td>\n",
       "      <td>0.894737</td>\n",
       "      <td>0.263158</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>CD5</td>\n",
       "      <td>11</td>\n",
       "      <td>60869867</td>\n",
       "      <td>60896324</td>\n",
       "      <td>CHB</td>\n",
       "      <td>0.571429</td>\n",
       "      <td>0.571429</td>\n",
       "      <td>0.714286</td>\n",
       "      <td>0.571429</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>CD36</td>\n",
       "      <td>7</td>\n",
       "      <td>79998891</td>\n",
       "      <td>80308593</td>\n",
       "      <td>YRI</td>\n",
       "      <td>0.158730</td>\n",
       "      <td>0.063492</td>\n",
       "      <td>0.111111</td>\n",
       "      <td>0.206349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>APOL1</td>\n",
       "      <td>22</td>\n",
       "      <td>36649117</td>\n",
       "      <td>36663576</td>\n",
       "      <td>YRI</td>\n",
       "      <td>0.250000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.250000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Gene  Chromosome      start        end  pop  btreeCmedi  bsfsCmedi  \\\n",
       "2      LCT           2  136545410  136594750  CEU    1.000000   0.545455   \n",
       "3  SLC45A2           5   33944721   33984835  CEU    0.888889   0.666667   \n",
       "4    HERC2          15   28356186   28567298  CEU    1.000000   0.279070   \n",
       "5  SLC24A5          15   48413169   48434869  CEU    0.800000   0.000000   \n",
       "0     EDAR           2  109510927  109605828  CHB    0.473684   0.473684   \n",
       "1      CD5          11   60869867   60896324  CHB    0.571429   0.571429   \n",
       "6     CD36           7   79998891   80308593  YRI    0.158730   0.063492   \n",
       "7    APOL1          22   36649117   36663576  YRI    0.250000   0.000000   \n",
       "\n",
       "   TajDCmedi  CollessCmedi  \n",
       "2   0.454545      1.000000  \n",
       "3   0.444444      0.888889  \n",
       "4   0.720930      1.000000  \n",
       "5   0.600000      1.000000  \n",
       "0   0.894737      0.263158  \n",
       "1   0.714286      0.571429  \n",
       "6   0.111111      0.206349  \n",
       "7   0.000000      0.250000  "
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats = ['btreeCmedi', 'bsfsCmedi', 'TajDCmedi', 'CollessCmedi']\n",
    "signif = 0.05\n",
    "pop = ''\n",
    "res = []\n",
    "for st in st2:\n",
    "    \n",
    "    if st['pop'] != pop:\n",
    "        df = pd.read_csv('merged/s'+str(st['pop'])+'Cmedi.csv').fillna(0)\n",
    "        df = df[locvars+stats]\n",
    "        for stat in stats:\n",
    "            if  stat[:7] == 'Colless':\n",
    "                df[stat] = -df[stat]\n",
    "            df['ep_'+stat] = df[stat].rank(pct=True)<signif        \n",
    "\n",
    "    pop = st['pop']\n",
    "    chrno = st['Chromosome']\n",
    "    rs, re = st['start'], st['end']\n",
    "    \n",
    "    nrs, nre = np.round(rs, -4), np.round(re, -4)\n",
    "    re = np.where(nre>re, nre, nre+1e4).astype('int')\n",
    "    rs = np.where(nrs<rs, nrs, nrs-1e4).astype('int')\n",
    "    \n",
    "    dx = df[df['Chromosome'] == chrno].copy().reset_index(drop = True).copy()\n",
    "    dx = dx[(dx['start']>=rs) & (dx['end']<=re)]\n",
    "    \n",
    "    for stat in stats:\n",
    "        st['pop'] = Pops[int(pop)]['pop']\n",
    "        st[stat] = dx['ep_'+stat].mean()\n",
    "    res.append(st)\n",
    "res = pd.DataFrame(res)\n",
    "res.sort_values(['pop', 'Chromosome', 'start'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Genome-scan p-values\n",
    "Now I plot the genome-scans of all statistics for each population and save the results along with the genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "if Train:\n",
    "\n",
    "    stats = ['SSCmedi', 'TajDCmedi', 'CollessCmedi', 'bsfsCmedi', 'btreeCmedi', \n",
    "             'FulDCmedi', 'ZngECmedi', 'FerLCmedi', 'FayHCmedi']\n",
    "\n",
    "    jobs = []\n",
    "    for stat in stats:\n",
    "        for pop in range(26):\n",
    "            job = 'python pvalscan.py ' +str(pop)+' '+ stat +' Cmedi'\n",
    "            jobs.append(srun.run(job))"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}