<p>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 <span class="math-container">$Y|\lambda\sim Poisson(\lambda), \lambda\sim Gamma(r_0,b_0)$</span> hierarchy.</p>
<p><span class="math-container">$$\begin{align}
P(Y=y)&amp;=\int_{0}^{\infty}P(Y=y|\lambda)f(\lambda)d\lambda\\
&amp;=\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\\
&amp;=\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\\
&amp;=\frac{b_0^{r_0}}{\Gamma(r_0)y!}\int_{0}^{\infty}e^{-\lambda(b_0+1)}\lambda^{y+r_0-1}d\lambda\\
&amp;=\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\\
&amp;=\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\\
&amp;=\frac{\Gamma(r_0+y)}{\Gamma(r_0)y!}\frac{b_0^{r_0}}{(b_0+1)^{r_0+y}}\\
&amp;=\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}$$</span> which is negative binomial distribution.</p>
<p>Negative binomial distribution depicts the number of failures <span class="math-container">$k$</span> before <span class="math-container">$r$</span>th success of a series of independent Bernoulli experiments, with success probability <span class="math-container">$p$</span>:
<span class="math-container">$$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$$</span></p>
<p><span class="math-container">$$\mu=\mathbb E(X)=\frac{r(1-p)}{p}$$</span>
<span class="math-container">$$\sigma^2=\mathbb V(X)=\frac{r(1-p)}{p^2}$$</span></p>
<p>Then
<span class="math-container">$$p=\frac{\mu}{\sigma^2}$$</span>
<span class="math-container">$$r=\frac{\mu^2}{\sigma^2-\mu}$$</span>
The probability mass function (PMF) is:</p>
<p><span class="math-container">$$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}}$$</span></p>
<p>The <span class="math-container">$\mu$</span> and <span class="math-container">$\sigma^2$</span> of each gene can be estimated from data.</p>
<p>Assume that the number of reads in sample <span class="math-container">$j$</span> that are assigned to gene <span class="math-container">$i$</span> can be modeled by a negative binomial (NB) distribution,</p>
<p><span class="math-container">$$K_{ij}\sim NB(\mu_{ij},\sigma^2_{ij})\tag1$$</span> which has two parameters, the mean <span class="math-container">$\mu_{ij}$</span> and the variance <span class="math-container">$\sigma^2_{ij}$</span>. The read counts <span class="math-container">$K_{ij}$</span> are non-negative integers. <span class="math-container">$q_{i, \rho(j)}$</span> is proportional to the expectation value of the true (but unknown) concentration of fragments from gene <span class="math-container">$i$</span> under condition <span class="math-container">$\rho(j)$</span>. The size factor <span class="math-container">$s_j$</span> represents the coverage, or sampling depth, of library <span class="math-container">$j$</span>.</p>
<ul>
<li>First, the mean parameter <span class="math-container">$\mu_{ij}$</span> , that is, the expectation value of the observed counts for gene <span class="math-container">$i$</span> in sample <span class="math-container">$j$</span>, is the product of a condition-dependent per-gene value <span class="math-container">$q_{i, \rho(j)}$</span>(where <span class="math-container">$\rho(j)$</span> is the experimental condition of sample <span class="math-container">$j$</span>) and a size factor <span class="math-container">$s_j$</span>,</li>
</ul>
<p><span class="math-container">$$\mu_{ij} = q_{i, \rho(j)}s_j\tag2$$</span></p>
<ul>
<li>Second, the variance <span class="math-container">$\sigma^2_{ij}$</span> is the sum of a shot noise term and a raw variance term,</li>
</ul>
<p><span class="math-container">$$\sigma^2_{ij}=\mu_{ij}+s_j^2v_{i,\rho(j)}\tag3$$</span></p>
<ul>
<li>Third, we assume that the per-gene raw variance parameter <span class="math-container">$v_{i,\rho(j)}$</span> is a smooth function of <span class="math-container">$q_{i, \rho(j)}$</span>,</li>
</ul>
<p><span class="math-container">$$v_{i,ρ(j)}=v_ρ(q_{i, \rho(j)}) \tag4$$</span></p>
<p>This assumption is needed because the number of replicates is typically too low to get a precise estimate of the variance for gene <span class="math-container">$i$</span> 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 <a href="https://genomebiology.biomedcentral.com/articles/10.1186/gb-2010-11-10-r106" rel="nofollow noreferrer" title="article">article</a>).</p>

CC BY-SA 4.0
2021-09-15T14:20:21.950Z
2021-09-17T13:59:11.043Z
2021-09-17T13:59:11.043Z
Dan Li
<p>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.</p>

