# Negative binomial modeling of RNA-Seq data

A common way to model RNA-seq data is using a negative binomial distribution, where each sample-gene pair is modeled by a different negative binomial distribution with mean $\mu_{ij}$ where $i$ and $j$ are indices for genes and samples, respectively (see article). My question is, given we only have a single observation ($X_{ij}$) from each sample-gene pair, how come we can learn the $\mu_{ij}$ from that single observation? We need multiple observations to estimate a mean, right? Maybe I am missing something critical.

• The negative binomial distribution is not fitted for a single gene-sample, but along a sample
– llrs
Feb 26 '18 at 8:05
• “how come we can learn the $\mu_{ij}$ from that single observation?” — we can’t. This fundamentally doesn’t work. Feb 26 '18 at 10:36
• @Llopis thanks, but if they are fitting for each gene separately, where is the sample-specific information coming from such that they can learn a different $\mu_{ij}$ for each pair? Feb 26 '18 at 13:43
• I wasn't sure of my first comments and I'm not with this one but, the value for each gene can't be estimated for a single sample. But it is estimated from all the samples. See the first section of the methods in the linked paper.
– llrs
Feb 26 '18 at 14:55
• @TomKelly It is, in fact, extremely common. The two most commonly used packages for differential expression use it (edgeR, DESeq[2]). Jun 18 '19 at 11:34

The $\mu_{ij}$ in that section of the DESeq manuscript is the expected value of sample $j$ in gene $i$ given its group association (with expected value $\mu_{ij}$). This is why $K_{ij}$ (this is presumably what you meant by $X_{ij}$) is a function of $\mu_{ij}$, rather than $\mu_{ij}$ itself. I think you're just swapping $K_{ij}$ and $\mu_{ij}$ in your head.

• "given its group association" -- what is the group here? a gene? Feb 26 '18 at 13:40
• @user5054 An experimental group (e.g., WT or treatment). Feb 26 '18 at 13:43
• Where exactly is this group modeled? Feb 26 '18 at 13:45
• @user5054 you can’t do differential expression without groups or covariates. Feb 26 '18 at 17:45
• You can however estimate the parameters of a negative binomial model by treating all samples as if they come from a single group. Mar 28 '18 at 17:53

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. However, the Poisson distribution assumes the mean equals the variance, which is usually not satisfied by the data. And for the RNA sequencing data, the variance is usually larger than the mean. So we can assume the data follows a $$Y|\lambda\sim Poisson(\lambda), \lambda\sim Gamma(r_0,b_0)$$ hierarchy.

\begin{align} P(Y=y)&=\int_{0}^{\infty}P(Y=y|\lambda)f(\lambda)d\lambda\\ &=\int_{0}^{\infty}\frac{e^{-\lambda}\lambda^y}{y!}\frac{b_0^{r_0}}{\Gamma(r_0)}\lambda^{r_0-1}e^{-b_0\lambda}d\lambda\\ &=\frac{b_0^{r_0}}{\Gamma(r_0)y!}\int_{0}^{\infty}e^{-\lambda}\lambda^y\lambda^{r_0-1}e^{-b_0\lambda}d\lambda\\ &=\frac{b_0^{r_0}}{\Gamma(r_0)y!}\int_{0}^{\infty}e^{-\lambda(b_0+1)}\lambda^{y+r_0-1}d\lambda\\ &=\frac{\Gamma(r_0+y)}{\Gamma(r_0)y!}\frac{b_0^{r_0}}{(b_0+1)^{r_0+y}}\int_{0}^{\infty}\frac{1}{\Gamma(r_0+y)}e^{-\lambda(b_0+1)}\lambda^{y+r_0-1}(b_0+1)^{r_0+y}d\lambda\\ &=\frac{\Gamma(r_0+y)}{\Gamma(r_0)y!}\frac{b_0^{r_0}}{(b_0+1)^{r_0+y}}\int_{0}^{\infty}\frac{1}{\Gamma(r_0+y)}e^{-\lambda(b_0+1)}[\lambda(b_0+1)]^{r_0+y}\frac{1}{\lambda}d\lambda\\ &=\frac{\Gamma(r_0+y)}{\Gamma(r_0)y!}\frac{b_0^{r_0}}{(b_0+1)^{r_0+y}}\\ &=\frac{\Gamma(r_0+y)}{\Gamma(r_0)y!}\left(\frac{1}{b_0+1}\right)^y\left(\frac{b_0}{b_0+1}\right)^{r_0}\\ \end{align} which is negative binomial distribution.

Negative binomial distribution depicts the number of failures $$k$$ before $$r$$th success of a series of independent Bernoulli experiments, with success probability $$p$$: $$f(k;r,p)\equiv Pr(X=k)=\frac{\Gamma(k+r)}{k!\Gamma(r)}p^r(1-p)^k\quad k=0,1,2,\cdots$$

$$\mu=\mathbb E(X)=\frac{r(1-p)}{p}$$ $$\sigma^2=\mathbb V(X)=\frac{r(1-p)}{p^2}$$

Then $$p=\frac{\mu}{\sigma^2}$$ $$r=\frac{\mu^2}{\sigma^2-\mu}$$ The probability mass function (PMF) is:

$$f(k;r,p)\equiv Pr(X=k)={k+\frac{\mu^2}{\sigma^2-\mu}-1\choose{\frac{\mu^2}{\sigma^2-\mu}-1}}\left(\frac{\sigma^2-\mu}{\sigma^2}\right)^k\left(\frac{\mu}{\sigma^2}\right)^{\frac{\mu^2}{\sigma^2-\mu}}$$

The $$\mu$$ and $$\sigma^2$$ of each gene can be estimated from data.

Assume that the number of reads in sample $$j$$ that are assigned to gene $$i$$ can be modeled by a negative binomial (NB) distribution,

$$K_{ij}\sim NB(\mu_{ij},\sigma^2_{ij})\tag1$$ which has two parameters, the mean $$\mu_{ij}$$ and the variance $$\sigma^2_{ij}$$. The read counts $$K_{ij}$$ are non-negative integers. $$q_{i, \rho(j)}$$ is proportional to the expectation value of the true (but unknown) concentration of fragments from gene $$i$$ under condition $$\rho(j)$$. The size factor $$s_j$$ represents the coverage, or sampling depth, of library $$j$$.

• First, the mean parameter $$\mu_{ij}$$ , that is, the expectation value of the observed counts for gene $$i$$ in sample $$j$$, is the product of a condition-dependent per-gene value $$q_{i, \rho(j)}$$(where $$\rho(j)$$ is the experimental condition of sample $$j$$) and a size factor $$s_j$$,

$$\mu_{ij} = q_{i, \rho(j)}s_j\tag2$$

• Second, the variance $$\sigma^2_{ij}$$ is the sum of a shot noise term and a raw variance term,

$$\sigma^2_{ij}=\mu_{ij}+s_j^2v_{i,\rho(j)}\tag3$$

• Third, we assume that the per-gene raw variance parameter $$v_{i,\rho(j)}$$ is a smooth function of $$q_{i, \rho(j)}$$,

$$v_{i,ρ(j)}=v_ρ(q_{i, \rho(j)}) \tag4$$

This assumption is needed because the number of replicates is typically too low to get a precise estimate of the variance for gene $$i$$ from just the data available for this gene. This assumption allows us to pool the data from genes with similar expression strength for the purpose of variance estimation (see article).