生物 DeSeq2

，这是 负二项分布（Negative Binomial distribution）在广义线性模型中的一个特性：

$$\text{Var}(Y)=\mu +\alpha {\mu }^2$$

• 其中：

  • $\mu =E[Y]$：基因表达的均值
  • $\alpha$：离散度参数（dispersion）

• 当 $\alpha =0$，就退化为泊松分布（Poisson）

所以，这个方差结构（mean + overdispersion）就是负二项建模的核心原因。

拟合 dispersion-mean 曲线用的是局部回归吗？怎么判断用线性还是多项式？

DESeq2 默认使用的是：

$$\log ({\alpha }_i)=a+\frac{b}{{\mu }_i}$$

这是一个经验拟合公式，不是标准局部回归（loess），但起到类似的“拟合趋势”的作用。

• 具体过程是基于类似非线性最小二乘方法拟合的 parametric 曲线，而不是非参数的 loess 或 spline。
• 也就是说，不是多项式，也不是滑动平均，而是预设的非线性函数形式。

得到平滑线后，只是把极端 dispersion 的基因往平滑线收缩吗？

这就是所谓的 dispersion shrinkage 或 shrink to trend：

• 对每个基因 i，先估计它自己的 raw dispersion ${\hat{\alpha }}_i^{\text{MLE}}$
• 然后根据它对应的 ${\mu }_i$ 在拟合的 trend curve 上找到对应的 ${\hat{\alpha }}_i^{\text{trend}}$
• 最终的 ${\hat{\alpha }}_i^{\text{final}}$ 是这两者的 加权平均（shrinkage）

这一步是为了防止 low-count 基因的 dispersion 被极端高估，增加稳定性。

---

使用所有点进行负二项建模，此时 log2FC 是否直接来自 β1？
$\log ({\mu }_{ij})={\beta }_0+{\beta }_{1}x$，β 是什么？为什么 β1 是 log2FC？

是的，在建模中：

$$\log ({\mu }_{ij})={\beta }_0+{\beta }_1{x}_j$$

• ${\mu }_{ij}$：第 i 个基因在第 j 个样本的期望表达量

• $x_j$：样本的条件（比如 0=control，1=treated）

• ${\beta }_0$：对照组的 log 表达基准

• ${\beta }_1$：处理组相对对照组的 log fold change

• ${\beta }_1={\log }_2(\text{FC})\cdot {\log }_{e}2$，是以 自然对数为底的 log fold change

• results() 输出的 log2FC 是对 ${\beta }_1$ 做换底处理后的结果：

$$\text{log2FC}=\frac{{\beta }_1}{{\log }_{e}2}$$

步骤	内容	数学含义
1. 归一化	计算 size factor	消除测序深度差异
2. dispersion estimation	估计 raw dispersion	${\alpha }_i^{MLE}$
3. dispersion shrinkage	shrink to trend line	${\alpha }_i^{final}$
4. GLM 拟合	使用负二项 GLM 拟合每个基因	$\log ({\mu }_{ij})={\beta }_0+{\beta }_1{x}_j$
5. Wald 检验	对 ${\beta }_1$ 做显著性检验	得到 log2FC 和 p 值

所以使用负二项分布建立了？

DESeq2 用负二项 GLM 模型对 count 数据建模，其方差形式是 $\mu +\alpha {\mu }^2$。通过对 log(μ) 建立线性模型估计 β 系数，其中 β₁ 就是处理组 vs 对照组的 log2 fold change。p 值来自对 β₁ 做 Wald 检验。

GLM（广义线性模型）是一整类模型，其核心形式是：

$$g(\mu )=X\beta$$

• $\mu =E[Y]$：响应变量的期望（可以是 count、概率、连续值都可以）

• $X\beta$：线性预测器（解释变量与系数相乘）

• $g(\cdot )$：连接函数（link function），连接期望值与线性部分

• 响应变量 Y 可以服从不同的分布：

  • 二项分布（logistic 回归）
  • 泊松分布（计数）
  • 负二项分布（overdispersion count）
  • 正态分布（普通线性回归）

所以 DESeq2 中的负二项回归，就是一个特定形式的 GLM：响应变量 Y 服从负二项分布，连接函数是 log()

📌 三、那 dispersion-mean 图是怎么来的？

在 RNA-seq 中，count 数据远远超出了 Poisson 模型假设的均值 = 方差，所以我们使用 负二项分布：

$$\text{Var}(Y)=\mu +\alpha {\mu }^2$$

