4. 差异分析与特征选择

在多组学药物重定位流程中，差异分析承担着“把全基因组/全蛋白/全代谢物测量压缩为候选特征”的关键任务。它的目标不是简单找出 p 值最小的分子，而是识别出在药物处理与对照之间具有可重复、方向明确、效应足够大的分子特征，为后续重定位打分提供可靠输入。本节围绕转录组差异表达（limma、DESeq2）、蛋白与代谢物筛选、多重检验校正、效应量评估以及特征稳定性检验展开。

一、转录组差异表达：limma 与 DESeq2 的适用边界

limma 与 DESeq2 是两类最常用的差异表达工具，但它们的统计假设不同，不能混用。limma 基于线性模型与经验贝叶斯方差收缩，最初面向 log2 荧光强度等近似正态的连续数据；当用于 RNA-seq counts 时，需先做 voom 变换，把 counts 转换为带权重的 log2 CPM，再进行线性建模。DESeq2 则直接对原始 counts 建模，假设其服从负二项分布，通过中位数比值法估计 size factor，并用离散度收缩稳定方差。

选择依据可以概括为：样本量小、counts 数据、需要严格处理离散度时优先 DESeq2；样本量较大、已有 log 尺度表达矩阵、或需要灵活设计矩阵（如多因素、交互项、批次校正）时 limma 更高效。下面给出一个可直接参考的 DESeq2 分析片段（示例参数，非实测结果）：

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = coldata,
                              design = ~ batch + condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "drug", "control"),
               alpha = 0.05, lfcThreshold = 0.5)
res <- lfcShrink(dds, coef = "condition_drug_vs_control", type = "apeglm")
summary(res)

其中 lfcThreshold = 0.5 表示只检验 |log2FC| > 0.5 的假设，lfcShrink 用于收缩低表达基因的效应量，避免“低 counts + 大 fold change”的假阳性。limma-voom 的对应写法为：

library(limma); library(edgeR)
dge <- DGEList(counts = counts)
dge <- calcNormFactors(dge)
design <- model.matrix(~ batch + condition, data = coldata)
v <- voom(dge, design)
fit <- lmFit(v, design)
fit <- eBayes(fit)
topTable(fit, coef = "conditiondrug", adjust.method = "BH", n = 20)

二、差异蛋白与代谢物筛选的特殊性

蛋白组与代谢组的数据结构与转录组差异明显。蛋白组常来自 DIA/LFQ 强度，缺失值比例高，且存在“缺失非随机”问题；代谢组则峰数量大、注释层级不一、浓度动态范围宽。因此不能直接套用 RNA-seq 的 counts 模型。

蛋白差异筛选通常先做 log2 转换与缺失值过滤（如保留至少 70% 样本中检出的蛋白），再用 limma 对 log2 强度建模，并加入 arrayWeights 或 duplicateCorrelation 处理样本质量差异。代谢物筛选则更依赖标准化：先做 PQN（概率商归一化）或内标校正，再对 log 转换后的强度做 limma 或非参数检验。对于代谢物，效应量常以 log2FC 与 VIP（OPLS-DA 变量投影重要性）联合判断，VIP > 1 且 FDR < 0.05 是常见组合阈值。

三、多重检验校正：从 p 值到 q 值

一次差异分析会同时检验上万个分子，若以 p < 0.05 为阈值，假阳性数量约为 0.05 × 检验数。例如 20000 个基因会产生约 1000 个假阳性。因此必须做多重检验校正。常用方法有三类：

• Bonferroni：p × m，控制族错误率 FWER，最严格，适合检验数少或必须零假阳性的场景。
• Benjamini-Hochberg（BH）：控制 FDR，最常用，q 值 = p × m / rank，适合高通量筛选。
• Storey q-value：在 BH 基础上估计真实零假设比例 π0，检验数大且信号丰富时更灵敏。

实际操作中，建议同时报告 p 值、BH 校正 q 值和效应量，而不是只给一个“差异基因列表”。判断阈值可设为 q < 0.05 且 |log2FC| > 1；若信号较弱，可放宽到 q < 0.1 并提高效应量门槛，但需在报告中明确说明。

四、效应量评估与特征稳定性检验

p 值回答“是否有差异”，效应量回答“差异有多大、是否值得关注”。常用效应量包括 log2FC、Cohen's d、Hedges' g 以及蛋白/代谢物的 VIP。对于药物重定位，效应量方向比显著性更重要：一个 q = 0.01 但 log2FC = 0.2 的基因，其生物学意义远不如 q = 0.03 且 log2FC = 2 的基因。

特征稳定性检验用于回答“这个差异特征是否可重复”。常用手段有三种：

1. Bootstrap 重采样：对样本有放回抽样 1000 次，统计每个特征被选为差异的频次，保留频次 > 70% 的特征。
2. 交叉验证：将样本分为训练/验证集，在训练集选特征，在验证集检查方向与显著性是否一致。
3. 置换检验：打乱分组标签 1000 次，构建零分布，计算真实差异数的经验 FDR。

下面给出一个 Bootstrap 稳定性检验的示例代码框架：

set.seed(42)
stab <- replicate(1000, {
  idx <- sample(seq_len(ncol(counts)), replace = TRUE)
  dds_b <- DESeqDataSetFromMatrix(counts[, idx], coldata[idx, ], ~ condition)
  dds_b <- DESeq(dds_b)
  r <- results(dds_b, alpha = 0.05)
  rownames(r)[which(r$padj < 0.05 & abs(r$log2FoldChange) > 1)]
})
freq <- table(unlist(stab)) / 1000
stable_features <- names(freq)[freq >= 0.7]

最终进入重定位打分的特征，应是“统计显著 + 效应量足够 + 稳定性达标”三者的交集。建议将差异分析结果整理为统一表格，字段包括：feature_id、log2FC、pvalue、padj、effect_size、stability_freq、direction、omics_layer。这样后续无论是做打分矩阵还是做多组学整合，都能直接复用。