SNP detection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SNP Detection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting the dependencies"  
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/senne/nanopore/SNP/known_SNP_sequence/SNP_sequence.fasta\n",
      "Number of SNPs in reference file:\n",
      "52\n"
     ]
    }
   ],
   "source": [
    "snpFile   = '/home/senne/nanopore/SNP/known_SNP_sequence/SNP_sequence.fasta'  \n",
    "readFile  = '/home/senne/nanopore/SNP/Nanopore_data/2_potential_snp_amplicons_3mism.fasta'\n",
    "resultDir = '/home/senne/nanopore/SNP/results_yannick'\n",
    "fastqFile = '/home/senne/nanopore/SNP/Nanopore_data/ligatedSNPs.fastq'\n",
    "bwa       = '/opt/tools/bwa-0.7.15'     # v0.7.5\n",
    "samtools  = '/opt/tools/samtools-1.3.1' # v1.3.1\n",
    "bcftools  = '/opt/tools/bcftools-1.3.1' # v1.3.1\n",
    "\n",
    "# Check\n",
    "!ls {snpFile}\n",
    "print('Number of SNPs in reference file:')\n",
    "!grep -c \">\" {snpFile}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Info on the fastq file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(14324, 507.74413571628037, 411.0, 5850, 324.59961624575982)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "lengthData = []\n",
    "\n",
    "with open(fastqFile, 'r') as inFile:\n",
    "  lineCount = 0\n",
    "  \n",
    "  for line in inFile:\n",
    "    lineCount += 1\n",
    "  \n",
    "    if lineCount % 4 == 2:\n",
    "      lengthData.append(len(line.rstrip()))\n",
    "\n",
    "# Display info\n",
    "len (lengthData), np.mean(lengthData), np.median(lengthData) , max(lengthData), np.std(lengthData)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Map reads to reference sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[bwa_index] Pack FASTA... 0.00 sec\n",
      "[bwa_index] Construct BWT for the packed sequence...\n",
      "[bwa_index] 0.00 seconds elapse.\n",
      "[bwa_index] Update BWT... 0.00 sec\n",
      "[bwa_index] Pack forward-only FASTA... 0.00 sec\n",
      "[bwa_index] Construct SA from BWT and Occ... 0.00 sec\n",
      "[main] Version: 0.7.15-r1140\n",
      "[main] CMD: /opt/tools/bwa-0.7.15 index /home/senne/nanopore/SNP/known_SNP_sequence/SNP_sequence.fasta\n",
      "[main] Real time: 0.016 sec; CPU: 0.013 sec\n",
      "[M::bwa_idx_load_from_disk] read 0 ALT contigs\n",
      "[M::process] read 23057 sequences (1988758 bp)...\n",
      "[M::mem_process_seqs] Processed 23057 reads in 1.191 CPU sec, 1.189 real sec\n",
      "[main] Version: 0.7.15-r1140\n",
      "[main] CMD: /opt/tools/bwa-0.7.15 mem -a /home/senne/nanopore/SNP/known_SNP_sequence/SNP_sequence.fasta /home/senne/nanopore/SNP/Nanopore_data/2_potential_snp_amplicons_3mism.fasta\n",
      "[main] Real time: 1.255 sec; CPU: 1.252 sec\n",
      "Done\n"
     ]
    }
   ],
   "source": [
    "# Map reads to reference sequences\n",
    "\n",
    "# Build index of the references\n",
    "!{bwa} index {snpFile}\n",
    "\n",
    "# Map reads\n",
    "!{bwa} mem -a {snpFile} {readFile} > {resultDir}/SNPforID.sam\n",
    "\n",
    "# Make sorted bam and index\n",
    "!{samtools} view -Sbu {resultDir}/SNPforID.sam | {samtools} sort -o {resultDir}/SNPforID_sorted.bam -\n",
    "!{samtools} index {resultDir}/SNPforID_sorted.bam {resultDir}/SNPforID_sorted.bam.bai\n",
    "print('Done')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "23061 + 0 in total (QC-passed reads + QC-failed reads)\r\n",
      "0 + 0 secondary\r\n",
      "4 + 0 supplementary\r\n",
      "0 + 0 duplicates\r\n",
      "7824 + 0 mapped (33.93% : N/A)\r\n",
      "0 + 0 paired in sequencing\r\n",
      "0 + 0 read1\r\n",
      "0 + 0 read2\r\n",
      "0 + 0 properly paired (N/A : N/A)\r\n",
      "0 + 0 with itself and mate mapped\r\n",
      "0 + 0 singletons (N/A : N/A)\r\n",
      "0 + 0 with mate mapped to a different chr\r\n",
      "0 + 0 with mate mapped to a different chr (mapQ>=5)\r\n"
     ]
    }
   ],
   "source": [
    "# Display some stats\n",
    "\n",
    "!{samtools} flagstat {resultDir}/SNPforID_sorted.bam\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Generate vcf file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: Neither --ploidy nor --ploidy-file given, assuming all sites are diploid\n",
      "[mpileup] 1 samples in 1 input files\n",
      "Done\n"
     ]
    }
   ],
   "source": [
    "# Generate vcf file from bam file. Needs the reference and its index file \n",
    "\n",
    "# Note: the commands below are for samtools and bcftools v1.3.1 and will not work on v0.1.19!\n",
    "\n",
    "# Reporting all positions\n",
    "!{samtools} mpileup -d 100000 -uf {snpFile} {resultDir}/SNPforID_sorted.bam | {bcftools} call -V indels -m - > {resultDir}/SNPforID_sorted.bam.vcf\n",
    "\n",
    "print('Done')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Generate the SNP profile "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Got data for 52 SNPs:\n",
      "  rs1005533 (97)  A (62.9 %)  G (37.1 %)\n",
      "  rs1015250 (63)  C (3.2 %)  G (96.8 %)\n",
      "  rs1024116 (98)  A (40.8 %)  G (59.2 %)\n",
      "  rs1028528 (134)  A (3.0 %)  G (97.0 %)\n",
      "  rs1029047 (29)  A (100.0 %)\n",
      "  rs1031825 (76)  A (26.3 %)  C (73.7 %)\n",
      "  rs10495407 (177)  A (60.5 %)  G (39.5 %)\n",
      "  rs1335873 (119)  A (1.7 %)  T (98.3 %)\n",
      "  rs1355366 (158)  A (5.1 %)  G (94.9 %)\n",
      "  rs1357617 (218)  A (100.0 %)\n",
      "  rs1360288 (183)  C (60.1 %)  T (39.9 %)\n",
      "  rs1382387 (200)  G (2.5 %)  T (97.5 %)\n",
      "  rs1413212 (164)  A (65.2 %)  G (34.8 %)\n",
      "  rs1454361 (195)  A (100.0 %)\n",
      "  rs1463729 (148)  A (97.3 %)\n",
      "  rs1490413 (224)  A (98.2 %)\n",
      "  rs1493232 (69)  A (73.9 %)  C (26.1 %)\n",
      "  rs1528460 (83)  C (2.4 %)  T (97.6 %)\n",
      "  rs1886510 (113)  C (6.2 %)  T (93.8 %)\n",
      "  rs1979255 (301)  C (98.0 %)\n",
      "  rs2016276 (69)  A (95.7 %)\n",
      "  rs2040411 (76)  A (0.0 %)  G (100.0 %)\n",
      "  rs2046361 (110)  A (3.6 %)  T (96.4 %)\n",
      "  rs2056277 (199)  C (47.7 %)  T (52.3 %)\n",
      "  rs2076848 (123)  A (73.2 %)  T (26.8 %)\n",
      "  rs2107612 (48)  A (0.0 %)  G (100.0 %)\n",
      "  rs2111980 (203)  A (62.1 %)  G (37.9 %)\n",
      "  rs251934 (123)  C (1.6 %)  T (98.4 %)\n",
      "  rs2830795 (72)  A (76.4 %)  G (23.6 %)\n",
      "  rs2831700 (290)  A (99.3 %)\n",
      "  rs354439 (127)  A (73.2 %)  T (26.8 %)\n",
      "  rs717302 (110)  A (9.1 %)  G (90.9 %)\n",
      "  rs719366 (51)  C (54.9 %)  T (45.1 %)\n",
      "  rs722098 (230)  A (99.1 %)\n",
      "  rs727811 (97)  A (7.2 %)  C (92.8 %)\n",
      "  rs729172 (222)  A (66.2 %)  C (33.8 %)\n",
      "  rs733164 (171)  A (21.6 %)  G (78.4 %)\n",
      "  rs735155 (146)  A (6.8 %)  G (93.2 %)\n",
      "  rs737681 (234)  C (100.0 %)\n",
      "  rs740910 (31)  A (0.0 %)  G (100.0 %)\n",
      "  rs763869 (158)  C (62.0 %)  T (38.0 %)\n",
      "  rs8037429 (66)  C (56.1 %)  T (43.9 %)\n",
      "  rs826472 (146)  C (100.0 %)\n",
      "  rs873196 (145)  C (17.9 %)  T (82.1 %)\n",
      "  rs876724 (117)  C (99.1 %)\n",
      "  rs891700 (39)  A (64.1 %)  G (35.9 %)\n",
      "  rs901398 (199)  C (58.3 %)  T (41.7 %)\n",
      "  rs907100 (153)  C (96.1 %)\n",
      "  rs914165 (156)  A (6.4 %)  G (93.6 %)\n",
      "  rs917118 (400)  C (100.0 %)\n",
      "  rs938283 (103)  C (72.8 %)  T (27.2 %)\n",
      "  rs964681 (106)  C (3.8 %)  T (96.2 %)\n"
     ]
    }
   ], 
   "source": [
    "# Get SNP profile\n",
    "\n",
    "snpData = {}\n",
    "\n",
    "with open(resultDir+'/SNPforID_sorted.bam.vcf') as f:\n",
    "    for l in f:\n",
    "        if l.startswith('#'):\n",
    "            continue\n",
    "            \n",
    "        snp, pos, id, ref, alt, qual, filter, info, d, dd = l.split()\n",
    "        \n",
    "        # Our SNP of interest is always at position 26 of the reference\n",
    "        if int(pos) != 26:\n",
    "            continue\n",
    "\n",
    "        par = {}\n",
    "        for p in info.split(';'):\n",
    "            pv = p.split('=')\n",
    "            par[pv[0]] = pv[1]\n",
    "        \n",
    "        snpData[snp] = {'pos': pos, 'ref': ref, 'alt': alt, 'qual': qual, 'filter': filter, 'info': par}\n",
    "\n",
    "# Display count\n",
    "print('Got data for {} SNPs:'.format(len(snpData)))\n",
    "\n",
    "# Save/print results\n",
    "with open(resultDir + '/ttt_profile.csv', 'w') as f:\n",
    "    # Table header\n",
    "    f.write('snp, coverage, ref_allele, ref_percent, alt_allele, alt_percent, genotype\\n')\n",
    "    \n",
    "    # Table data\n",
    "    for s in sorted(snpData.keys()):\n",
    "        totalDepth = int(snpData[s]['info']['DP'])\n",
    "        depthList  = [int(d) for d in snpData[s]['info']['DP4'].split(',')]\n",
    "        refDepth   = sum(depthList[0:2])\n",
    "        altDepth   = sum(depthList[2:4])\n",
    "        \n",
    "        # Estimate the diploid genotype: when the minor allele is more than 10 times weaker than the major allele,\n",
    "        # we should ignore it for a pure sample?\n",
    "        if refDepth > altDepth and altDepth/refDepth < 0.1:\n",
    "            genotype = snpData[s]['ref'] + snpData[s]['ref']\n",
    "        elif altDepth > refDepth and refDepth/altDepth < 0.1:\n",
    "            genotype = snpData[s]['alt'] + snpData[s]['alt']\n",
    "        else:\n",
    "            genotype = snpData[s]['ref'] + snpData[s]['alt']\n",
    "        \n",
    "        if snpData[s]['alt'] == '.':\n",
    "            # Only 1 allele was observed\n",
    "            f.write(','.join([s, str(totalDepth), snpData[s]['ref'], '{:.1f}'.format(100*refDepth/totalDepth), '', '', snpData[s]['ref']+snpData[s]['ref']]) + '\\n')\n",
    "            print('  {} ({})  {} ({:.1f} %)'.format(s, totalDepth, snpData[s]['ref'], 100*refDepth/totalDepth))\n",
    "        else:\n",
    "            # Two alleles were observed\n",
    "            f.write(','.join([s, str(totalDepth), snpData[s]['ref'], '{:.1f}'.format(100*refDepth/totalDepth), snpData[s]['alt'], '{:.1f}'.format(100*altDepth/totalDepth), genotype]) + '\\n')\n",
    "            print('  {} ({})  {} ({:.1f} %)  {} ({:.1f} %)'.format(s, totalDepth, snpData[s]['ref'], 100*refDepth/totalDepth, snpData[s]['alt'], 100*altDepth/totalDepth))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.3.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}