Framing the negative binomial distribution for DNA sequencing

Question

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?

Answer 1

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.

Answer 2

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.

Answer 3

• 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.

Answer 4

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.

Answer 5

I'll try to give a simplistic mechanistic interpretation that I found useful when thinking about this.

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.