ipython.notebooks/Cell2015_N0B_A5SS_Fastq_to_Spliced_Reads.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy\n",
    "import scipy.sparse\n",
    "import scipy.stats\n",
    "import os\n",
    "import scipy.io as sio\n",
    "import dnatools\n",
    "from collections import Counter\n",
    "%matplotlib inline\n",
    "from pylab import *\n",
    "# Plotting Params:\n",
    "rc('mathtext', default='regular')\n",
    "fsize=14"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting a mapping between intronic sequences and 3' UTR barcodes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This ipython notebook first gets all of the intronic sequence to 3'UTR barcode mappings. It then counts the number of spliced reads for each sequence at each position. Some important information:<br>\n",
    "1. We only want to look at reads with the first 6 nt of the index read matching 'CACTGT'<br>\n",
    "2. Read2 reads the 3'UTR barcode and then into the fixed sequence after<br>\n",
    "NNNNNNNNNNNNNNNNNNNNGCAGGTAATG<br>\n",
    "I only count reads that at least include caggtaat in the first 10nt. This is a little less stringent than requiring a perfect match.<br>\n",
    "3. Read 1 actually starts 1nt past the exon-exon junction<br>\n",
    "Plasmid reads will look like this:<br>\n",
    "TGCTTGGNNNNNNNNNNNNNNNNNNNNNNNNNGGTCGACCCAGGTTCGTGNNNNNNNNNNNNNNNNNNNNNNNNNGAGGTATTCTTATCACCTTCGTGGCT<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "resultsdir = '../results/N0_A5SS_Fastq_to_Splice_Reads/'\n",
    "if not os.path.exists(resultsdir):\n",
    "    os.makedirs(resultsdir)\n",
    "figdir = '../figures/N0_A5SS_Fastq_to_Splice_Reads/'\n",
    "if not os.path.exists(figdir):\n",
    "    os.makedirs(figdir)\n",
    "    \n",
    "#Choose if you want to actually save the plots:\n",
    "SAVEFIGS = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 1 | 1000000 347487 |\n"
     ]
    }
   ],
   "source": [
    "f = {}\n",
    "f[0] = open('../fastq/A5SS_dna_R1.fq','r')\n",
    "f[1] = open('../fastq/A5SS_dna_R2.fq','r')\n",
    "tagcount = {}\n",
    "count = 0\n",
    "tagqual = {}\n",
    "header = {}\n",
    "seq = {}\n",
    "strand = {}\n",
    "quality ={}\n",
    "while True:\n",
    "    for i in range(2):\n",
    "        header[i] = f[i].readline()[:-1]\n",
    "        seq[i] = f[i].readline()[:-1]\n",
    "        strand[i] = f[i].readline()[:-1]\n",
    "        quality[i] = f[i].readline()[:-1]\n",
    "    tag = dnatools.reverse_complement(seq[1][:30])\n",
    "\n",
    "    if len(quality[1])==0:\n",
    "        break # End of File\n",
    "    if (tag[:10].count('ATTACCTG')>0):\n",
    "        qual = np.array([ord(i)-66 for i in quality[1][:20]])\n",
    "        try:\n",
    "            tagcount[tag] = tagcount[tag] + 1.\n",
    "        except:\n",
    "            tagcount[tag] = 1\n",
    "            tagqual[tag] = qual\n",
    "        else:\n",
    "            tagqual[tag]+=qual\n",
    "    if(count%1000000)==0:\n",
    "        print count,len(tagcount),'|',\n",
    "    count = count +1\n",
    "f[0].close()\n",
    "f[1].close()\n",
    "\n",
    "tagqual = pd.Series(tagqual)\n",
    "tagcount = pd.Series(tagcount)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get sequences for tags:\n",
    "If a sequence for a tag is observed twice it is assumed to be the true mapping.  For tags, which have no\n",
    "sequence found twice, the highest quality read sequence is chosen as the mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1810.0"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tagcount.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filter 3' barcode tag sequences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I will only keep barcodes that have a minimum average pred of 21 or more. I will also only keep barcodes that were observed at least twice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "filter_thresh = 21\n",
    "count_thresh = 2\n",
    "DNA_index='CACTGT'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of Tags 304288\n"
     ]
    }
   ],
   "source": [
    "#Filter tags\n",
    "# I summed the tag quality at every base, so to get the average, divide\n",
    "# by the number of tag counts.\n",
    "tag_qual_min = (tagqual/tagcount).apply(min)\n",
    "\n",
    "tags_filtered = tag_qual_min[(tag_qual_min>=filter_thresh)&(tagcount>=count_thresh)]\n",
    "print 'Number of Tags', len(tags_filtered)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a map between 3' barcodes and Intronic Sequences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the mapping between intronic randomized sequences and the corresponding 3' barcodes, we will find the first sequence that occurs twice with the same barcode. This might not be the optimal solution, but it is efficient. For the barcodes that do not have an intronic sequence occuring together twice, we will take the sequence with the highest minimum base quality score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tag_seqs = {}\n",
    "mapped_tag2seq = {}\n",
    "for tag in tags_filtered.index:\n",
    "    tag_seqs[tag] = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000000 182443 | 302519\n"
     ]
    }
   ],
   "source": [
    "f = {}\n",
    "f[0] = open('../fastq/A5SS_dna_R1.fq','r')\n",
    "f[1] = open('../fastq/A5SS_dna_R2.fq','r')\n",
    "count = 0\n",
    "header = {}\n",
    "seq = {}\n",
    "strand = {}\n",
    "quality ={}\n",
    "while True:\n",
    "    for i in range(2):\n",
    "        header[i] = f[i].readline()[:-1]\n",
    "        seq[i] = f[i].readline()[:-1]\n",
    "        strand[i] = f[i].readline()[:-1]\n",
    "        quality[i] = f[i].readline()[:-1]\n",
    "    if len(quality[1])==0:\n",
    "        break # End of File\n",
    "    tag = dnatools.reverse_complement(seq[1][:30])\n",
    "    try:\n",
    "        mapped_tag2seq[tag]\n",
    "    except:\n",
    "        try:\n",
    "            tag_seqs[tag][seq[0]]\n",
    "        except:\n",
    "            try:\n",
    "                tag_seqs[tag]\n",
    "            except:\n",
    "                pass\n",
    "            else:\n",
    "                tag_seqs[tag][seq[0]] = quality[1]\n",
    "        else:\n",
    "            mapped_tag2seq[tag] = seq[0]\n",
    "            del tag_seqs[tag]\n",
    "\n",
    "    count = count + 1\n",
    "    if (count%1000000)==0:\n",
    "        print count,len(mapped_tag2seq),'|',\n",
    "\n",
    "f[0].close()\n",
    "f[1].close()\n",
    "\n",
    "# Map tags that did not have any sequences with two observations based\n",
    "# on minimum base quality score:\n",
    "for tag in tag_seqs:\n",
    "    max_min_qual = 0\n",
    "    for seq in tag_seqs[tag]:\n",
    "        cur_min_qual = min([ord(s)-66 for s in tag_seqs[tag][seq]])\n",
    "        if cur_min_qual>max_min_qual:\n",
    "            mapped_tag2seq[tag] = seq\n",
    "            max_min_qual = cur_min_qual\n",
    "print len(mapped_tag2seq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tag2seqs = pd.Series(mapped_tag2seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove sequences with . or N:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tag2seqs = tag2seqs[~(tag2seqs.str.contains('\\.') | tag2seqs.str.contains('N'))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove sequences with deletions or insertions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tag2seqs = tag2seqs[tag2seqs.apply(lambda s:(s[32:49]=='GGTCGACCCAGGTTCGT')&(s[:7]=='TGCTTGG')&(s[75:100]=='GAGGTATTCTTATCACCTTCGTGGC') & (not ('N' in s)))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So unfortunately I did not save the tag2seqs in alphabetical order. Subsequent analysis must be done with my ordering. The barcodes and intronic sequences are identical to the sequences produced here, just in a different order. There is no way to reproduce this ordering, because it came from saving a dictionary. Instead we will just load my file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tag2seqs = pd.read_csv('../data/A5SS_Seqs.csv',index_col=0).Seq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Counting the Spliced Reads at Each Position For Each Plasmid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This contains the full intronic sequence plus part of the second exon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "intronseq = 'TGCTTGGNNNNNNNNNNNNNNNNNNNNNNNNNGGTCGACCCAGGTTCGTGNNNNNNNNNNNNNNNNNNNNNNNNNGAGGTATTCTTATCACCTTCGTGGCTACAGAGTTTCCTTATTTGTCTCTGTTGCCGGCTTATATGGACAAGCATATCACAGCCATTTATCGGAGCGCCTCCGTACACGCTATTATCGGACGCCTCGCGAGATCAATACGTATACCAGCTGCCCTCGATACATGTCTTGGCATCGTTTGCTTCTCGAGTACTACCTGGTTCCTCTTCTTTCTTTCTCTTCTCTTTCAGGACGGCAGCGTGCAGCTCGCCGACCACTACCAGCAGAACACCCCCATCGGCGACGGCCCCGTGCTGCTGCCCGACAACCACTACCTGAGCTACCAGTCCGCCCTGAGC'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To map the position of the spliced exon-exon junction, I simple map the last 20 nt of the read. For example, if the read was spliced at the first SD, then the 81-101 nt of the read will map 81-101 nt into the second exon. If the read was spliced 44nt 3' of the first SD, the 81-101 nt of the read will map 37-57 nt into the second exon. Unspliced reads will have the 81-101 nt of the read mapping within the intron. However, requiring an exact match for these 20 nt is very stringent and we may lose reads. So I will allow one mismatch. To do this efficiently, I precompute a mapping between all potential 1nt errors and the position the read should map to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Make a map of all 1nt mismatches to actuall positions in intron\n",
    "intronseq_map = {}\n",
    "bases = ['A','T','C','G']\n",
    "for b in xrange(len(intronseq)):\n",
    "    seq = intronseq[b:b+19]\n",
    "    intronseq_map[seq]=b\n",
    "    for pos in range(19):\n",
    "        for b1 in bases:\n",
    "            mut_seq = seq[:pos]+b1+seq[pos+1:]\n",
    "            intronseq_map[mut_seq] = b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This computes the position of the exon junction as described above. I generate the read count matrix of plasmids (rows) and exon junction positions (columns) by using scipy's sparse matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0 | 1000000 638192 | 2000000 1259923 | 3000000 1864253 | 4000000 2465283 | 5000000 3093399 | 6000000 3696937 | 7000000 4278083 | 8000000 4887881 | 9000000 5512395 | 10000000 6115250 | 11000000 6694008 | 12000000 7324952 | 13000000 7943569 | 14000000 8543184 | 15000000 9135218 | 16000000 9764341 | 17000000 10369947 | 18000000 10949916 | 19000000 11507487 | 20000000 12134084 | 21000000 12734212 | 22000000 13309266 |\n"
     ]
    }
   ],
   "source": [
    "tag_index = dict(zip(tag2seqs.index,range(len(tag2seqs))))\n",
    "\n",
    "bases = ['A','T','C','G']\n",
    "watsoncrick = {'N':'N','.':'.','C':'G','G':'C','A':'T','T':'A'}\n",
    "\n",
    "\n",
    "f = {}\n",
    "f[0] = open('../fastq/A5SS_rna_R1.fq','r')\n",
    "f[1] = open('../fastq/A5SS_rna_R2.fq','r')\n",
    "tag_list = []\n",
    "ss_list = []\n",
    "count=0\n",
    "used_reads=0\n",
    "header = {}\n",
    "seq = {}\n",
    "strand = {}\n",
    "quality ={}\n",
    "while True:\n",
    "    for i in range(2):\n",
    "        header[i] = f[i].readline()[:-1]\n",
    "        seq[i] = f[i].readline()[:-1]\n",
    "        strand[i] = f[i].readline()[:-1]\n",
    "        quality[i] = f[i].readline()[:-1]\n",
    "    tag = dnatools.reverse_complement(seq[1][:30])    \n",
    "    found = False\n",
    "    if len(quality[1])==0:\n",
    "        break # End of File\n",
    "\n",
    "    try:\n",
    "        plasmid_num = tag_index[tag]\n",
    "    except:\n",
    "        pass\n",
    "\n",
    "    else:\n",
    "        found = True\n",
    "    if found:\n",
    "        last20 = str(seq[0][82:101])\n",
    "        try:\n",
    "            splice_pos = intronseq_map[last20]\n",
    "        except:\n",
    "            pass\n",
    "        else:\n",
    "            splice_pos = 385-splice_pos #get splice position\n",
    "            if(splice_pos>=0)&(splice_pos<304):\n",
    "                tag_list.append(tag_index[tag])\n",
    "                ss_list.append(splice_pos)\n",
    "                used_reads += 1\n",
    "\n",
    "    if ((count%1000000)==0):\n",
    "        print count,used_reads,'|',\n",
    "    count += 1\n",
    "f[0].close()\n",
    "f[1].close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make the sparse matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "splices = {'A5SS':scipy.sparse.csr_matrix((list(np.ones_like(ss_list))+[0],\n",
    "                                           (tag_list+[len(tag2seqs)-1],ss_list+[303])))}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sio.savemat('../data/A5SS_Reads.mat',splices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filter reads that do not have the canonical GT/GC splice donor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know that splice donor have either GT or GC in the +1 to +2 positions (in the intron). Other \"spliced\" reads are probably noise. So we will enforce this condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 10000 20000 30000 40000 50000 60000 70000 80000 90000 100000 110000 120000 130000 140000 150000 160000 170000 180000 190000 200000 210000 220000 230000 240000 250000 260000\n"
     ]
    }
   ],
   "source": [
    "GT_GC_matrix = np.zeros_like(np.array(splices['A5SS'].todense()))\n",
    "intronic_sequences = tag2seqs.values\n",
    "for i in range(len(tag2seqs)):\n",
    "    cur_seq = 'G'+intronic_sequences[i]\n",
    "    for j in range(len(cur_seq)):\n",
    "        cur_dinuc = cur_seq[j:j+2]\n",
    "        if (cur_dinuc=='GT') | (cur_dinuc=='GC'):\n",
    "            GT_GC_matrix[i,j]=1\n",
    "    # Allow no splicing:\n",
    "    GT_GC_matrix[i,-1]=1\n",
    "    if(i%10000)==0:\n",
    "        print i,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "splices_GC_GT = scipy.sparse.csr_matrix((np.array(splices['A5SS'].todense())*GT_GC_matrix),\n",
    "                                        dtype=np.float64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save final matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combine with the A3SS data into a single file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Reads = sio.loadmat('../data/Reads.mat')\n",
    "Reads['A5SS'] = splices_GC_GT\n",
    "sio.savemat('../data/Reads.mat',Reads)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}