While reading old DESeq paper (Anders and Huber 2010) I came across following line.
If reads were independently sampled from a population with given, fixed fractions of genes, the read counts would follow a multinomial distribution, which can be approximated by the Poisson distribution.
I am unable to grasp this line. Why would read counts follow multinomial distribution? My guess would have been Gaussian distribution. Can anyone explain this?
Textbook example of multinomial distribution is multiple dice toss. The fair dice has following probabilities:
side 1: 16,66666666%
side 2: 16,66666666%
side 3: 16,66666666%
side 4: 16,66666666%
side 5: 16,66666666%
side 6: 16,66666666%
Lets roll a typical fair dice n = 20 times:
6 5 1 1 3 5 4 3 1 2 4 4 6 6 2 2 5 6 5 2
So this particular outcome of above multinomial variable is:
side 1: 3 rolls
side 2: 4 rolls
side 3: 2 rolls
side 4: 3 rolls
side 5: 4 rolls
side 6: 4 rolls
Multinomial distribution is not restricted to fair dice - the probabilities can be rigged. It also is not restricted to 6 sides - there can be any number of categories. Another textbook example of multinomial distribution is puling colored balls from wery (infinitely) large bag. Probabilities correspond to proportions of colored balls in the bag:
Blue: 53.283%
Green: 19.956%
Orange: 8.336%
Purple: 5.213%
Red: 4.374%
Silver: 3.920%
White: 2.751%
Yellow: 2.167%
Lets simulate n = 200 balls pulled from such a bag in R:
sample(c("B","G","O","P","R","S","W","Y"), replace = T, size = 200, prob = c(0.53283,0.19956,0.08336,0.05213,0.04374,0.03920,0.02751,0.02167))
B S O B B G G P B O B B P G B B G G G G B P G B B P B G B G B S W B B O O G B B O G G B S B O G B B B B B O O O B O B B B B B B O O G B B Y B R G B B B G P W Y G B P W B S R B W G Y B W B G O G B R B B G B B B B B B P R G P B B B G Y G S B B B G P B B B B B G B O G P B G B G P B B G B R B P R B G W B B B B O P G B B B B B B B B G B G G R W B B G G B G O B B B B B G B B B P B B B O R G B B B O B G
So this particular outcome multinomial variable is:
Blue: 102 balls
Green: 42 balls
Orange: 18 balls
Purple: 14 balls
Red: 8 balls
Silver: 5 balls
White: 7 balls
Yellow: 4 balls
In RNA-seq we are "pulling" reads (balls) out of a large set of suitable cDNA molecules fragments in a sample (bag). Each read belongs to a gene (color). We assume fixed fractions of genes in a sample:
gene 1: 0.05217%
gene 2: 0.00319%
gene 3: 0.00073%
...
But what we get from RNA-seq (n = milions of reads) are integer counts of reads:
gene 1: 492 reads
gene 2: 44 reads
gene 3: 5 reads
...
To comment on your idea of using Gausian distribution: Gausian distribution is continuous so using it we would assume we can get fractions of reads per gene which is not the case.
