I am a biologist and use different packages like DESeq, ... to normalize my data and find deferential expressed genes.
Recently I have started to learn probability and statistics and I have studied distributions quite well. But I still have a problem: I think I do not very well understand why we really use this distributions to infer expression levels for genes, normalization, find differential expressed genes?

Why do we need e.g. Poisson model, negative binomial, ... for obtaining an approximate expression level? or in a package called mmseq: "Expression levels are inferred for each transcript using the mmseq program by modelling mappings of reads or read pairs (fragments) to sets of transcripts"!! why modeling? why do we need to estimate expression level while we can directly count the number of reads per gene?

Or why is it appropriate to model read counts as a e.g. Poisson process?


Is it only due to the fact that knowing the distribution (e.g negative binomial which can very well explain the observed counts, considering noise, ...) help us to apply the right properties like mean, var, ..., on data or there are more things to learn from the distributions?

Sorry if my question is primitive but it is a long time that I am struggling with that


2 Answers 2

  1. The data is count data because it's the number of counts aligned to a gene. It's not continuous and therefore can't be modelled as say a normal distribution.
  2. Poisson distribution is designed for modelling count data.
  3. However, the Poisson distribution assumes the first and second moments (mean and variance) are equal. This is not true for RNA-Seq. Lowly expressed genes have much higher variance than highly expressed genes.
  4. To account for the variability, we use the negative-binomial model which is really an extension of Poisson. The NB model has an extra parameter to model for the variance. It can be proven as the variance approach to the mean, the NB model becomes the Poisson model.


To answer your comments:

  1. Normalization is usually necessary to model the different sequencing depth between libraries. However, you don't need to do it yourself if you use DESeq2 or edgeR. They have their own normalization algorithm (Trim-mean-valued and upper-quartiles).

  2. Those packages do the normalization for you. Fit your data to a NB model, estimate the dispersion (ie: variance). Once they have a model, they can use whatever statistical method required (I think it's the Wald test for DESeq2) to check for differential expressed genes. The results depend on how much they express and also how much variance they have.

  • $\begingroup$ So, we first normalize the data with whatever approach and then fit it to Poisson/NB ... to find diff expressed genes? $\endgroup$
    – Sia_a
    Jul 10, 2016 at 13:22
  • $\begingroup$ @Sia_a EDITED. Short answer: YES, but you don't have to do it yourself. $\endgroup$
    – SmallChess
    Jul 10, 2016 at 13:26
  • $\begingroup$ Thanks for your answers. I routinely use DESeq2 for my stuff but I wanted to know more about what is going on under the scene. I was a bit confused how distribution can be used to model read count!! $\endgroup$
    – Sia_a
    Jul 10, 2016 at 13:28
  • $\begingroup$ @Sia_a DESeq is very hard to read. But just remember when an expressed gene is a gene that is sufficiently away from the fitted NB model. How to define "sufficient" statistically makes the package popular. Everyone has their own definitions, for example, I can say it should be 1 standard-deviation away. What makes DESeq2 popular is how they define their measures. $\endgroup$
    – SmallChess
    Jul 10, 2016 at 13:35
  • 1
    $\begingroup$ +1; very clear and concise. Just want to emphasize that the negative binomial modeling is required for getting estimates of the variance of the expression levels to allow significance testing. Expression levels could in principle just be estimated as the OP suggests, based on the number of reads mapped to a gene, and corrected for sequencing depth. But you wouldn't know how reliable the estimates would be. The modeling also allows information on variance to be pooled among all reads in parallel, avoiding the less reliable variance estimates that would come from analysis one gene at a time. $\endgroup$
    – EdM
    Jul 10, 2016 at 13:53

It is usual to model counts taking into account that:

1 - they only take on integer values
2 - they are always non-negative
3 - (most importantly) their variability increases with their mean

The negative binomial is used when the dispersion of the counts is greater than would be expected from the Poisson. In both cases it is usual to model using a log link.

  • $\begingroup$ Thanks for your reply. did I get it correctly: "once we have row counts, we first normalize the data, then, we fit our normalized data to a distribution like NB to use its properties like mean, var, ... to find differential expressed genes"? BUT I do not understand this sentence: "Poisson/NB model is often used to infer expression levels for genes using read counts." ? why infer?? while we can directly count the number of reads per gene? $\endgroup$
    – Sia_a
    Jul 10, 2016 at 12:02
  • $\begingroup$ @Sia_a The models are just estimation. We fit the data to the NB model and then calculate "distance" between a given gene with the NB model. The distance measure could be something like a null hypothesis that the gene is not differential expressed. Once we have the model, we can use it to compute Wald test or log-likehood ratio test. $\endgroup$
    – SmallChess
    Jul 10, 2016 at 13:18
  • $\begingroup$ @Sia_a I've edited my answer to respond to your now deleted comments. $\endgroup$
    – SmallChess
    Jul 10, 2016 at 13:25
  • $\begingroup$ @mdewey could you please comment on the concept of "variability increases with their (the counts) mean". This seems contradictory with the top answer that says: "Lowly expressed genes have much higher variance than highly expressed genes." $\endgroup$
    – Sergio.pv
    Jul 28, 2017 at 12:47

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.