06 验证与对比：怎么知道结果不是 artifact

跑完主线分析拿到一堆显著结果，先别高兴太早。组学分析中假阳性的来源很多：批次效应、样本污染、参数选择、多重检验不充分……这一章讲的是如何用系统性的验证策略确认你的核心发现是可靠的。一个经过充分验证的结果，审稿人很难拒绝。

为什么需要验证

组学数据的维度极高（20000+ 基因），样本量相对有限（TCGA-LIHC 371 例 tumor）。在这种"p >> n"的场景下，过拟合几乎是必然的——如果你用同一份数据既选特征又评估性能，得到的结果一定比真实性能好。这就是为什么审稿人几乎必问"有没有验证"。

验证的核心逻辑是：用独立的数据来检验你在训练数据上得到的结论。"独立"可以有不同层次：同一数据集的不同子集（内部验证）、不同队列的数据（外部验证）、不同平台的数据（跨平台验证）。层次越高，验证的说服力越强。

对于 TCGA-LIHC 项目，我们的验证策略是：先用内部验证确认结果的稳健性，再用 ICGC-LIRI（日本肝癌队列，~240 例）做外部验证，最后用多工具对比确认核心差异基因不依赖于特定工具的选择。

内部验证：Bootstrap 与交叉验证

内部验证不能证明结果可推广，但能证明结果不是由少数极端样本驱动的。两种最常用的方法是 bootstrap 重采样和 K-fold 交叉验证。

Bootstrap 验证 LASSO signature 的稳定性：

# ============================================================
# 内部验证 1: Bootstrap 稳定性检验
# ============================================================
library(glmnet)
library(survival)

# 对 LASSO 做 100 次 bootstrap，看哪些基因稳定入选
set.seed(123)
n_boot <- 100
gene_freq <- rep(0, ncol(x))
names(gene_freq) <- colnames(x)

for (i in 1:n_boot) {
  # 有放回抽样
  boot_idx <- sample(1:nrow(x), replace = TRUE)
  x_boot <- x[boot_idx, ]
  y_boot <- y[boot_idx, ]
  
  # 重新跑 LASSO
  cv_boot <- cv.glmnet(x_boot, y_boot, family = "cox", alpha = 1, nfolds = 10)
  coef_boot <- coef(cv_boot, s = "lambda.min")
  selected <- rownames(coef_boot)[coef_boot[, 1] != 0]
  
  gene_freq[selected] <- gene_freq[selected] + 1
}

# 入选频率 > 70% 的基因被认为是稳定的
stable_genes <- names(gene_freq[gene_freq > 70])
cat(sprintf("Bootstrap 稳定基因 (>70%%): %d / %d\n", 
            length(stable_genes), length(active_genes)))

# 可视化基因入选频率
barplot(sort(gene_freq[gene_freq > 20], decreasing = TRUE),
        las = 2, cex.names = 0.7,
        main = "Gene Selection Frequency in 100 Bootstrap Runs",
        ylab = "Frequency")
abline(h = 70, col = "red", lty = 2)

如果你的 LASSO signature 有 15 个基因，但 bootstrap 显示只有 8 个基因的入选频率超过 70%，那么这 8 个基因是核心基因，其余 7 个可能是噪声。你可以考虑只用这 8 个基因重新构建 signature。

K-fold 交叉验证评估预测性能：

# ============================================================
# 内部验证 2: 5-fold 交叉验证
# ============================================================
library(survcomp)  # for concordance.index

set.seed(42)
n_folds <- 5
fold_ids <- sample(rep(1:n_folds, length.out = nrow(x)))
cv_cindex <- numeric(n_folds)

for (k in 1:n_folds) {
  train_idx <- fold_ids != k
  test_idx <- fold_ids == k
  
  # 在训练集上拟合 LASSO
  cv_train <- cv.glmnet(x[train_idx, ], y[train_idx, ], 
                         family = "cox", alpha = 1, nfolds = 10)
  
  # 在测试集上计算 risk score
  pred_test <- predict(cv_train, newx = x[test_idx, ], s = "lambda.min")
  
  # 计算 C-index
  ci <- concordance.index(
    x = pred_test,
    surv.time = y[test_idx, 1],
    surv.event = y[test_idx, 2],
    method = "noether"
  )
  cv_cindex[k] <- ci$c.index
}

cat(sprintf("5-fold CV C-index: %.3f ± %.3f\n", 
            mean(cv_cindex), sd(cv_cindex)))

C-index 的解读：0.5 表示随机预测，0.7-0.8 表示中等预测能力，>0.8 表示优秀。对于组学 signature 来说，CV C-index 在 0.65-0.75 之间是比较正常的。如果 CV C-index 和 training C-index 差距很大（比如 training 0.85 但 CV 只有 0.60），说明严重过拟合。

外部验证：ICGC-LIRI 数据集

外部验证是最有说服力的验证方式。ICGC-LIRI（Liver Cancer - RIKEN, Japan）包含约 240 例日本肝癌患者的 RNA-seq 数据和临床信息，是验证 TCGA-LIHC 发现的理想数据集。

