{
"metadata": {
"name": "",
"signature": "sha256:3911bc4d6ad61bdfd7f9ebb93e6bc1128d38541e33d113272ffd4dbe26d78a62"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# EE126 Lab10: RNA sequencing\n",
"\n",
"The direct sequencing of RNA transcripts, known as RNA Sequencing (http://en.wikipedia.org/wiki/RNA-Seq), has many applications including genome annotation, comprehensive identification of fusions in cancer, discovery of novel isoforms of genes, and genome sequence assembly [Lior Pachter 2011] (http://arxiv.org/abs/1104.3889, http://genomemedicine.com/content/3/11/74/abstract). The problem of RNA sequencing is to figure out how much and what type of RNA is present in a genome at a given moment in time.\n",
"\n",
"For our purposes, we'll phrase the problem as follows: Given a set of short reads that are sampled from a set of larger genes, how can we find the relative abundance of each gene. That is, given just the short reads, how do we know how frequently each original gene occurs. This process is depicted in Figure 1. (Aside: in the actual paper, these \"genes\" are actually \"transcripts,\" but that's not relevant for us)\n",
"\n",
"\n",
"\n",
"####Figure 1: We want to use the reads to guess the underlying proportion.\n",
"\n",
"Turns out, we can use some methods we learned in class (MLE, EM) to solve this problem! Let's try it out.\n",
"\n",
"We assume that all genes are of the same length, $\\ell_t$, and all reads are of the same length, $\\ell_x$. Then, we assume the following Bayesian generative model:\n",
"\n",
"1) A read comes from a randomly chosen gene\n",
"\n",
"2) A read's starting point is randomly chosen among all possible $\\ell_t$ starting positions in that gene\n",
"\n",
"Again, our goal is to estimate the abundance of each gene."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's begin by coming up with a way to use the tools we've learned so far to tackle this task. First, given a set of reads $X = \\{X_1, X_2, \\ldots, X_n\\}$, we want to figure out what distribution over the genes $\\rho$ was most likely to give us $X$. To do this, all we need to do is maximize the following likelihood function:\n",
"\n",
"$ L(\\rho) = \\log{ P( X|\\rho) } = \\log{ \\prod_{i=1}^{N}{ P(X_i|\\rho) } } $\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#Simple model\n",
"\n",
"To make life easy, we's first going to assume no read is ambiguous\u2014given a read, we can immediately tell which gene it came from (we know the color coding of the reads in in figure 1)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###Q1. Given the assumption above, each $X_i$ is aligned with only one gene, $t_i$. What is $P(X_i|\\rho)$? Your solution should take both gene and read lengths into consideration.\n",
"Hint: How many possible starting positions exist for $X_i$? \n",
"\n",
"You may denote the probability of seeing a specific gene as $\\rho_{t_i}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###A1. Your Solution Here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###Q2. Assume that you have two genes, and $x$ reads are compatible with gene $1$, and $n-x$ reads are compatible with the other gene, gene $2$. Using your solution above, find the maximum likelihood estimator of $\\rho$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###A2. Your Solution Here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Q3. (Even more optional) Assume that you have $M$ genes. Out of $n$ reads, $x_i$ reads are compatible with gene $i$, where $\\sum_{i=1}^{M}{x_i} = n$. Find the maximum likelihood estimator of $\\rho$. \n",
"\n",
"You might just be able to make an \"educated guess\" from your answer above. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### A3. Your Solution Here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Harder Model\n",
"\n",
"Life gets harder when you allow for ambiguous reads. Going back to the figure above, we can imagine it as not knowing exactly which color a read belongs to, but we have a rough idea of the possible colors it could have been. See the modified figure below.\n",
"\n",
"\n",
"#### Figure 2 \n",
"In this portion, we'll consider a general problem where a read is aligned possibly more than one gene. We first define a compatibility matrix $A \\in \\{0,1\\}^{N \\times M} = \\{a_{i,j}\\}$, where $a_{i,j}$ is $1$ if read $i$ is aligned with gene $j$, $0$ otherwise. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###Q4. Let's assume we have $3$ genes and $5$ reads, aligned as follows:\n",
"(ref: pp16-17 Pachter 2011)\n",
"\n",
"\n",
"#### Figure 3 \n",
"###Find the compatibility matrix $A$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###A4. Your Solution Here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###Q5. Given an arbitrary compatibility matrix, write an expression to find the likelihood function for $\\rho$. \n",
"\n",
"You'll have to tweak your solution to the last portion by carefully considering how you can represent $P(X_i|\\rho)$ given the compatibility matrix."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###A5. Your Solution Here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###Q6. (Even more optional. Not crucial for rest.) Going back to the model with $3$ genes and $5$ reads in Q4, find the maximum likelihood estimator of $\\rho$.\n",
"\n",
"Use the likelihood function you defined above. Clever algebraic manipulations can be used to bring it down to a form that you can maximize easily (a maximization over a single variable instead of three)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###A6. Your Solution Here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Great! So now we were able to find the most likely $\\rho$ given ambiguous reads...in that one case. As you probably realized when doing the problem above, things can get ugly really fast. If you didn't do the problem above, the solutions will show you how bad it got. Let's look at another way to estimate the most likely distribution instead of trying to directly calculate it.\n",
"\n",
"\n",
"#EM Algorithm! :D\n",
"\n",
"In general, it is not easy to find the exact maximum likelihood estimator of $\\rho$ if ambiguous reads exist. Instead, we can rely on iterative methods that we hope will converge to the true value. One way to go about this is via expectation maximization (EM), as you have seen in class. Here is an EM algorithm that can be used for this specific problem. You'll be implementing it shortly, so read it carefully and understand why it works. \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"####1. Initialize $\\rho = (\\frac{1}{M}, \\frac{1}{M}, \\ldots, \\frac{1}{M})$, all genes are equally probable.\n",
"####2. Find the compatibility matrix $A$.\n",
"####3. Repeat the following until $\\rho$ converges:\n",
"####3-1. For each gene $i$, find non-zero columns of $i$th row of $A$. Call these $\\mathcal{I}_i$.\n",
"####3-2. Choose the values of $\\rho$ in $\\mathcal{I}_i$, normalize the vector, and replace $i$th row of $A$ with the normalized vector.\n",
"####3-3. Replace $\\rho$ such that $\\rho_i = \\frac{1}{N}\\sum_{j=1}^{N}{A_{i,j}}$ \n",
"\n",
"The following figure is a visual representation of the algorithm (Ref: p17 Pachter 2011). It's a lot to digest at once, so only look at the first step to begin with.\n",
"\n",
"There are three possible genes to analyze (red, green, blue). There are five reads (a,b,c,d,e) that you can use to help you. One maps to all three, one only to red, and the other three to each of the three pairs. Initially, you assume a uniform prior (rho_guess).\n",
"\n",
"During the expectation (E) step, reads are proportionately assigned to each of the genes they could have come from according to their abundances (rho_guess, currently uniform). \n",
"\n",
"Next, during the maximization (M) step, the abundances are recalculated from the proportionately assigned read counts. \n",
"\n",
"For example, the abundance of red after the first M step is estimated by\n",
"0.47 = (0.33 + 0.5 + 1 + 0.5)/(2.33 + 1.33 + 1.33)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"#### Figure 4 "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Q6. Explain the above EM algorithm in your own words. Will it always find the maximum likelihood estimator of $\\rho$? \n",
"Hint: Is the likelihood function convex? concave?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###A6. Your Solution Here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###Q7. Implement the above EM algorithm. Run it with the given set of reads & genes. What is your estimated $\\rho$?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###A7. Complete the skeleton provided below. How close do you get? :D\n",
"\n",
"Feel free to alter any of the code in the fourth code box."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from __future__ import division\n",
"s = ['ATCTCGACGCACTGC', 'GAGTTCGAACTCTTC', 'AGAGTTCCAGTGTCA', 'AAAGCTCACTGCGGA', 'AGCGATATCAGAGTD']\n",
"M = len(s) # Number of transcripts"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"rho = rand(M)\n",
"rho /= sum(rho)\n",
"print rho\n",
"# This rho is unknown to us"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def rand_choice( pmf ):\n",
" v = rand()\n",
" tmp = 0\n",
" for i in range(len(pmf)):\n",
" each_val = pmf[i]\n",
" tmp += each_val\n",
" if v <= tmp:\n",
" return i\n",
"\n",
"def random_read( s, rho, L ):\n",
" chosen_seq = s[rand_choice( rho )]\n",
" start_idx = randint( len(chosen_seq) - L )\n",
" end_idx = start_idx + L\n",
" return chosen_seq[ start_idx:end_idx ]\n",
" \n",
"N = 1000 # Number of reads\n",
"L = 5\n",
"\n",
"reads = []\n",
"for i in range(N):\n",
" reads.append( random_read(s, rho, L) )\n",
" \n",
"print 'First 20 reads...', reads[0:20] "
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"N_iter = 100 #number of E/M iterations\n",
"\n",
"def find_all_alignments(s, read):\n",
" tmp = []\n",
" for j in range(len(s)):\n",
" if read in s[j]:\n",
" tmp.append(j)\n",
" return tmp\n",
"\n",
"# A: Compatibility Matrix. \n",
"#Note: You can choose to represent this matrix in whatever form is easiest for your calculations\n",
"#The form we are pushing you towards here is just how we chose to implement the algorithm\n",
"A = []\n",
"for each_read in reads:\n",
" A.append( find_all_alignments(s, each_read) )\n",
"\n",
"hidden_state_prior = zeros([N, M])\n",
"rho_est = []\n",
"\n",
"#Your EM code here. Print your final answer as rho_est\n",
"\n",
"print '(real) rho', rho\n",
"print '(est.) rho', rho_est"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#Congratulations! You've just finished your last lab in 126! Take a minute to appreciate what you've accomplished through the last semester. This class wasn't easy, but hopefully you've learned a lot (and let us know if you didn't >.<). Come see us if you'd like some help figuring out where to go from here. \n",
"\n",
"#Happy Holidays!"
]
}
],
"metadata": {}
}
]
}