• $\mu$：期望表达（每个基因平均表达量）
• $\alpha$：基因的 dispersion（生物学变异）

DESeq2 做的是：

1. 对每个基因，计算一个粗略估计的 dispersion。

2. 画出所有基因的 dispersion vs mean（均值）图。

3. 拟合一条平滑趋势线（mean-dispersion 关系）

4. 使用 shrinkage 把每个基因的 raw dispersion 拉向这条趋势线

   • 避免在低表达下离散度估计过于极端

最后每个基因都有一个稳定的 ${\mu }_i$、${\alpha }_i$，用来进行负二项回归建模。

1. 计算标准因子
2. 标准化 count
3. 画出离散度和均值的关系图
4. 拟合负二项模型
5. 对离散度做 shrinkage
6. 得到平滑曲线（dispersion ~ mean）
7. 使用 beta 系数的 t 检验，得到 log2FC 和 p 值

步骤	是否正确	说明
1. 计算标准因子	✅ 对	用 estimateSizeFactors() 计算每个样本的库大小因子
2. 标准化 count	✅ 对	用 size factors 将 raw counts 标准化
3. 画 dispersion vs mean	✅ 对	展示原始 dispersion 和平均表达的关系
4. 拟合 dispersion 模型	✅ 对	使用局部回归拟合 mean-dispersion 关系（用于 shrinkage）
5. shrinkage 收缩离散度	✅ 对	将基因的 raw dispersion 向全局趋势曲线收缩，避免过拟合
6. 得到平滑曲线	✅ 对	没错，最后看到的就是 smooth 的 mean-dispersion 曲线
7. 用 β 系数做 t 检验 得到 log2FC & p 值	✅ 基本正确	不完全是 t 检验，但思路类似（见下）

✅ beta 系数（log2FoldChange）

在 DESeq2 中，我们对每个基因建立一个 广义线性模型（GLM）：

$$\log ({\mu }_{ij})=\log (s_j)+{\beta }_0+{\beta }_1\cdot x_j$$

其中：

• ${\mu }_{ij}$ 是第 i 个基因在第 j 个样本的期望表达量
• $s_j$ 是样本的 size factor（用于标准化）
• $x_j$ 是实验组变量（如 control = 0，treatment = 1）
• ${\beta }_1$ 就是我们关心的 log2 fold change（log2FC）

---

✅ p 值是怎么来的？

DESeq2 对 β 系数进行假设检验：

• 零假设：${\beta }_1=0$（即无差异表达）
• 使用 Wald 检验（Wald test） 或 LRT（似然比检验）
• Wald 检验的统计量：

$$W=\frac{\widehat{{\beta }_1}}{\text{SE}(\widehat{{\beta }_1})}$$

• 这个类似 t 检验，计算 log2FC 相对于其标准误的比值，然后查正态分布得 p 值。

我已经有了表达量和观察到的方差，也通过局部回归得到了平滑的理论关系（方差 vs. 表达均值曲线），那么：
为什么还要做“负二项分布建模”？我不是已经知道 alpha（离散度）了吗？怎么又突然进入线性建模了？

得到 方差 vs. 均值的关系（用于估计 α） 是为了建模 count 的变异性；而 线性模型的目标是估计组间差异（log2FC，β）。
这两个模型一起构成了完整的 负二项 GLM，分别建模：

• 期望值（均值 μ） → 用于检测组间差异（log2FC）
• 方差（α） → 用于估计不确定性（决定统计显著性）

核心：你已经知道了 μ 和 α，那还建模什么？

✅ 第一步：你知道的“理论方差”来自哪？

你说的“理论值”是通过如下方式来的：

• 对每个基因，算出平均表达量（μ_i）和离散度（α_i，或者是原始观察值）
• 再通过 局部回归（LOESS），拟合出表达均值 vs. α 的关系（即平滑线）

这一条“平滑线”叫做 期望离散度曲线（expected dispersion trend）

这只是一个经验性估计，用来缓解每个基因的 raw dispersion 的波动性（特别是在低表达量时容易估计不准）

你最终得到的，是：

• 每个基因的 μ_i（表达均值）
• 每个基因的 α_i（经 shrink 后的 dispersion）

👉 它只告诉你“这个基因的表达大致有多不稳定”，但不能告诉你组间有没有差异！

第二步：负二项 GLM 是建模“表达均值与组变量的关系”

你真正关心的问题是：“这个基因在两个组之间表达是否有差异？”

所以你必须使用：

$$\log ({\mu }_{ij})=\log (s_j)+{\beta }_0+{\beta }_1{x}_j$$