# ============================================================
# 外部验证: ICGC-LIRI
# ============================================================

# 读取 ICGC-LIRI 表达数据（已归一化为 FPKM）
icgc_expr <- read.csv("ICGC-LIRI_expression.csv", row.names = 1)
icgc_clin <- read.csv("ICGC-LIRI_clinical.csv", row.names = 1)

# 注意：ICGC 用的是 gene symbol，TCGA 用的是 Ensembl ID
# 需要做 ID 转换，确保基因能对应上
library(org.Hs.eg.db)
gene_symbol_map <- bitr(active_genes, 
                        fromType = "ENSEMBL", 
                        toType = "SYMBOL",
                        OrgDb = org.Hs.eg.db)

# 找到在 ICGC 数据中也存在的 signature 基因
common_genes <- intersect(gene_symbol_map$SYMBOL, rownames(icgc_expr))
cat(sprintf("Signature 基因在 ICGC 中匹配: %d / %d\n",
            length(common_genes), length(active_genes)))

# 用 TCGA 训练的系数计算 ICGC 的 risk score
# 注意：需要对 ICGC 数据做和 TCGA 相同的预处理（log2 转换）
icgc_log2 <- log2(icgc_expr[common_genes, ] + 1)

# 提取对应基因的 LASSO 系数
coef_for_validation <- coef_min[gene_symbol_map$ENSEMBL[
  match(common_genes, gene_symbol_map$SYMBOL)], ]

# 计算 risk score
icgc_risk <- as.numeric(t(icgc_log2) %*% coef_for_validation)

# KM 验证
icgc_group <- ifelse(icgc_risk > median(icgc_risk), "High-risk", "Low-risk")

icgc_surv <- data.frame(
  time = icgc_clin$OS.time,
  status = icgc_clin$OS,
  risk_group = factor(icgc_group, levels = c("Low-risk", "High-risk"))
)

fit_icgc <- survfit(Surv(time, status) ~ risk_group, data = icgc_surv)
ggsurvplot(fit_icgc, data = icgc_surv,
           pval = TRUE, risk.table = TRUE,
           title = "External Validation: ICGC-LIRI",
           palette = c("#2E9FDF", "#E7B800"))

外部验证的几个注意事项：

第一，不要期望完全一致的效果。TCGA-LIHC 以美国患者为主，ICGC-LIRI 是日本患者，两个队列的遗传背景、环境因素、治疗方案都不同。如果 TCGA 中 p = 1e-5 而 ICGC 中 p = 0.03，这已经是很好的验证结果了。关键是方向一致（高风险组预后更差）且统计显著。

第二，处理平台差异。TCGA 用的是 Illumina HiSeq，ICGC 可能用不同的测序平台。归一化方法也可能不同（TCGA 提供 raw counts，ICGC 可能只提供 FPKM）。在计算 risk score 时，确保两个数据集经过相同的预处理（都做 log2 转换）。

第三，报告匹配率。如果你的 signature 有 15 个基因但在 ICGC 中只匹配到 12 个，要在文章中说明这一点，并讨论缺失基因对结果的潜在影响。

多工具对比

同一份数据用不同工具分析，结果一致的部分可信度更高。对于差异分析，最常见的对比是 DESeq2 vs edgeR vs limma-voom。

# ============================================================
# 多工具对比: DESeq2 vs edgeR vs limma-voom
# ============================================================
library(edgeR)
library(limma)

# --- edgeR ---
dge <- DGEList(counts = counts_filtered, group = col_data$condition)
dge <- calcNormFactors(dge)
design <- model.matrix(~ condition, data = col_data)
dge <- estimateDisp(dge, design)
fit_edger <- glmQLFit(dge, design)
res_edger <- glmQLFTest(fit_edger, coef = 2)
deg_edger <- topTags(res_edger, n = Inf)$table %>%
  filter(abs(logFC) > 1, FDR < 0.05) %>%
  rownames()

# --- limma-voom ---
v <- voom(dge, design, plot = FALSE)
fit_limma <- lmFit(v, design)
fit_limma <- eBayes(fit_limma)
res_limma <- topTable(fit_limma, coef = 2, number = Inf)
deg_limma <- res_limma %>%
  filter(abs(logFC) > 1, adj.P.Val < 0.05) %>%
  rownames()

# --- 对比 ---
library(VennDiagram)
venn_list <- list(
  DESeq2 = deg$gene_id,
  edgeR = deg_edger,
  limma = deg_limma
)

# 三工具交集
consensus_deg <- Reduce(intersect, venn_list)
cat(sprintf("三工具共同差异基因: %d\n", length(consensus_deg)))
cat(sprintf("DESeq2: %d, edgeR: %d, limma: %d\n",
            length(venn_list$DESeq2), length(venn_list$edgeR), length(venn_list$limma)))

# Venn 图
venn.diagram(venn_list, filename = "DEG_three_tools_venn.png",
             fill = c("#E41A1C", "#377EB8", "#4DAF4A"),
             alpha = 0.5, cat.cex = 1.2)

