ACKR1_gene_chimpanzee.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ACKR1_gene_chimpanzee\n",
    "\n",
    "Data obtained from http://biologiaevolutiva.org/greatape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pathlib\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from collections import Counter\n",
    "from selectiontest import selectiontest\n",
    "import pysam\n",
    "from vcf import Reader        # https://pypi.org/project/PyVCF/\n",
    "from Bio import SeqIO\n",
    "from selectiontest import selectiontest\n",
    "import gzip, pickle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set hg18 (NCBI36 Ensembl release 54) coordinates for gene ACKR1 from http://may2009.archive.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000213088"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "chrom = 'chr1'\n",
    "ACKR1_start_hg18 = 157439721\n",
    "ACKR1_end_hg18 =   157442914 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate SFS, ignoring variants with missing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total SNPs processed        =  61\n",
      "SNPs without missing data   =  46\n",
      "SFS  [0 3 0 4 0 4 2 0 1 0 0 0 0 0 1 1 0 0 0 3 0 0 0 0]\n",
      "Number of segregating sites =  19\n",
      "Sample size                 =  25\n",
      "SNPs not segregating in sample =  27\n"
     ]
    }
   ],
   "source": [
    "def vcf2sfs2(vcf_file, panel, coord, start, end, select_chr=True):\n",
    "    n = panel.shape[0] \n",
    "    if not select_chr:\n",
    "        n = 2 * n\n",
    "    snps = vcf_file.fetch(str(coord), start, end)\n",
    "    count, valid_snp_count = 0, 0\n",
    "    allele_counts = list()\n",
    "    non_seg_snps = list()\n",
    "    for record in snps:\n",
    "        allele_count = 0\n",
    "        if record.is_snp:\n",
    "            count += 1\n",
    "            missing_data = False\n",
    "            for sample in record.samples:\n",
    "                if sample.sample in panel.index:\n",
    "                    if sample.gt_bases is None:\n",
    "                        missing_data = True\n",
    "                    else:\n",
    "                        gt = sample.gt_alleles\n",
    "                        if select_chr:\n",
    "                            allele_count += int(gt[0])\n",
    "                        else:\n",
    "                            allele_count += int(gt[0]) + int(gt[1])\n",
    "            if missing_data:\n",
    "                continue\n",
    "            valid_snp_count += 1    \n",
    "            if 0 < allele_count < n:    #Some SNPs may not segregate in some subpopulations.\n",
    "                allele_counts.append(allele_count)\n",
    "            else:\n",
    "                non_seg_snps.append(record.POS)\n",
    "    sfs_c = Counter(allele_counts)\n",
    "    del sfs_c[0]\n",
    "    sfs = np.zeros(n - 1, int)\n",
    "    for i in sfs_c.keys():\n",
    "        sfs[i - 1] = sfs_c[i]\n",
    "    print('Total SNPs processed        = ', count)\n",
    "    print('SNPs without missing data   = ', valid_snp_count)\n",
    "    return sfs, n, non_seg_snps, count, valid_snp_count\n",
    "\n",
    "\n",
    "vcf_filename = '/Users/helmutsimon/OneDrive - Australian National University/Data sets/Pan_troglodytes/Pan_troglodytes.vcf.gz'\n",
    "fname = pathlib.Path(vcf_filename)\n",
    "assert fname.exists(), f'No such file: {fname}'  # check that the file exists\n",
    "vcf_file = Reader(filename=vcf_filename, compressed=True, encoding='utf-8')\n",
    "panel_name = '/Users/helmutsimon/OneDrive - Australian National University/Data sets/Pan_troglodytes/Pan_trog_panel.csv'\n",
    "panel = pd.read_csv(panel_name, sep=',', index_col=0)\n",
    "sfs, n, non_seg_snps, count, valid_snp_count = vcf2sfs2(vcf_file, panel, chrom, ACKR1_start_hg18, ACKR1_end_hg18)\n",
    "\n",
    "print('SFS ', sfs)\n",
    "print('Number of segregating sites = ', sum(sfs))\n",
    "print('Sample size                 = ', n)\n",
    "print('SNPs not segregating in sample = ', len(non_seg_snps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate rho and Tajima's D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rho =  -1.7842814054136191\n",
      "Tajimas D =  0.9777780110383649\n"
     ]
    }
   ],
   "source": [
    "print('Rho = ', selectiontest.test_neutrality(sfs, reps=1000000))\n",
    "print('Tajimas D = ', selectiontest.calculate_D(sfs))"
   ]
  },
  {
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}