2
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ I guess the author just assumes uniform distribution instead of poisson. $\endgroup$ – haci Nov 19 at 14:20
  • $\begingroup$ possible.. but why 50:50 ? $\endgroup$ – StupidWolf Nov 19 at 14:20
  • 2
    $\begingroup$ 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. $\endgroup$ – Jesse Nov 19 at 15:47
  • $\begingroup$ 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? $\endgroup$ – Jesse Nov 19 at 15:57
  • $\begingroup$ 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. $\endgroup$ – burger Nov 19 at 19:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.