# Framing the negative binomial distribution for DNA sequencing

The negative binomial distribution has become a popular model for count data (specifically the expected number of sequencing reads within a given region of the genome from a given experiment) in bioinformatics. Explanations vary:

• Some explain it as something that works like the Poisson distribution but has an additional parameter, allowing more freedom to model the true distribution, with a variance not necessarily equal to the mean
• Some explain it as a weighted mixture of Poisson distributions (with a gamma mixing distribution on the Poisson parameter)

Is there a way to square these rationales with the traditional definition of a negative binomial distribution as modeling the number of successes of Bernoulli trials before seeing a certain number of failures? Or should I just think of it as a happy coincidence that a weighted mixture of Poisson distributions with a gamma mixing distribution has the same probability mass function as the negative binomial?

• It is also a compound Poisson distribution where you sum a Poisson-distributed number of logarithmic random variables. Sep 22, 2012 at 0:09

IMOH, I really think that the negative binomial distribution is used for convenience.

So in RNA Seq there is a common assumption that if you take an infinite number of measurements of the same gene in an infinite number of replicates then the true distribution would be lognormal. This distribution is then sampled via a Poisson process (with a count) so the true distribution reads per gene across replicates would be a Poisson-Lognormal distribution.

But in packages that we use such as EdgeR and DESeq this distribution modeled as a negative binomial distribution. This is not because the guys that wrote it didn't know about a Poisson Lognormal distribution.

It is because the Poisson Lognormal distribution is a terrible thing to work with because it requires numerical integration to do the fits etc. so when you actually try to use it sometimes the performance is really bad.

A negative binomial distribution has a closed form so it is a lot easier to work with and the gamma distribution (the underlying distribution) looks a lot like a lognormal distribution in that it sometimes looks kind of normal and sometimes has a tail.

But in this example (if you believe the assumption) it can't possibly be theoretically correct because the theoretically correct distribution is the Poisson lognormal and the two distributions are reasonable approximations of one another but are not equivalent.

But I still think the "incorrect" negative binomial distribution is often the better choice because empirically it will give better results because the integration performs slowly and the fits can perform badly, especially with distributions with long tails.

I looked through a few web pages and couldn't find an explanation, but I came up with one for integer values of $r$. Suppose we have two radioactive sources independently generating alpha and beta particles at the rates $\alpha$ and $\beta$, respectively.

What is the distribution of the number of alpha particles before the $r$th beta particle?

1. Consider the alpha particles as successes, and the beta particles as failures. When a particle is detected, the probability that it is an alpha particle is $\frac{\alpha}{\alpha+\beta}$. So, this is the negative binomial distribution $\text{NB}(r,\frac{\alpha}{\alpha+\beta})$.

2. Consider the time $t_r$ of the $r$th beta particle. This follows a gamma distribution $\Gamma(r,1/\beta).$ If you condition on $t_r = \lambda/\alpha$, then the number of alpha particles before time $t_r$ follows a Poisson distribution $\text{Pois}(\lambda).$ So, the distribution of the number of alpha particles before the $r$th beta particle is a Gamma-mixed Poisson distribution.

That explains why these distributions are equal.

• Some explain it as something that works like the Poisson distribution but has an additional parameter, allowing more freedom to model the true distribution, with a variance not necessarily equal to the mean
• Some explain it as a weighted mixture of Poisson distributions (with a gamma mixing distribution on the Poisson parameter)

Mathematically one obtains Negative binomial by integrating the Poisson distribution over Gamma-distributed weights, see Gamma-Poisson mixture. This mathematical fact remains regardless of whether we accept it as the justification for using the distribution or not.

Poisson distribution is a rather natural choice when talking about counting reads arising from DNA sequencing (one could use binomial, but given that one sequences only a small fraction of reads/DNA obtained from the sample, the difference is negligible, and we can use whatever seems more convenient.) We are also sure that the parameter of this Poisson distribution varies, although the reason for this variation depends on the exact nature of the experiment - e.g., it can be variation due to

• replicating the same experiment several times
• the reads originating from different cells with somewhat different properties
• comparing the numbers of reads corresponding different genes
• genes having different chemical structure and therefore amplified differently by PCR or some reads more likely to make their way to the sequencing machine
• the student/postdoc preparing the libraries not being very careful/consistent
• etc.

In other words, we are sure that the variation exists (and we do observe it experimentally), but we don't know exactly where it comes from, and we cannot directly know what probability distribution describes it. We couldn't model it using the normal distribution, since the Poisson parameter should be positive, so we use the Gamma distribution, because it is "almost like normal", but with non-negative support... but we could have also used log-normal or something else. As long as we are not looking for the fine biological effects that could turn out to be artifacts of the particular distribution we use, anything that is computationally convenient is good.

Note that, besides the flexibility provided by an extra parameter, negative binomial has a thicker tail than the Poisson distribution, making it less sensitive to outliers. This provides an additional motivation for using this distribution: it allows more robust inference.

I can only offer intuition, but the gamma distribution itself describes (continuous) waiting times (how long does it take for a rare event to occur). So the fact that a gamma-distributed mixture of discrete poisson distributions would result in a discrete waiting time (trials until N failures) does not seem too surprising. I hope someone has a more formal answer.

Edit: I always justified the negative binomial dist. for sequencing as follows: The actual sequencing step is simply sampling reads from a large library of molecules (poisson). However that library is made from the original sample by PCR. That means that the original molecules are amplified exponentially. And the gamma distribution describes the sum of k independent exponentially distributed random variables, i.e. how many molecules in the library after amplifying k sample molecules for the same number of PCR cycles.

Hence the negative binomial models PCR followed by sequencing.

• That makes sense, but in the context of measuring the number of sequencing reads in a genome is there an intuitive explanation for what the waiting period in the negative binomial distribution represents? In this case there is no waiting period - he's just measuring counts of sequencing reads. Sep 21, 2012 at 20:31
• See my edit. I don't see how thinking of it in terms of waiting times fits the sequencing setting. The gamma poisson mixture is easier to interpret. But in the end they are the same thing. Sep 21, 2012 at 20:34
• Ok - then perhaps the real question is by what coincidence does modeling k successes + r failures in Bernoulli trials follow a gamma Poisson mixture? Maybe a negative binomial modeling k successes + r failures can be thought of as an overdispersed Poisson dbn due to the many possible permutations of success and failure trials resulting in the exactly k observed successes and r observed failures, which can be described as a collection of separate dbns? Sep 21, 2012 at 21:02

Assume we have a perfect uniform coverage of the genome before library prep, and we observe $\mu$ reads covering a site on average. Say that sequencing is a process that picks an original DNA fragment, puts it through a stochastic process that goes through PCR, subsampling, etc, and comes up with a base from the fragment at frequency $p$, and a failure otherwise. If sequencing proceeds until $\mu\frac{1-p}{p}$ failures, it can be modeled with a negative binomial distribution, $NB(\mu\frac{1-p}{p}, p)$.
Calculating the moments of this distribution, we get expected number of successes $\mu\frac{1-p}{p}\frac{p}{1-p} = \mu$ as required. For variance of the number of successes, we get $\sigma^2 = \mu(1-p)^{-1}$ - the rate at which the library prep fails for a fragment increases the variance in the observed coverage.
While the above is a slightly artificial description of the sequencing process, and one could make a proper generative model of the PCR steps etc, I think it gives some insight into the origin of the overdispersion parameter $(1-p)^{-1}$ directly from the negative binomial distribution. I do prefer the Poisson model with rate integrated out as an explanation in general.