多工具对比的解读原则：

情况	解读	处理方式
三工具都显著	高可信度	作为核心结果报告
两个工具显著	中等可信度	可以报告，注明一致性
只有一个工具显著	低可信度	谨慎对待，可能是工具特异性 artifact
三工具方向不一致	可能有问题	检查数据质量和预处理步骤

需要注意的是，三个工具的统计模型有相似之处（都基于负二项分布或类似假设），所以"三个都显著"不等于"绝对正确"——它们可能犯同样的错误。但如果连这三个主流工具都不一致，那结果的可信度确实要打折扣。

统计学陷阱与 Artifact 识别

以下是组学分析中最常见的统计学陷阱，每个都配有 TCGA-LIHC 中的具体例子：

陷阱一：Circular Analysis（循环分析）。 用同一份数据既选特征又验证。例如：在全部 371 个 tumor 样本上做 LASSO 选基因，然后在同样的 371 个样本上画 KM 曲线——这不是验证，这是 resubstitution。正确做法：至少做 train/test split 或交叉验证。

陷阱二：Data Leakage（数据泄露）。 在特征选择之前就用了测试集的信息。例如：先对全部样本做 variance filtering（去掉低变异基因），然后再 split 成 train/test——filtering 步骤已经"看到"了 test 数据。正确做法：先 split，再在 train set 上做 filtering。

陷阱三：Confounding（混杂因素）。 TCGA-LIHC 的 normal 样本全部来自 stage I-II 患者（因为晚期患者很少做手术取癌旁组织）。如果你发现某个基因在 tumor 中高表达，它可能不是"tumor vs normal"的差异，而是"late stage vs early stage"的差异。解决方法：在 design formula 中纳入 stage 作为协变量，或者只用 matched pair 做分析。

陷阱四：Overfitting（过拟合）。 LASSO 选了 30 个基因构建 signature，在训练集上 C-index = 0.85，但在外部验证集上只有 0.58。这说明模型记住了训练数据的噪声。解决方法：用更严格的 lambda（lambda.1se 而不是 lambda.min），或者减少候选基因数量。

陷阱五：Multiple Testing 不充分。 对 20000 个基因做单因素 Cox，用 p < 0.05 筛选得到 1500 个"预后相关基因"。但如果所有基因都和预后无关，按 5% 的假阳性率也会有 1000 个假阳性。解决方法：对单因素 Cox 的 p 值也做 BH 校正，用 FDR < 0.05 而不是 p < 0.05。

常见坑

• 只做内部验证不做外部验证：内部验证（CV、bootstrap）只能说明结果在当前数据上是稳定的，不能说明它在新数据上也成立。审稿人几乎必问外部验证。如果实在找不到合适的外部数据集，至少要在文章中把这作为 limitation 讨论。

• 验证数据集选择不当：用一个只有 20 例的 GEO 数据集做验证，样本量太小，即使真实效应存在也可能检测不到。或者用一个不同癌种的数据集做验证（比如用肺癌数据验证肝癌 signature），这在生物学上说不通。

• 把"三个工具都显著"当成金标准：DESeq2、edgeR、limma-voom 的统计假设类似，它们可能对同一种 artifact（如 library size 极端值）都不敏感。真正的"独立验证"应该来自不同的数据，而不是不同的工具。

• 验证失败就放弃整个项目：外部验证中 15 个 signature 基因有 12 个方向一致、3 个不一致，这不是"验证失败"。关键是核心发现（高风险组预后更差）是否在外部数据中重现。个别基因的不一致可以在 Discussion 中讨论。

• 不报告验证失败的结果：如果你试了 3 个外部数据集，2 个验证通过、1 个没通过，应该如实报告所有结果。只报告成功的验证是一种 reporting bias。

下一步

接着深入： 验证通过之后，下一步是把统计结果翻译成生物学语言。07 整合解读：把结果翻译成生物学语言 会讲如何从显著基因和通路中提炼出有意义的生物学故事。

横向延伸： 如果你想深入理解"为什么结果可能是错的"，可以看 01 数据可信度：怎么不被自己骗，那里有 5 个"看上去对但其实错"的真实场景。

参考资源

• Harrell FE Jr (2015) Regression Modeling Strategies. Springer. 关于 bootstrap 验证和 C-index 的经典教材。
• Steyerberg EW et al. (2010) "Assessing the performance of prediction models: a framework for some traditional and novel measures." Epidemiology 21:128-138. 预测模型验证的框架。
• Castaldi PJ et al. (2011) "An empirical assessment of validation practices for molecular classifiers." Briefings in Bioinformatics 12:189-202. 分子分类器验证的实证研究。
• ICGC/TCGA Pan-Cancer Analysis of Whole Genomes Consortium (2020) "Pan-cancer analysis of whole genomes." Nature 578:82-93. ICGC 数据的来源和使用方法。
• Ritchie ME et al. (2015) "limma powers differential expression analyses for RNA-sequencing and microarray studies." Nucleic Acids Research 43:e47. limma-voom 方法论文。