这个 广义线性模型（GLM） 去拟合表达均值（μ）与组别（x）之间的关系。

• 得到的 β₁ 就是你要的 log2 fold change
• 你关心的是：β₁ 是否显著 ≠ 0 → 差异显著性

这个模型是通过最大化如下的 负二项似然函数 来估计参数 β：

$$L(\beta )=\prod _{i,j}P(Y_{ij}\mid {\mu }_{ij}(\beta ),{\alpha }_i)$$

你必须知道 α 才能估计 β 的置信区间（用来计算 p 值），否则你就不知道这个 β 靠谱不靠谱！

内容	用途	建模方式	输出
μ 和 α 的估计	建模表达期望和变异	经验估计 + loess 平滑	每个基因一个 α
负二项 GLM	建模表达与组变量关系	最大似然估计 β	log2FC 和 p 值
两者关系	α 提供了方差结构	β 提供了差异估计	缺一不可

得到了方差和 μ，为什么不能直接得到线性拟合器？

你可以直接拟合，但如果没有 α，你：

• 没法计算误差项（standard error）
• 没法估计 β 是否显著（p 值）
• 没法控制误判率（FDR）

这就像你做线性回归，但不给你 residual variance，你就不知道哪个 β 是显著的、哪个只是 noise。

你不是为了得到 α 而建模；你是为了得到 β（log2FC）和其显著性，而必须通过知道 μ 和 α 来完成完整的负二项 GLM 拟合。

是的，你可以在没有 α 的情况下，拟合出 β —— 也就是所谓的“拟合值”。
但如果你想要p 值（显著性判断），必须要知道 α，因为 α 决定了每个数据点的变异性/误差项。

那为什么还要知道 α 呢？

因为你不只是要 β，你还要知道：

• β 是否显著 ≠ 0（也就是 p 值）
• β 的标准误差是多少（误差越大，置信度越低）

要做这些事情时，你必须知道误差的方差结构，而在 DESeq2 的负二项模型中，误差不是固定的，而是：

$$Var(Y_{ij})={\mu }_{ij}+{\alpha }_i{\mu }_{ij}^2$$

这里的 α 是必须的，否则你无法计算：

$$SE(\hat{\beta })\Rightarrow \text{无法计算 z-score 或 Wald test 的 p 值}$$

你要的是什么	α 是否必要？	为什么
点估计 β（log2FC）	❌ 不需要 α	可通过拟合表达 vs. 组别
显著性检测（p 值）	✅ 需要 α	需要估计误差和变异性
Wald test / LRT	✅ 需要 α	必须知道 Var(Y) 的精确形式

1. 估计 size factors（标准化）
2. 估计 raw dispersion（α raw）
3. 使用 LOESS 平滑得到 α trend
4. shrink α 得到稳定估计
5. 构建 GLM： log(μ) = Xβ
6. 使用 α + μ 共同决定 variance
7. 最大似然估计 β，基于负二项分布
8. Wald 检验或 LRT 得到 p 值

---

🧠 类比理解

想象你测量两组人身高差异（β 是身高差异），你可以估出平均差（点估计），但你不知道每组身高波动有多大（α，方差），就无法判断这个差异是否 显著，因为可能只是噪声。

怎么从负二项分布的 log-likelihood 中用 α 来推导出 β 的误差项（Fisher 信息 / z-score）计算方式。是否需要？
非常好！你现在想要完整理解：

在 DESeq2 中，如何使用负二项分布的似然函数来估计 β 和 α？特别是 α 和 Y_ij 是如何共同出现在似然函数中的？

我们使用 DESeq2 中的 负二项 GLM（Negative Binomial Generalized Linear Model）：

对于每个基因 $i$ 和样本 $j$：

• $Y_{ij}$ 是 count（观测值）
• ${\mu }_{ij}$ 是期望表达量（均值）
• ${\alpha }_i$ 是基因特异的 dispersion（离散度）
• $s_j$ 是样本的 size factor（归一化因子）
• $X_j$ 是组别矩阵，$\beta$ 是我们要估的系数，最核心是 β₁（log2FC）

目标是使用 负二项分布 建立如下模型：

---

📌 1. 负二项分布的定义

DESeq2 使用的负二项分布形式如下：

$$Y_{ij}\sim \text{NB}({\mu }_{ij},{\alpha }_i)$$

这个负二项分布的均值和方差分别为：

