batch correction

如果我们只把 Batch Correction 当成一种数据降噪/平滑技术（例如调用 ComBat、limma::removeBatchEffect 或 Harmony），而忽视了它背后的因果图结构（Causal DAG），就极易陷入致命的因果谬误。

批次矫正的因果问题，核心可以总结为：“你矫正掉的，到底纯粹是技术噪音，还是被你无意中连带抹杀（或凭空创造）的生物学信号？”

---

一、 核心因果图（DAG）与三大灾难场景

在因果推断框架中，我们关注三个核心变量：

• $X$ (Treatment / Biology)：你关心的核心生物学变量（如：健康 vs 癌症、野生型 vs 基因敲除）。
• $B$ (Batch / Technical Noise)：技术批次（如：测序板、抽提日期、实验人员、不同医院）。
• $Y$ (Omics Measurement)：测得的组学数据（表达量、蛋白质丰度、代谢物峰面积）。

我们可以通过三种不同的因果结构，来剖析批次矫正的常见陷阱：

---

3. 灾难三：中介效应与生物批次（Mediation & Biological Batches）

• 因果结构：$X\to B\to Y$（$B$ 不是技术噪音，而是生物学中介变量/上下文）。

• 场景：

• 你研究不同医院（Batch）收集的肿瘤样本。

• 但医院 A 是顶尖肿瘤专科医院，接收的都是晚期/重症患者；医院 B 是社区医院，接收的多是早期患者。

• 此时，“医院（Batch）”这个变量不仅包含了技术噪音（测序仪差异），还包含了真正的生物学病情严重程度（Disease Severity）。

• 因果陷阱：

• 如果你粗暴地把“医院”作为 Batch 矫正掉，你不仅去除了仪器噪音，还顺带把“早期 vs 晚期”的真实生物学差异作为中介效应一并抹杀了。

---

二、 为什么非线性/单细胞批次矫正（如 Harmony、MNN、Seurat Integration）更容易产生因果欺骗？

在单细胞 RNA-seq（scRNA-seq）时代，因果问题变得更加复杂：

1. 强行对齐（Over-integration / Over-clustering）：

• 假设 Batch 1（对照组）里有 T细胞 和 B细胞；Batch 2（严重感染组）里由于免疫耗竭，**完全没有 B 细胞，取而代之产生了一群全新的 耗竭 T 细胞**。
• 如果使用基于流形对齐**（Manifold Alignment）**的非线性批次矫正算法（如 MNN 或 Seurat v4），算法的因果假设是“两组数据应该具有相同的细胞类型结构”。
• 后果：算法会强行把 Batch 2 的 耗竭 T 细胞 扭曲并重叠拉扯到 Batch 1 的 B 细胞 上，宣称它们是同一种细胞。这种“强行对齐”直接毁灭了真正的生物学发现（产生假的细胞类型合并）。

---

三、 防范批次矫正因果陷阱的工程与哲学准则

为了防止在批次矫正中被数学算法“欺骗”，资深生信专家遵循以下因果准则：

1. 实验设计优先（Design-level Causal Control）

• 平衡设计（Block Randomization）：在做实验前，确保每一个 Batch 里都按比例均匀混合对照组和处理组（如：Batch 1 包含 5 个 Ctrl + 5 个 Treat，Batch 2 同样包含 5 个 Ctrl + 5 个 Treat）。
• 这是保证 $X\perp B$（$X$ 与 $B$ 相互独立）的唯一因果路径。

2. “拟合时纳入”胜过“显式减去”（Model-in, don’t subtract）

• 不推荐的做法：先用 removeBatchEffect() 抹掉批次，拿出一个“干净”的矩阵，再传给下游做差异分析。
• 因果严谨的做法：保持原始数据矩阵不变，直接在广义线性模型（GLM，如 DESeq2 / limma）的设计矩阵中将 Batch 作为协变量纳入：

$\text{Expression}\sim \text{Batch}+\text{Treatment}$

• 因果含义：在控制了 Batch 变量的前提下，评估 Treatment 的净因果效应（Conditional Independence）。这避免了将矫正后的矩阵导出导致的自由度失真和信息泄露。

3. 区分“可视化矩阵”与“假说检验矩阵”

• PCA/UMAP 展示：可以使用经过 Batch 矫正的矩阵，用于无监督的视觉降噪。
• 差异表达 / 统计推断：永远使用原始未经矫正的 Counts/表达量，并在统计模型内部通过协变量（Covariates）来控制批次。

批次2

在差异表达分析中，真正决定一个生信分析师水平的，不是“会不会运行 DESeq2”，而是能否将复杂的生物学假设精准翻译为线性模型（Design Matrix），并从模型系数（Coefficients）中提炼出无污染的因果效应。

一、 复杂实验场景设定

假设我们正在进行一项药物免疫反应研究，实验设计如下：

• 个体/生物学重复（Individual / Mouse）：共 6 只小鼠（M1 ~ M6），具有基因背景差异。
• 处理组（Treatment）：给药组（Drug）vs 对照组（Vehicle）。
• 时间序列（Time）：分别在给药后 0h（基线） 和 24h 采集同只小鼠的血液样本（配对设计 / Paired Design）。
• 技术批次（Batch）：由于样本处理量限制，M1~M3 的 RNA 抽提与建库由 Batch 1 完成，M4~M6 由 Batch 2 完成。

二、 核心任务 1：设计矩阵（Design Matrix）构建

为了回答“药物在 24h 诱导的特异性转录响应（排除时间本身带来的应激，并控制个体及批次差异）”，我们需要构建广义线性模型（GLM）。

1. 变量属性与冲突排除（因果推断第一步）

在写公式之前，必须先检查变量之间的共线性（Multicollinearity）：

