text: user1331: 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.
  
  user77: 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.
  
  user9698: 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).
  
original_text: 1331: <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>
   77: <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>
   9698: <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>
  