• $\mathbb{E}[Y_{ij}]={\mu }_{ij}$
• $\mathrm{Var}(Y_{ij})={\mu }_{ij}+{\alpha }_i{\mu }_{ij}^2$

这里的 α 控制 count 的“额外离散度”，α 越大，方差越高。

---

📌 2. 概率质量函数（PMF）

但在 DESeq2 中，参数改写成 均值 $\mu$ 和 dispersion $\alpha$，更方便解释。

记：

$$r=\frac{1}{{\alpha }_i},p=\frac{r}{r+{\mu }_{ij}}=\frac{1}{1+{\alpha }_i{\mu }_{ij}}$$

代入后可得：

$$P(Y_{ij}=k)=\frac{\Gamma (k+\frac{1}{{\alpha }_i})}{k!\Gamma (\frac{1}{{\alpha }_i})}{\left(\frac{1}{1+{\alpha }_i{\mu }_{ij}}\right)}^{\frac{1}{{\alpha }_i}}{\left(\frac{{\alpha }_i{\mu }_{ij}}{1+{\alpha }_i{\mu }_{ij}}\right)}^k$$

---

📌 3. 负二项对数似然函数（log-likelihood）

我们将上面的概率对所有观测样本 $j$ 取对数求和：

$$\log L_i(\beta ,{\alpha }_i)=\sum _j\left[\log \Gamma (Y_{ij}+\frac{1}{{\alpha }_i})-\log \Gamma (Y_{ij}+1)-\log \Gamma (\frac{1}{{\alpha }_i})+\frac{1}{{\alpha }_i}\log \left(\frac{1}{1+{\alpha }_i{\mu }_{ij}}\right)+Y_{ij}\log \left(\frac{{\alpha }_i{\mu }_{ij}}{1+{\alpha }_i{\mu }_{ij}}\right)\right]$$

记住这个公式非常关键。它是 DESeq2 估计 β 和 α 的核心。

---

📌 4. ${\mu }_{ij}$ 是如何建模的？

$${\mu }_{ij}=s_j\cdot q_{ij},\text{其中}\log q_{ij}=x_j^T\beta$$

也就是说：

$$\log {\mu }_{ij}=\log s_j+x_j^T\beta$$

从而把设计矩阵 $X$ 和待估参数 β 带入模型。

---

📌 5. 求解方式：最大化对数似然

✅ 步骤一：给定 α，求解 β（log2FC）

DESeq2 首先用经验贝叶斯的方法估计好每个基因的 α（使用 LOESS trend + shrinkage），然后 固定 α，最大化上面 log-likelihood，求解 β（log2 fold change）：

$${\hat{\beta }}_i=\arg \underset{\beta }{\max }\log L_i(\beta \mid {\alpha }_i)$$

此时，α 是常数，只影响标准误（下面会讲）。

---

✅ 步骤二：已知 β，估计 α

DESeq2 在前期也会尝试估计每个基因的 α。这时候是：

$${\hat{\alpha }}_i=\arg \underset{\alpha }{\max }\log L_i({\beta }_i,\alpha )$$

不过实际流程中，一般先估 α，再用 α 去估 β（否则会高维难优化）。

---

📌 6. 为什么 α 影响标准误？

我们估计出 β 后，还要知道它的标准误差（standard error），才能做显著性检验。

在负二项 GLM 中，标准误差的计算依赖于观察值的方差：

$$\mathrm{Var}(Y_{ij})={\mu }_{ij}+{\alpha }_i{\mu }_{ij}^2$$

这会被用于构建 Fisher 信息矩阵（或 Hessian）：

$$\text{Fisher Info}(\beta )=X^{T}WX\text{其中 }W\text{ 是根据 Var(Y) 构建的权重矩阵}$$

这意味着，越离散（α 越大）→ β 的误差越大 → 越难得出显著性差异（p 值越大）。

---

✅ 最终：beta 是什么？

在这个模型中：

$$\log {\mu }_{ij}=\log s_j+{\beta }_0+{\beta }_1{x}_j$$

所以：

• ${\beta }_1$ 表示 组别对 log(μ) 的影响
• 也就是 log2FC（log fold change，底数为 e）
• 可转为 log2 fold change：${\log }_{2}FC={\beta }_1/\log (2)$

项目	作用	是否需要 α
β（log2FC）	表达差异	✅ α 决定方差结构
负二项似然函数	联合估计 β 和 α	两者耦合出现
α 的估计	控制变异性，防止假阳性	用于构造正确的统计模型
Fisher 信息（标准误）	用于计算 p 值	必须用到 α 和 μ