ACKR1 Gene (FY, DARC, Duffy).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ACKR1 Gene (FY, DARC, Duffy)\n",
    "\n",
    "For gene details, including GRCh37 coordinates, see \n",
    "https://uswest.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000213088;ph=31051;r=1:159203307-159206500\n",
    "\n",
    "For specific information on the 'Duffy-null allele' rs2814778 see Hodgson et. al at https://royalsocietypublishing.org/doi/full/10.1098/rspb.2014.0930"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "selectiontest version:  0.3.22\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pysam\n",
    "from vcf import Reader        # https://pypi.org/project/PyVCF/\n",
    "from Bio import SeqIO\n",
    "from pyliftover import LiftOver\n",
    "from cogent3 import make_table\n",
    "from selectiontest import selectiontest\n",
    "import gzip, pickle\n",
    "import pandas as pd\n",
    "\n",
    "\n",
    "path = \"/Users/helmutsimon/\"\n",
    "if not os.getcwd() == path:\n",
    "    os.chdir(path)\n",
    "    \n",
    "print('selectiontest version: ', selectiontest.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We examine the hypothesis that there is a signal of selection in populations for which the 'Duffy-null allele' rs2814778 segregates. In this cell we compute Tajima's D and $\\rho $ for the ACKR1 gene for all populations. Note that the 1KG data uses GRCh37 coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chrom = 1\n",
    "fname = 'Data sets/1KG variants full/integrated_call_samples_v3.20130502.ALL.panel'\n",
    "panel_all = pd.read_csv(fname, sep=None, engine='python', skipinitialspace=True, index_col=0)\n",
    "vcf_filename = 'Data sets/1KG variants full/ALL.chr' + str(chrom) \\\n",
    "                + '.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz'\n",
    "vcf_file = Reader(filename=vcf_filename, compressed=True, encoding='utf-8')\n",
    "pops = list(set(panel_all['pop']))\n",
    "result_dict = dict()\n",
    "for pop in pops:\n",
    "    panel = panel_all[panel_all['pop'] == pop]\n",
    "    #Use GRCh37 coordinates for ACKR1 gene.\n",
    "    sfs, n, non_seg_snps = selectiontest.vcf2sfs(vcf_file, panel, 1, 159173097, 159176290)\n",
    "    tajd = selectiontest.calculate_D(sfs)\n",
    "    rho = selectiontest.test_neutrality(sfs, reps=200000)\n",
    "    result_dict[pop] = [tajd, rho]\n",
    "results = pd.DataFrame(result_dict, index=['tajd', 'rlnt'])\n",
    "results.to_csv('Google Drive/Genetics/Bayes SFS/Neutrality test/gene_ACKR1_1chrom.csv')\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find populations with selection according to $\\rho $. All are east Asian, African or African ancestry (Barbados)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = 0.5\n",
    "selection_tab = results.loc[:,results.loc['rlnt'] > threshold]\n",
    "select_pops = selection_tab.columns\n",
    "print(len(pops), len(select_pops))\n",
    "select_pops"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find populations in which the 'Duffy-null allele' rs2814778 occurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#First find hg19 coords for the variant\n",
    "lo = LiftOver('hg38', 'hg19')\n",
    "rs2814778_hg19loc = lo.convert_coordinate('chr1', 159204893)\n",
    "print(rs2814778_hg19loc)\n",
    "\n",
    "chrom = 1\n",
    "vcf_filename = 'Data sets/1KG variants full/ALL.chr' + str(chrom) \\\n",
    "                + '.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz'\n",
    "vcf_file = Reader(filename=vcf_filename, compressed=True, encoding='utf-8')\n",
    "snps = vcf_file.fetch(str(chrom), rs2814778_hg19loc[0][1] - 5, rs2814778_hg19loc[0][1] + 5)\n",
    "for snp in snps:\n",
    "    if snp.ID == 'rs2814778':\n",
    "        break\n",
    "print(snp.ID, snp)\n",
    "probands = list()\n",
    "for proband in snp.samples:\n",
    "    gt = proband.gt_alleles\n",
    "    if int(gt[0]) + int(gt[1]) > 0:\n",
    "    #print(proband)\n",
    "        probands.append(proband.sample)\n",
    "fname = 'Data sets/1KG variants full/integrated_call_samples_v3.20130502.ALL.panel'\n",
    "panel_all = pd.read_csv(fname, sep=None, engine='python', skipinitialspace=True, index_col=0)\n",
    "panel_all.head()\n",
    "seg_pops = list(set(panel_all.loc[probands]['pop']))\n",
    "print(len(seg_pops))\n",
    "seg_pops"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is relation between populations undergoing selection and segregating for the 'Duffy-null allele' rs2814778. We find that rs2814778 occurs in all African populations undergoing selection, but in none of the east Asian populations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(set(seg_pops) - set(select_pops))\n",
    "print(set(select_pops) - set(seg_pops))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pool Asian and African populations and attempt to narrow down where selection signal occurs, looking at 800kb segments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analyse_region(chrom, start_hg19, interval, pops):\n",
    "    fname = 'Data sets/1KG variants full/integrated_call_samples_v3.20130502.ALL.panel'\n",
    "    panel_all = pd.read_csv(fname, sep=None, engine='python', skipinitialspace=True, index_col=0)\n",
    "    vcf_filename = 'Data sets/1KG variants full/ALL.chr' + str(chrom) \\\n",
    "                    + '.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz'\n",
    "    vcf_file = Reader(filename=vcf_filename, compressed=True, encoding='utf-8')\n",
    "    panel = panel_all[panel_all['pop'].isin(pops)]\n",
    "    reps = 10000\n",
    "    tajd_results = list()\n",
    "    rho_results = list()\n",
    "    for segment in range(4):\n",
    "        seg_start = start_hg19 + segment * interval\n",
    "        seg_end = seg_start + interval\n",
    "        sfs, n, non_seg_snps = selectiontest.vcf2sfs(vcf_file, panel, 1, seg_start, seg_end)\n",
    "        tajd = selectiontest.calculate_D(sfs)\n",
    "        tajd_results.append(tajd)\n",
    "        rho = selectiontest.test_neutrality(sfs, reps=200000)\n",
    "        rho_results.append(rho)\n",
    "        print(pop, rho, tajd)\n",
    "    return rho_results, tajd_results\n",
    "            \n",
    "chrom = 1\n",
    "start_hg19 = 159173097     # 159176290\n",
    "interval = 800\n",
    "results = pd.DataFrame()\n",
    "pops = ['YRI', 'LWK', 'GWD', 'MSL', 'ESN', 'ASW', 'ASW', 'ACB']\n",
    "#pops = ['LWK', 'GWD', 'ESN', 'ACB']\n",
    "rho_results, tajd_results = analyse_region(chrom, start_hg19, interval, pops)\n",
    "results['afr_rho'] = rho_results\n",
    "results['afr_tajd'] = tajd_results\n",
    "print('\\n')\n",
    "pops = ['JPT', 'BEB', 'CHS', 'KHV', 'CDX', 'CHB']\n",
    "#pops = ['JPT', 'CHS',  'CDX']\n",
    "rho_results, tajd_results = analyse_region(chrom, start_hg19, interval, pops)\n",
    "results['asia_rho'] = rho_results\n",
    "results['asia_tajd'] = tajd_results\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Format result as Latex table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = results.copy()\n",
    "result.columns = pd.MultiIndex.from_arrays([['Africa', 'Africa', 'Asia', 'Asia'], \\\n",
    "                                            ['rho', 'tajd', 'rho', 'tajd']])\n",
    "newix = [str((start_hg19 + interval * i)) for i in range(4)]\n",
    "result.index = newix\n",
    "print(result)\n",
    "result.insert(loc=0, column='Start of 800-bps segment', value=newix)\n",
    "t = make_table(data_frame=result, title=\"caption\", \\\n",
    "               header = ['a\\\\b', 'a\\\\b', 'a\\\\b', 'a\\\\b'], digits=2)\n",
    "t.write(\"Downloads/duffy_sep20.tex\", label=\"tab:duffy\", justify=\"lcccc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate appropriate thresholds for $\\rho$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = selectiontest.compute_threshold(600, 10)\n",
    "t2 = selectiontest.compute_threshold(600, 5)\n",
    "print(t1, t2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Identify SNPs in east Asian populations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = 'Data sets/1KG variants full/integrated_call_samples_v3.20130502.ALL.panel'\n",
    "panel_all = pd.read_csv(fname, sep=None, engine='python', skipinitialspace=True, index_col=0)\n",
    "vcf_filename = 'Data sets/1KG variants full/ALL.chr' + str(chrom) \\\n",
    "                + '.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz'\n",
    "vcf_file = Reader(filename=vcf_filename, compressed=True, encoding='utf-8')\n",
    "pops = ['CDX', 'CHS', 'JPT']\n",
    "panel = panel_all[panel_all['pop'].isin(pops)]\n",
    "snps = vcf_file.fetch('1', 159173897, 159174697)\n",
    "count = 0\n",
    "seg_snps = list()\n",
    "for record in snps:\n",
    "    if record.is_snp:\n",
    "        count += 1\n",
    "    for proband in record.samples:\n",
    "        if proband.sample in panel.index:\n",
    "            gt = proband.gt_alleles\n",
    "            if int(gt[0]) + int(gt[1]) > 0:\n",
    "                seg_snps.append(record.ID)\n",
    "            \n",
    "seg_snps_asia = list(set(seg_snps))\n",
    "print(count)\n",
    "seg_snps_asia"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Identify SNPs in African populations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seg_snps_asia = seg_snps\n",
    "pops = ['GWD', 'LWK', 'ESN', 'ACB']\n",
    "panel = panel_all[panel_all['pop'].isin(pops)]\n",
    "snps = vcf_file.fetch('1', 159173897, 159174697)\n",
    "count = 0\n",
    "seg_snps = list()\n",
    "for record in snps:\n",
    "    if record.is_snp:\n",
    "        count += 1\n",
    "    for proband in record.samples:\n",
    "        if proband.sample in panel.index:\n",
    "            gt = proband.gt_alleles\n",
    "            if int(gt[0]) + int(gt[1]) > 0:\n",
    "                seg_snps.append(record.ID)\n",
    "seg_snps_afr = list(set(seg_snps))\n",
    "print(count)\n",
    "print(seg_snps_afr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(set(seg_snps_asia) - set(seg_snps_afr))\n",
    "print(set(seg_snps_afr) - set(seg_snps_asia))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We check that the ancestral allele T for rs2814778 is a high-confidence call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos37 = 159174683 - 1\n",
    "anc_filename = 'Data sets/human_ancestor_GRCh37_e59/human_ancestor_1.fa'\n",
    "for seq_record in SeqIO.parse(anc_filename, \"fasta\"):\n",
    "    print(seq_record.id)\n",
    "    print(seq_record[pos37 - 5: pos37 +5].seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look for signal of selection in African populations for the gene region chr1:159173097-159176290. We compute likelihood ratio for two models M_0 being a selective model as in Hamlin (2000) and M_1 a neutral model. We use an approximation for the sample from M_0, from SFS samples generated by roc_simulation.py."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CHS 159173097 105 10 -0.9137 0.1220\n"
     ]
    }
   ],
   "source": [
    "chrom = 1\n",
    "fname = 'Data sets/1KG variants full/integrated_call_samples_v3.20130502.ALL.panel'\n",
    "panel_all = pd.read_csv(fname, sep=None, engine='python', skipinitialspace=True, index_col=0)\n",
    "vcf_filename = 'Data sets/1KG variants full/ALL.chr' + str(chrom) \\\n",
    "                    + '.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz'\n",
    "vcf_file = Reader(filename=vcf_filename, compressed=True, encoding='utf-8')\n",
    "   \n",
    "#Use GRCh37 coordinates for FY gene.\n",
    "start = 159173097\n",
    "end = 159176290\n",
    "reps = 100000\n",
    "for pop, popx in zip(['CHS'], ['CHS']):\n",
    "    panel = panel_all[panel_all['pop'] == pop]\n",
    "    n = panel.shape[0]\n",
    "    fname = '/Users/helmutsimon/Google Drive/Genetics/Software/msms/lib/data/sfs_non_neutral_duffy' + \\\n",
    "                popx + '11.pklz'\n",
    "    with gzip.open(fname, 'rb') as q0:  \n",
    "        q0 = pickle.load(q0)\n",
    "    row_sums = q0.sum(axis=1)\n",
    "    variates0 = q0 / row_sums[:, np.newaxis]\n",
    "    \n",
    "    q1 = np.empty((reps, n - 1), dtype=float)\n",
    "    for i, q in enumerate(selectiontest.sample_wf_distribution(n, reps)):\n",
    "                q1[i] = q\n",
    "\n",
    "    sfs, n, non_seg_snps = selectiontest.vcf2sfs(vcf_file, panel, 1, start, end)\n",
    "    odds_ratio = selectiontest.test_neutrality(sfs, variates0=variates0, variates1=q1)\n",
    "    print(pop, start, n, sum(sfs), \"%.4f\" % odds_ratio, \"%.4f\" % 10 ** odds_ratio)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "93"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pop = 'CDX'\n",
    "fname = 'Data sets/1KG variants full/integrated_call_samples_v3.20130502.ALL.panel'\n",
    "panel_all = pd.read_csv(fname, sep=None, engine='python', skipinitialspace=True, index_col=0)\n",
    "panel = panel_all[panel_all['pop'] == pop]\n",
    "panel.shape[0]"
   ]
  },
  {
   "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.7.3"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}