Negative binomial modeling of RNA-Seq data

Question

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.

Answer 1

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

Answer 2

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.