CC BY-SA 3.0
2018-02-26T13:40:35.337Z
"given its group association" -- what is the group here? a gene?
user5054
CC BY-SA 3.0
2018-02-26T13:43:03.947Z
@user5054 An experimental group (e.g., WT or treatment).
Devon Ryan
CC BY-SA 3.0
2018-02-26T13:45:32.437Z
Where exactly is this group modeled?
user5054
CC BY-SA 3.0
2018-02-26T13:46:58.267Z
See the context surrounding equation 2 in the methods.
Devon Ryan
CC BY-SA 3.0
2018-02-26T17:44:14.650Z
Then you cannot estimate $\mu_{ij}$ if you don't have any information about the grouping of samples?
user5054
CC BY-SA 3.0
2018-02-26T17:45:37.170Z
@user5054 you can’t do differential expression without groups or covariates.
Devon Ryan
CC BY-SA 3.0
2018-03-28T17:53:02.057Z
You can however estimate the parameters of a negative binomial model by treating all samples as if they come from a single group.
Ian Sudbery
CC BY-SA 3.0
2018-02-26T07:48:14.643Z
2018-02-26T07:48:14.643Z
Devon Ryan
<p>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 <a href="https://genomebiology.biomedcentral.com/articles/10.1186/gb-2010-11-10-r106" rel="nofollow noreferrer" title="article">article</a>). 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.</p>

CC BY-SA 4.0
2019-06-18T11:34:44.200Z
@TomKelly It is, in fact, extremely common. The two most commonly used packages for differential expression use it (edgeR, DESeq[2]).
Konrad Rudolph
CC BY-SA 4.0
2023-05-20T22:51:40.610Z
Is the model wrong ? I thought the NB has as parameters the mean and dispersion (and not the variance).
Picasa
CC BY-SA 3.0
2018-02-26T08:05:22.873Z
The negative binomial distribution is not fitted for a single gene-sample, but along a sample
llrs
CC BY-SA 3.0
2018-02-26T10:36:32.177Z
“how come we can learn the $\mu_{ij}$ from that single observation?” — we can’t. This fundamentally doesn’t work.
Konrad Rudolph
CC BY-SA 3.0
2018-02-26T13:43:05.620Z
@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?
user5054
CC BY-SA 3.0
2018-02-26T14:55:44.017Z
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
CC BY-SA 4.0
2018-06-02T16:18:12.977Z
I don’t think this is common. I’ve seen Poisson distributions suggested for modelling discrete read counts and normal distribution for continuous log-transformed data. Negative binomial would be a relatively uncommon approach. With such a wide dynamic range (even for old arrays) it’s safe to assume the central limit theorem: the most common way to simulate log-expression values is with a multivariate-normal distribution.
Tom Kelly ケリー・トム
<rna-seq><deseq2>
Negative binomial modeling of RNA-Seq data
# Negative binomial modeling of RNA-Seq data

(Asked by: username_0 on 2018-02-26)

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](https://genomebiology.biomedcentral.com/articles/10.1186/gb-2010-11-10-r106 "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.

**Tags:** rna-seq, deseq2
> **Comments:**
> **username_3**: The negative binomial distribution is not fitted for a single gene-sample, but along a sample  
> **username_1**: “how come we can learn the $\mu_{ij}$ from that single observation?” — we can’t. This fundamentally doesn’t work.  
> **username_0**: @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?  
> **username_3**: 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.  
> **username_4**: I don’t think this is common. I’ve seen Poisson distributions suggested for modelling discrete read counts and normal distribution for continuous log-transformed data. Negative binomial would be a relatively uncommon approach. With such a wide dynamic range (even for old arrays) it’s safe to assume the central limit theorem: the most common way to simulate log-expression values is with a multivariate-normal distribution.  
> **username_1**: @TomKelly It is, in fact, extremely common. The two most commonly used packages for differential expression use it (edgeR, DESeq[2]).  
> **username_2**: Is the model wrong ? I thought the NB has as parameters the mean and dispersion (and not the variance).  

---

## 2 answers

### Answer by username_5 on 2021-09-15. Score: 2

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](https://genomebiology.biomedcentral.com/articles/10.1186/gb-2010-11-10-r106 "article")).


---

### Answer by username_6 on 2018-02-26. Score: 1

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.

> **Comments:**
> **username_0**: "given its group association" -- what is the group here? a gene?  
> **username_6**: @user5054 An experimental group (e.g., WT or treatment).  
> **username_0**: Where exactly is this group modeled?  
> **username_6**: See the context surrounding equation 2 in the methods.  
> **username_0**: Then you cannot estimate $\mu_{ij}$ if you don't have any information about the grouping of samples?  
> **username_6**: @user5054 you can’t do differential expression without groups or covariates.  
> **username_7**: You can however estimate the parameters of a negative binomial model by treating all samples as if they come from a single group.  

---