• 小鼠个体（Mouse）与批次（Batch）是什么关系？M1~M3 属于 Batch1，M4~M6 属于 Batch2。这意味着 Batch 变量完全嵌套（Nested）在 Mouse 变量内部！
• 因果定律：如果把 ~ Batch + Mouse + ... 一起写入模型，模型矩阵会发生完全共线性（Rank Deficient），DESeq2 会直接报错崩溃。因为控制了 Mouse，就已经在数学上自动控制了 Batch。

2. 正确的 Design Matrix 写法

我们采用配对设计（Paired Design）+ 交互作用（Interaction）模型：

# 在 colData 中，确保 Mouse、Treatment、Time 均为因子 (Factor)
colData$Mouse     <- factor(colData$Mouse)
colData$Treatment <- factor(colData$Treatment, levels = c("Vehicle", "Drug"))
colData$Time      <- factor(colData$Time,      levels = c("0h", "24h"))
 
# 推荐写法：嵌套模型 (Nested / Grouped Design)
# 公式：~ Mouse + Treatment + Time + Treatment:Time
design(dds) <- ~ Mouse + Treatment * Time
 

---

三、 核心任务 2：模型系数（Coefficients）的生物学含义解构

在运行 dds <- DESeq(dds) 后，线性模型会为每一个基因拟合如下的对数表达量方程：

${\log }_2(\text{Expression})={\beta }_0+{\beta }_{\text{Mouse}}\cdot X_{\text{Mouse}}+{\beta }_{\text{Drug}}\cdot X_{\text{Drug}}+{\beta }_{\text{24h}}\cdot X_{\text{24h}}+{\beta }_{\text{Drug:24h}}\cdot (X_{\text{Drug}}\cdot X_{\text{24h}})$

每个系数（Coefficient）代表的精准生物学含义如下：

系数 (Coefficient)	命名 (DESeq2 results)	数学表达	真实的生物学含义
${\beta }_0$	Intercept	基线表达量	Vehicle 组在 0h 的参考基线表达量（结合 Mouse1 的个体基线）。
${\beta }_{\text{Mouse}}$	Mouse_M2_vs_M1 等	${\text{Base}}_{\text{Mouse}}-{\text{Base}}_{\text{M1}}$	个体/批次基线差异：吸收小鼠间的固化固有表达差异及批次效应。
${\beta }_{\text{Drug}}$	Treatment_Drug_vs_Vehicle	${\text{Drug}}_{\text{0h}}-{\text{Vehicle}}_{\text{0h}}$	给药前（0h）的初始组间偏置：理论上应接近 0。若偏离 0 很大，说明随机分组合班时基线不平。
${\beta }_{\text{24h}}$	Time_24h_vs_0h	${\text{Vehicle}}_{\text{24h}}-{\text{Vehicle}}_{\text{0h}}$	纯粹的时间/生理节律/对照组应激效应：小鼠不给药只给溶剂，随着时间推移 24h 自然发生的基因表达变化。
${\beta }_{\text{Drug:24h}}$	TreatmentDrug.Time24h	$({\text{Drug}}_{\text{24h}}-{\text{Drug}}_{\text{0h}})-({\text{Vehicle}}_{\text{24h}}-{\text{Vehicle}}_{\text{0h}})$	真正的药物特异性响应（Interaction Effect / DID）：减去了时间应激后，药物在 24h 诱导的净响应。

如何通过 results() 提取不同生物学问题的差异基因？

1. 问题 A：我想找“真正由药物诱导的差异基因”（排除了 24h 正常生理波动/应激）：

# 直接提取交互项系数
res_drug_effect <- results(dds, name = "TreatmentDrug.Time24h")
 

• 因果含义：这是双重差分法（Difference-in-Differences, DID）。它算的是：(Drug24h - Drug0h) - (Vehicle24h - Vehicle0h)。

2. 问题 B：我想找“药物组在 24h 相较于自己 0h”改变了哪些基因（包含药物效应 + 24h 基础生理波动）：

# 联合系数：\beta_{24h} + \beta_{Drug:24h}
res_drug_total_24h <- results(dds, contrast = list(c("Time_24h_vs_0h", "TreatmentDrug.Time24h")))
 

---

四、 核心任务 3：预测可能的混杂因素与实验设计漏洞

• 危机场景：Drug 组毒性过大，在 24h 采样前，6 只小鼠死了 3 只。

1. 后果：极大量的昼夜节律基因（如 PER1, CLOCK）会发生巨大的表达量翻转，其方差会轻松压过药物响应的方差。 **做法：**统一时间点
2. 采样批次： 但是 Monday 和 Tuesday 离心机的温度设定不一致，或者 RNA 提取试剂盒换了一盒。
3. 存活性
   • 后果：能活到 24h 被采血的 3 只小鼠，可能本身就具有某种特殊的抗性基因型。这种选择性失踪（Collider Bias）会严重扭曲统计推断。
   • 防范方案：
   • 在 Design Matrix 中记录死亡/濒死状态作为协变量。
   • 实验设计时需考虑药效学 LD50，避免过高致死剂量；如果发生失踪，不能简单把缺失样本删除，需进行敏感性分析。

---

总结：高级实验设计的思维模型

[原始组学变异 (Total Variance)]
   │
   ├──> 可控生物变异 (Treatment, Time, Treatment:Time) ---> 写入设计矩阵核心项
   │
   ├──> 结构化混杂 (Mouse ID, Batch, Technical Runs)  ---> 作为协变量 (Blocking Variables) 纳入模型
   │
   └──> 隐性/随机噪音 (Unobserved Confounders, Stress) ---> 靠 随机化 (Randomization) + 盲法 抵消