Answers.Body: <p>I actually disagree with that.. I guess I write this down as a discussion.</p>
  <p>The expected value for an average gene is 1200. For this gene at 1/1200 expression, you expect 1 read.</p>
  <p>However, because of sampling, sometimes you get 0, 1 or 2. If we assume the sampling follows a poisson distribution, we can calculate the probability of getting 0, 1 , 2 ... reads:</p>
  <pre><code>plot(dpois(0:10,lambda=1),ylab=&quot;Probability&quot;,xlab=&quot;No of reads&quot;)
  </code></pre>
  <p><a href="https://i.sstatic.net/R3q8V.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/R3q8V.png" alt="enter image description here" /></a></p>
  <p>The probabilities are like:</p>
  <pre><code>data.frame(n=0:5,p=dpois(0:5,lambda=1))
  
    n           p
  1 0 0.367879441
  2 1 0.367879441
  3 2 0.183939721
  4 3 0.061313240
  5 4 0.015328310
  6 5 0.003065662
  </code></pre>
  <p>The probability of getting exactly 0 or 1 read should be ~0.37. The probability of getting at least 1 read would be 1 - P(n=0) = 1- 0.367879441 = 0.6321206</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-19T14:20:05.340Z
Answers.Text: I guess the author just assumes uniform distribution instead of poisson.
Answers.UserDisplayName: haci
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-19T14:20:38.690Z
Answers.Text: possible.. but why 50:50 ?
Answers.UserDisplayName: StupidWolf
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-19T15:47:46.373Z
Answers.Text: I was curious about this question but didn't know what the distribution would be, so I left a very clumsy monte-carlo simulation running yesterday.  After about 8,500 iterations I get exactly these values.
Answers.UserDisplayName: Jesse
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-19T15:57:52.867Z
Answers.Text: One other thing: Could the author have taken the probability determined for the lambda=1200 value of getting >=1200 (about 50%) and incorrectly assumed it'd stay the same for lambda=1?
Answers.UserDisplayName: Jesse
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-19T19:29:41.343Z
Answers.Text: There are a lot of back-of-the-envelope calculations in the original statement. I don't think 50:50 is supposed to be an exact number. The point is that you cannot be very confident in detecting that gene.
Answers.UserDisplayName: burger
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-19T14:06:02.783Z
Answers.LastActivityDate: 2020-11-19T14:06:02.783Z
Answers.UserDisplayName: StupidWolf
Body: <p>I am reading a book about RNAseq analysis and it says</p>
  <blockquote>
  <p>&quot;To calculate the probability that a read will map to a specific gene, we can assume an average gene size of 4000 nt (100 M nt divided by 25,000 genes). At 30 M reads equivalent to 30× coverage, at single read 100 nt (or paired-end read 50 nt) length, we can expect a single read to map to the average expressed and length gene 4000 nt× 30 coverage/100 nt 1200 times. Thus, if the gene is expressed at a level of 1/1200 compared to the average gene, then we have a 50:50 probability to have a read map to it.&quot;</p>
  </blockquote>
  <p>I understand that a single read would map to the average expressed 1200 times but I don't understand how the 50:50 probability is calculated. Please could you explain how this probability is determined.</p>
  
Tags: <rna-seq><reads><coverage>
Title: If a gene is expressed at a level of 1/1200 compared to the average gene, how is probability 50:50 that we have a read mapped to it?
ThreadText: # If a gene is expressed at a level of 1/1200 compared to the average gene, how is probability 50:50 that we have a read mapped to it?
  
  (Asked by: username_0 on 2020-11-18)
  
  I am reading a book about RNAseq analysis and it says
  
  > "To calculate the probability that a read will map to a specific gene, we can assume an average gene size of 4000 nt (100 M nt divided by 25,000 genes). At 30 M reads equivalent to 30× coverage, at single read 100 nt (or paired-end read 50 nt) length, we can expect a single read to map to the average expressed and length gene 4000 nt× 30 coverage/100 nt 1200 times. Thus, if the gene is expressed at a level of 1/1200 compared to the average gene, then we have a 50:50 probability to have a read map to it."
  
  I understand that a single read would map to the average expressed 1200 times but I don't understand how the 50:50 probability is calculated. Please could you explain how this probability is determined.
  
  **Tags:** rna-seq, reads, coverage
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2020-11-19. Score: 3
  
  I actually disagree with that.. I guess I write this down as a discussion.
  
  The expected value for an average gene is 1200. For this gene at 1/1200 expression, you expect 1 read.
  
  However, because of sampling, sometimes you get 0, 1 or 2. If we assume the sampling follows a poisson distribution, we can calculate the probability of getting 0, 1 , 2 ... reads:
  
  ```
  plot(dpois(0:10,lambda=1),ylab="Probability",xlab="No of reads")
  
  ```
  
  [![enter image description here](https://i.sstatic.net/R3q8V.png)](https://i.sstatic.net/R3q8V.png)
  
  The probabilities are like:
  
  ```
  data.frame(n=0:5,p=dpois(0:5,lambda=1))
  
    n           p
  1 0 0.367879441
  2 1 0.367879441
  3 2 0.183939721
  4 3 0.061313240
  5 4 0.015328310
  6 5 0.003065662
  
  ```
  
  The probability of getting exactly 0 or 1 read should be ~0.37. The probability of getting at least 1 read would be 1 - P(n=0) = 1- 0.367879441 = 0.6321206
  
  > **Comments:**
  > **username_2**: I guess the author just assumes uniform distribution instead of poisson.  
  > **username_1**: possible.. but why 50:50 ?  
  > **username_3**: I was curious about this question but didn't know what the distribution would be, so I left a very clumsy monte-carlo simulation running yesterday.  After about 8,500 iterations I get exactly these values.  
  > **username_3**: One other thing: Could the author have taken the probability determined for the lambda=1200 value of getting >=1200 (about 50%) and incorrectly assumed it'd stay the same for lambda=1?  
  > **username_4**: There are a lot of back-of-the-envelope calculations in the original statement. I don't think 50:50 is supposed to be an exact number. The point is that you cannot be very confident in detecting that gene.  
  
  ---
  
  