01 多组学整合实践教程

单一组学常常回答不了完整问题。转录组看到 mRNA 变化，却不知道蛋白层面有没有跟着改变；表观组解释了染色质开放，却不直接告诉你哪些基因表达了；基因组变异落在了某个启动子上，也需要表达和修饰数据才能说明它的效应。多组学整合的目标不是把所有数据塞进一个矩阵，而是让每一层提供一段证据，共同支撑一个生物学判断。

本专栏讲如何把多层数据放进同一个分析框架：从样本匹配、特征对齐，到常用整合方法（相关、网络、因子模型、机器学习）的适用场景和取舍。

一个项目大致是什么样子

假设你手上有一个队列：50 个样本，每个样本都有转录组、甲基化和蛋白质组数据，外加临床表型（亚型、分期、生存时间）。从这里到"一个能解释某个亚型的多组学特征"，大致要走：

在"整合建模"这一步，常见做法分三类：

• 早期整合（early integration）：把各层特征拼成一个大矩阵，后续按单组学思路处理。简单，但容易被维度高的那一层主导。
• 中期整合（intermediate integration）：在模型层面同时考虑多组学，典型代表是 MOFA 和 mixOmics。
• 晚期整合（late integration）：各层独立分析，最后在结果层面汇总，比如合并 p 值或做共识聚类。

没有谁最好。样本少、特征差异大时 MOFA 更稳；样本充足、特征维度可比时早期整合也能跑通；目标是临床预测时机器学习整合（Random Forest、神经网络）更常见。

常见工具栈

MOFA2 是个值得单独提的工具。它用概率矩阵分解识别"跨组学共同变异"的潜因子，对缺失值友好，对样本量小也能用，现在几乎是多组学整合教学的标配。

推荐公开数据集

多组学公开数据集数量不多，但下面几个足够做教学练习：

MOFA2 的 CLL 数据集特别适合教学：4 个组学层、200 个样本，包里 load 进来就能跑。

最小可跑的例子

下面用 MOFA2 自带的 CLL 数据做一次最简单的整合分析，从加载数据到看出潜因子的生物学含义，不到 15 行代码：

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("MOFA2")

library(MOFA2)

# 载入示例数据：慢淋（CLL）多组学
file <- system.file("extdata", "CLL_data.hdf5", package = "MOFA2")
model <- load_model(file)

# 先看整体：每个组学在模型里的方差贡献
plot_variance_explained(model, x = "view", y = "factor")

# 看某个因子在各层数据上解释了多少方差
plot_variance_explained(model, max_r2 = 10)

# 某个因子上权重最高的特征（基因/甲基化位点等）
plot_top_weights(model, view = "mRNA", factor = 1, nfeatures = 10)

# 把因子得分投射到样本元信息上（比如亚型）
metadata <- samples_metadata(model)
head(metadata)

plot_variance_explained 会输出一张热图，每行是潜因子、每列是组学层。有些因子只在一层里有高贡献（组学特异），有些跨多层都有贡献（跨组学共同变异），后者通常更有生物学意义。

理解了这个例子，再去读 MOFA 论文或 CPTAC 的整合分析会容易很多。

专栏模块规划

所有 8 个模块已上线。Module 05 带可跑脚本和真实数据图。

推荐前置知识

• bulk RNA-seq 实践教程
• 表观组学实践教程
• 蛋白质组学实践教程
• R 数据整理与 ggplot2 可视化

多组学分析的每一层预处理都涉及该组学的专属知识，建议先熟悉至少两个单组学的流程再来做整合。

常见坑

坑 1：期望整合一定比单组学好

很多人认为"数据越多结论越好"，但如果某一层数据质量差或与研究问题无关，强行整合反而引入噪声。跑完整合模型后一定要和单组学结果做对照实验，看是否真的有提升。没有提升本身也是有价值的结论。

坑 2：把相关当因果

跨层相关性分析发现甲基化和表达负相关，不代表甲基化"导致"了沉默——两者可能同时受上游信号（如转录因子结合）驱动。需要因果推断时应考虑 Mendelian randomization 或 mediation analysis，而不是仅凭 Pearson r 下结论。

坑 3：忽略各层数据量级差异导致主导效应

早期整合（拼接矩阵）时，如果 RNA-seq 有 20000 个基因而蛋白组只有 3000 个，模型会被转录组主导。解决方案是在拼接前做特征数平衡（每层取 top N）或用中期整合方法（MOFA2、mixOmics）让各层权重自动学习。

坑 4：样本量不够就上复杂模型

MOFA2 建议至少 15 个配对样本，SNF 在 30 以下不稳定。10 个样本就想跑 multi-omics factor analysis，结果的因子不可复现。样本少时先做简单的 cis 相关和单基因多层验证，不要强上降维模型。

下一步

接着深入：

• 01 多组学项目设计与样本匹配 — 从实验设计角度理解整合的起点
• 05 MOFA2 因子分析实战 — 动手跑最常用的中期整合方法

横向延伸：

• 蛋白质组学实践教程 — 整合中最常搭配的第二层组学
• 表观组学实践教程 — 甲基化数据的预处理和分析基础

参考资源

• MOFA2 官方教程
• mixOmics 文档
• SNFtool 说明
• MultiAssayExperiment Bioconductor
• TCGA GDC Portal
• CPTAC

02 多组学项目设计与样本匹配

速答

Q: 多组学整合第一步做什么？ A: 不是跑模型，而是确认"哪些样本在所有组学层都有数据"。样本 ID 不一致、缺失模式不清楚，后面分析全是空中楼阁。

Q: 不同平台的样本 ID 怎么统一？ A: TCGA 用 barcode（如 TCGA-A1-A0SK-01A-11R-A084-07），蛋白组可能只保留前 12 位，甲基化用 Sentrix ID。需要建一张映射表把各平台 ID 对齐到同一个样本。

Q: 样本缺失模式怎么看？ A: 用 UpSet plot（ComplexUpset 包），比韦恩图更适合 >3 个组学层的情况。能看出哪些样本在所有层都有数据、哪些只在一两层有。

多组学整合的第一步不是跑模型，而是确认"哪些样本在所有组学层都有数据"。样本 ID 不一致、批次不对齐、缺失模式不清楚——这些问题不解决，后面的分析全是空中楼阁。

样本 ID 的现实

不同组学平台对同一个样本的命名方式几乎不会一样。TCGA 用 barcode（TCGA-A1-A0SK-01A-11R-A084-07），蛋白组可能只保留前 12 位，甲基化数据用 Sentrix ID。你需要一张映射表把它们统一起来。

library(dplyr)

# 假设三张表各有自己的 ID 列
rna_samples   <- colnames(rna_mat)
meth_samples  <- colnames(meth_mat)
prot_samples  <- colnames(prot_mat)

# 用映射表统一
id_map <- read.csv("sample_id_mapping.csv")
# 列：sample_id, rna_id, meth_id, prot_id

# 找到三层都有的样本
common <- id_map %>%
  filter(rna_id %in% rna_samples,
         meth_id %in% meth_samples,
         prot_id %in% prot_samples)
cat("三层共有样本数:", nrow(common), "\n")

缺失模式可视化

在决定用哪些样本之前，先画一张缺失模式图。UpSet plot 比 Venn 图更适合三层以上的情况：

library(UpSetR)

sample_list <- list(

  RNA       = id_map$sample_id[!is.na(id_map$rna_id)],
  Methylation = id_map$sample_id[!is.na(id_map$meth_id)],
  Protein   = id_map$sample_id[!is.na(id_map$prot_id)]
)
upset(fromList(sample_list), order.by = "freq")

如果某一层缺失比例很高，需要考虑是否把它排除在整合之外，或者用支持缺失值的方法（如 MOFA2）。

MultiAssayExperiment 容器

Bioconductor 的 MultiAssayExperiment（MAE）是多组学数据的标准容器。它把多层数据、样本映射和临床信息绑在一起，后续分析函数可以直接操作。

library(MultiAssayExperiment)

# 构建 ExperimentList
exp_list <- ExperimentList(
  RNA   = SummarizedExperiment(assays = list(counts = rna_mat)),
  Meth  = SummarizedExperiment(assays = list(beta = meth_mat)),
  Prot  = SummarizedExperiment(assays = list(abundance = prot_mat))
)

# 样本映射表：每行说明某个 primary sample 在某层对应哪个 colname
smap <- listToMap(list(
  RNA  = data.frame(primary = common$sample_id, colname = common$rna_id),
  Meth = data.frame(primary = common$sample_id, colname = common$meth_id),
  Prot = data.frame(primary = common$sample_id, colname = common$prot_id)
))

# 临床信息
col_data <- DataFrame(row.names = common$sample_id,
                      subtype = common$subtype,
                      stage   = common$stage)

mae <- MultiAssayExperiment(experiments = exp_list,
                            colData     = col_data,
                            sampleMap   = smap)
mae

设计阶段的几个决策

1. 样本量：多组学整合对样本量要求不低。MOFA2 官方建议至少 15 个样本；SNF 在 30 个以下效果不稳定。如果某一层只有 10 个样本，考虑是否值得纳入。
2. 批次效应：不同组学的批次效应需要分别处理。RNA-seq 用 ComBat-seq，甲基化用 ComBat，蛋白组用 limma removeBatchEffect。不要在合并矩阵之后再做批次校正。
3. 配对 vs 非配对：如果不是所有样本都有全部组学数据，需要决定是只用完全配对的子集，还是用能处理缺失的方法。MOFA2 支持部分缺失，SNF 不支持。

常见坑

坑 1：样本 ID 截断不一致导致匹配失败

TCGA 的 barcode 有 28 位，但蛋白组只保留前 12 位、甲基化用 Sentrix ID。直接 intersect 会得到 0 个匹配。解决方案是先建一张统一映射表（sample_id_mapping.csv），所有下游操作都基于这张表的 primary ID，永远不要在代码里硬编码截断规则。

坑 2：批次与生物学变量完全混杂

如果所有 RNA-seq 样本在 batch 1 做、所有蛋白组在 batch 2 做，那批次效应和组学差异完全混杂，无法校正。项目设计时必须让同一个样本的多层数据在相似时间处理，或者至少在每个 batch 内都有多个条件的样本。发现混杂后唯一的办法是承认这个局限性，而不是强行 ComBat。

坑 3：只用完全配对样本导致样本量暴跌

三层数据取交集后样本量可能从 100 降到 30。如果某一层缺失率高，考虑用支持部分缺失的方法（MOFA2）而不是一刀切只要完全配对。另一个策略是分步整合：先做两两整合，最后汇总结论。

坑 4：忽略缺失模式的非随机性

UpSet plot 显示某些样本组合的缺失不是随机的（比如所有 stage IV 患者都缺蛋白组数据），这意味着整合结果会对晚期患者产生偏差。在报告中必须说明缺失模式，并讨论对结论的影响。

下一步

接着深入：

• 02 各组学数据清洗与特征对应 — 每层数据进入整合前的预处理要点
• 03 跨层相关性探索 — 用简单相关分析检验数据是否"有信号"

横向延伸：

• bulk RNA-seq 01 数据获取与质控 — 转录组预处理的标准流程
• 蛋白质组 01 质谱结果表与实验设计 — 蛋白组数据的起点

参考资源

• MultiAssayExperiment Bioconductor
• MultiAssayExperiment 使用手册
• UpSetR 包
• TCGA barcode 说明

03 各组学数据清洗与特征对应

每一层组学数据在进入整合分析之前，都需要独立完成预处理。整合方法不会帮你处理批次效应、低质量探针或未归一化的计数值——垃圾进去，垃圾出来。

各层预处理要点

转录组（RNA-seq）

library(DESeq2)

dds <- DESeqDataSetFromMatrix(countData = rna_counts,
                              colData   = sample_info,
                              design    = ~ 1)
# 过滤低表达基因
keep <- rowSums(counts(dds) >= 10) >= 5
dds <- dds[keep, ]

# 方差稳定化（用于下游整合，不是差异分析）
vsd <- vst(dds, blind = TRUE)
rna_mat <- assay(vsd)

对于整合分析，通常用 VST 或 rlog 变换后的矩阵，而不是原始 counts 或 TPM。

甲基化（450K / EPIC）

library(minfi)

# 从 IDAT 读入
rgSet <- read.metharray.exp(targets = targets)
# 预处理：Noob 背景校正 + dye-bias 校正
mSet <- preprocessNoob(rgSet)
# 获取 beta 值
beta <- getBeta(mSet)

# 过滤：去掉 SNP 探针、cross-reactive 探针、性染色体探针
beta <- beta[!rownames(beta) %in% snp_probes, ]
beta <- beta[!rownames(beta) %in% cross_reactive, ]
beta <- beta[!grepl("^ch\\.", rownames(beta)), ]

蛋白质组

蛋白质组数据通常已经是 log2 intensity。主要处理缺失值和批次效应：

# 过滤缺失率 > 50% 的蛋白
miss_rate <- rowMeans(is.na(prot_mat))
prot_mat <- prot_mat[miss_rate < 0.5, ]

# KNN 填补
library(impute)
prot_mat <- impute.knn(prot_mat)$data

# 批次校正
library(limma)
prot_mat <- removeBatchEffect(prot_mat, batch = batch_info)

特征对应：从探针到基因

整合分析需要不同组学的特征能对应到同一个实体（通常是基因）。甲基化探针和蛋白 ID 都需要映射到 gene symbol。

甲基化探针 → 基因

library(IlluminaHumanMethylation450kanno.ilmn12.hg19)

anno <- getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)
probe_gene <- data.frame(
  probe = anno$Name,
  gene  = anno$UCSC_RefGene_Name,
  stringsAsFactors = FALSE
)
# 一个探针可能对应多个基因，取第一个或做 promoter 过滤
probe_gene$gene <- sapply(strsplit(probe_gene$gene, ";"), `[`, 1)
probe_gene <- probe_gene[probe_gene$gene != "", ]

基因级别聚合

一个基因可能对应多个甲基化探针。常见做法是取 promoter 区域探针的均值：

# 只保留 TSS200 和 TSS1500 区域的探针
promoter_probes <- anno$Name[grepl("TSS", anno$UCSC_RefGene_Group)]

beta_promoter <- beta[rownames(beta) %in% promoter_probes, ]

# 按基因聚合（取均值）
probe_gene_sub <- probe_gene[probe_gene$probe %in% rownames(beta_promoter), ]
beta_gene <- aggregate(beta_promoter[probe_gene_sub$probe, ],
                       by = list(gene = probe_gene_sub$gene),
                       FUN = mean)
rownames(beta_gene) <- beta_gene$gene
beta_gene$gene <- NULL

蛋白 → 基因

library(biomaRt)

ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
prot_ids <- rownames(prot_mat)  # UniProt ID

mapping <- getBM(attributes = c("uniprotswissprot", "hgnc_symbol"),
                 filters    = "uniprotswissprot",
                 values     = prot_ids,
                 mart       = ensembl)
# 合并
prot_mat_gene <- prot_mat[mapping$uniprotswissprot, ]
rownames(prot_mat_gene) <- mapping$hgnc_symbol

对齐后的检查

预处理和映射完成后，确认三层数据的维度和样本一致：

# 确保列（样本）顺序一致
common_samples <- intersect(intersect(colnames(rna_mat), colnames(beta_gene)),
                            colnames(prot_mat_gene))
rna_final  <- rna_mat[, common_samples]
meth_final <- beta_gene[, common_samples]
prot_final <- prot_mat_gene[, common_samples]

cat("样本数:", length(common_samples), "\n")
cat("RNA 特征:", nrow(rna_final), "\n")
cat("甲基化特征:", nrow(meth_final), "\n")
cat("蛋白特征:", nrow(prot_final), "\n")

常见坑

坑 1：归一化顺序搞反

正确顺序是"各层独立归一化 → 特征映射 → 整合"。如果先把多层数据拼接再做统一归一化，会让不同量纲的数据互相污染。RNA-seq 用 VST、甲基化用 Noob、蛋白组用 VSN，各走各的管线，最后才合并矩阵。

坑 2：探针到基因的一对多映射引入重复

一个基因可能对应 50 个甲基化探针，如果不聚合就直接进入整合矩阵，这个基因在甲基化层的权重会是其他基因的 50 倍。解决方案是先按 promoter 区域过滤，再对同基因探针取均值或取最变异的一个。

坑 3：biomaRt 版本不一致导致 ID 映射丢失

Ensembl 每年更新两次，不同版本的 gene symbol 映射会有差异。如果 RNA-seq 用 release 105 的注释，蛋白组用 release 110 的 biomaRt 查询，会丢失一批映射。建议固定一个 Ensembl release 版本，或者用 org.Hs.eg.db 这种相对稳定的本地注释包。

坑 4：对齐后没检查维度和样本顺序

cbind 或 merge 后如果列顺序不一致，会把 sample A 的转录组和 sample B 的蛋白组配在一起，后续所有结果都是错的。对齐后必须用 identical(colnames(mat1), colnames(mat2)) 验证，并打印样本数和特征数确认。

下一步

接着深入：

• 03 跨层相关性探索 — 清洗完的数据做第一轮信号检验
• 05 MOFA2 因子分析实战 — 对齐好的矩阵直接送入整合模型

横向延伸：

• bulk RNA-seq 06 高级建模 — 转录组预处理的更多细节
• 蛋白质组 05 缺失值策略 — 蛋白组特有的缺失值处理

参考资源

• DESeq2 vignette
• minfi 预处理流程
• biomaRt 使用指南
• impute 包

04 跨层相关性探索

在跑复杂的整合模型之前，先做一轮简单的相关性分析。它能告诉你：哪些基因在不同组学层之间有一致的信号，哪些层之间整体相关性高，以及数据里有没有明显的异常样本。

全局相关性：RV 系数

RV 系数衡量两个矩阵之间的整体相似度，取值 0-1，类似于矩阵层面的 Pearson 相关。

library(FactoMineR)

# 计算两两 RV 系数
rv_rna_meth <- coeffRV(t(rna_final), t(meth_final))$rv
rv_rna_prot <- coeffRV(t(rna_final), t(prot_final))$rv
rv_meth_prot <- coeffRV(t(meth_final), t(prot_final))$rv

cat("RNA vs Methylation RV:", round(rv_rna_meth, 3), "\n")
cat("RNA vs Protein RV:", round(rv_rna_prot, 3), "\n")
cat("Methylation vs Protein RV:", round(rv_meth_prot, 3), "\n")

通常 RNA 和蛋白的 RV 系数在 0.3-0.6 之间，甲基化和表达的相关性更低（因为甲基化主要是负调控，且只有 promoter 区域的甲基化与表达强相关）。

Cis 相关分析

Cis 相关指的是同一个基因在不同组学层之间的相关性。最经典的例子是：某个基因的 promoter 甲基化水平与它自身的 mRNA 表达之间的 Pearson 相关。

# 找到两层都有的基因
common_genes <- intersect(rownames(rna_final), rownames(meth_final))

# 逐基因计算相关
cis_cor <- sapply(common_genes, function(g) {
  cor(as.numeric(rna_final[g, ]),
      as.numeric(meth_final[g, ]),
      method = "pearson", use = "complete.obs")
})

# 分布
summary(cis_cor)
hist(cis_cor, breaks = 50, main = "Cis correlation: RNA vs Methylation",
     xlab = "Pearson r", col = "steelblue")

预期结果：大部分基因的 cis 相关接近 0，少部分呈现强负相关（promoter 甲基化抑制表达）。如果你看到大量正相关，可能是映射到了 gene body 而不是 promoter。

跨层相关热图

用 ComplexHeatmap 画一张跨层相关矩阵，展示 top 变异基因在不同组学层之间的关系：

library(ComplexHeatmap)
library(circlize)

# 选 top 500 变异基因（在 RNA 层）
rv <- apply(rna_final, 1, var)
top_genes <- names(sort(rv, decreasing = TRUE))[1:500]

# 只保留三层都有的
top_common <- Reduce(intersect, list(top_genes,
                                     rownames(meth_final),
                                     rownames(prot_final)))

# 计算 RNA vs Protein 的逐基因相关
rna_prot_cor <- sapply(top_common, function(g) {
  cor(as.numeric(rna_final[g, ]),
      as.numeric(prot_final[g, ]),
      use = "complete.obs")
})

# 热图：行是基因，列是相关类型
cor_mat <- cbind(
  RNA_Meth = cis_cor[top_common],
  RNA_Prot = rna_prot_cor[top_common]
)

col_fun <- colorRamp2(c(-1, 0, 1), c("blue", "white", "red"))
Heatmap(cor_mat, name = "Pearson r",
        col = col_fun,
        cluster_columns = FALSE,
        show_row_names = FALSE,
        column_title = "Cross-layer cis correlation")

Trans 相关的注意事项

Trans 相关（一个基因的甲基化与另一个基因的表达）在计算上是 O(n^2) 的，特征数多时不现实。常见策略：

• 只算已知调控关系的 pair（比如 TF 和 target）
• 先用 WGCNA 找模块，再算模块间相关
• 用 MOFA2 的因子来间接捕捉 trans 效应

相关性不等于因果

高相关不代表一层驱动了另一层。甲基化和表达的负相关可能是因果（甲基化沉默基因），也可能是共同受上游调控。如果需要因果推断，考虑 Mendelian randomization 或 mediation analysis。

常见坑

坑 1：cis 相关做了上万次检验却不校正

逐基因计算 cis 相关时，如果有 5000 个基因，就做了 5000 次检验。不做多重校正（BH/Bonferroni）就筛 p < 0.05，假阳性率会非常高。正确做法是对所有 cis 相关的 p 值统一做 FDR 校正，再设阈值。

坑 2：混杂变量导致虚假相关

年龄、性别、肿瘤纯度等混杂因素同时影响多层数据。不校正这些变量，计算出的 RNA-甲基化相关可能只是反映了"肿瘤纯度高的样本 RNA 高、甲基化也变了"。建议在相关分析前先用线性模型回归掉已知混杂，或者在 cor.test 中用偏相关。

坑 3：甲基化映射到 gene body 导致正相关"反直觉"

Promoter 甲基化与表达通常负相关，但 gene body 甲基化与表达正相关。如果在映射时没有限制到 TSS200/TSS1500 区域，混入了 gene body 探针，cis 相关的分布会出现大量正相关，容易误判为"甲基化激活基因表达"。

坑 4：样本量太小时 Pearson r 的置信区间巨大

20 个样本计算 Pearson 相关，即使 r = 0.5，95% CI 也是 [0.07, 0.78]。在小样本下不要过度解读单基因的相关系数，应该关注整体分布模式（比如 cis 相关的中位数是否显著偏离 0）。

不想本地装环境？在 BioF3 上跑

05 WGCNA 与共表达模块

WGCNA（Weighted Gene Co-expression Network Analysis）通过构建基因共表达网络，把上万个基因压缩成几十个模块，每个模块代表一组协同变化的基因。在多组学整合中，WGCNA 的模块可以作为"特征压缩"的手段，也可以用来检验某个模块在不同数据集或不同组学层之间是否保守。

基本流程

library(WGCNA)
allowWGCNAThreads()

# 输入：行是样本，列是基因（注意转置）
datExpr <- t(rna_final)

# 1. 选择软阈值（soft threshold）
powers <- c(1:20)
sft <- pickSoftThreshold(datExpr, powerVector = powers, verbose = 3)

# 画 scale-free topology fit 图
plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
     xlab = "Soft Threshold (power)",
     ylab = "Scale Free Topology Model Fit (signed R^2)",
     type = "n", main = "Scale independence")
text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
     labels = powers, col = "red")
abline(h = 0.85, col = "red")

选择 R^2 首次超过 0.85 的 power 值。通常 RNA-seq 数据在 6-12 之间。

构建网络与模块识别

# 2. 一步法构建网络
net <- blockwiseModules(datExpr,
                        power = 8,
                        TOMType = "unsigned",
                        minModuleSize = 30,
                        reassignThreshold = 0,
                        mergeCutHeight = 0.25,
                        numericLabels = TRUE,
                        pamRespectsDendro = FALSE,
                        verbose = 3)

# 模块颜色
moduleColors <- labels2colors(net$colors)
table(moduleColors)

# 画聚类树 + 模块颜色
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
                   "Module colors", dendroLabels = FALSE,
                   hang = 0.03, addGuide = TRUE, guideHang = 0.05)

模块与表型的关联

每个模块用 module eigengene（ME，即模块内基因表达的第一主成分）来代表。然后计算 ME 与临床表型的相关：

MEs <- net$MEs
# 计算 ME 与表型的相关
trait_data <- data.frame(
  subtype = as.numeric(factor(col_data$subtype)),
  stage   = as.numeric(factor(col_data$stage))
)

cor_ME_trait <- cor(MEs, trait_data, use = "p")
pval_ME_trait <- corPvalueStudent(cor_ME_trait, nrow(datExpr))

# 热图展示
library(ComplexHeatmap)
Heatmap(cor_ME_trait, name = "Correlation",
        cell_fun = function(j, i, x, y, w, h, fill) {
          if (pval_ME_trait[i, j] < 0.05)
            grid.text("*", x, y)
        })

模块保守性分析

如果你有两个独立数据集（比如 TCGA 和自己的队列），可以检验某个模块在第二个数据集中是否保守。这用 modulePreservation 函数：

# 准备第二个数据集
datExpr2 <- t(rna_validation)

# 设置多集数据
multiExpr <- list(
  Discovery  = list(data = datExpr),
  Validation = list(data = datExpr2)
)
multiColor <- list(Discovery = moduleColors)

# 计算保守性
mp <- modulePreservation(multiExpr, multiColor,
                         referenceNetworks = 1,
                         nPermutations = 200,
                         randomSeed = 1,
                         verbose = 3)

# Zsummary > 10 表示高度保守，2-10 中等，< 2 不保守
stats <- mp$preservation$Z$ref.Discovery$inColumnsAlsoPresentIn.Validation
print(stats[, c("moduleSize", "Zsummary.pres")])

在多组学整合中的角色

WGCNA 模块在整合分析中有两个用途：

1. 特征降维：用 ME 代替上千个基因，减少后续整合模型的输入维度。比如把 20 个模块的 ME 和蛋白组数据一起送进 MOFA2。
2. 跨层验证：在 RNA 层发现的模块，检查对应基因在蛋白层或甲基化层是否也有一致的模式。

常见坑

坑 1：soft threshold 选得太低或太高

R² 没达到 0.85 就随便选了个 power = 4，网络不满足 scale-free 假设，模块划分没有生物学意义。反之选 power = 20 会让网络太稀疏，大多数基因被归入灰色模块。一定要看 pickSoftThreshold 的拟合曲线，选 R² 首次超过 0.85 的最小 power。

坑 2：模块稳定性差却不做验证

WGCNA 的模块划分对参数（minModuleSize、mergeCutHeight）敏感。换个 seed 或调一下 mergeCutHeight 模块就完全不同，说明结果不稳定。应该用 modulePreservation 在独立数据集上验证，或者用 bootstrap 检验模块的可重复性。Zsummary < 2 的模块不应该报告。

坑 3：样本量不足导致相关矩阵不可靠

WGCNA 依赖基因间的相关矩阵。样本 < 20 时，相关系数的估计噪声很大，构建出的网络是假信号。官方建议至少 15 个样本，实际经验 30 以上效果才稳定。样本太少时考虑用其他降维方法（PCA、MOFA）替代。

坑 4：module eigengene 与表型的相关当成因果

ME 和 subtype 显著相关不代表"这个模块驱动了亚型形成"。ME 可能只是 confounding（比如该模块反映的是增殖速率，而增殖本身就和分期高度相关）。需要结合富集分析和文献验证，给出合理的生物学解释。

不想本地装环境？在 BioF3 上跑

06 MOFA2 因子分析实战

MOFA2（Multi-Omics Factor Analysis v2）用概率矩阵分解从多层组学数据中提取潜因子。每个因子捕捉一种跨组学的共同变异模式，可以理解为"多组学版的 PCA"。本节用 MOFA2 自带的 CLL（慢性淋巴细胞白血病）数据集走一遍完整流程。

library(MOFA2)
library(ggplot2)

# 载入 CLL 示例数据
file <- system.file("extdata", "CLL_data.hdf5", package = "MOFA2")
model <- load_model(file)

# 查看模型基本信息
model

# 从矩阵列表创建 MOFA 对象
data_list <- list(
  mRNA       = rna_final,
  Methylation = meth_final,
  Protein    = prot_final
)

mofa_obj <- create_mofa(data_list)

# 设置参数
data_opts   <- get_default_data_options(mofa_obj)
model_opts  <- get_default_model_options(mofa_obj)
model_opts$num_factors <- 15

train_opts  <- get_default_training_options(mofa_obj)
train_opts$convergence_mode <- "slow"
train_opts$seed <- 42

mofa_obj <- prepare_mofa(mofa_obj,
                         data_options    = data_opts,
                         model_options   = model_opts,
                         training_options = train_opts)

# 训练（需要 Python 环境中安装 mofapy2）
model <- run_mofa(mofa_obj, outfile = "mofa_model.hdf5")

# 热图：行是因子，列是组学层
plot_variance_explained(model, x = "view", y = "factor")

# 每个组学层被所有因子解释的总方差
plot_variance_explained(model, plot_total = TRUE)[[2]]

# Factor 1 在 mRNA 层的 top 权重特征
plot_top_weights(model, view = "mRNA", factor = 1, nfeatures = 15)

# 用 Factor 1 和 Factor 2 画散点图，按亚型着色
plot_factor(model, factors = c(1, 2), color_by = "IGHV")

# 因子之间的相关性（理想情况下应接近 0）
plot_factor_cor(model)

# 展示各层数据的缺失模式和样本覆盖
plot_data_overview(model)

library(msigdbr)
library(fgsea)

# 获取 Factor 1 在 mRNA 层的权重
weights <- get_weights(model, views = "mRNA", factors = 1, as.data.frame = TRUE)
gene_ranks <- setNames(weights$value, weights$feature)
gene_ranks <- sort(gene_ranks, decreasing = TRUE)

# 用 Hallmark 基因集做 GSEA
hallmark <- msigdbr(species = "Homo sapiens", category = "H")
pathways <- split(hallmark$gene_symbol, hallmark$gs_name)

fgsea_res <- fgsea(pathways, gene_ranks, minSize = 15, maxSize = 500)
head(fgsea_res[order(pval), ], 10)

07 SNF 相似性网络与亚型识别

SNF（Similarity Network Fusion）的思路和 MOFA 完全不同：它不做矩阵分解，而是为每一层组学数据分别构建样本相似性网络，然后把这些网络融合成一个统一的网络。融合后的网络可以用谱聚类（spectral clustering）来识别亚型。

算法直觉

1. 对每一层组学数据，计算样本间的距离矩阵（通常用欧氏距离）
2. 把距离矩阵转换成相似性矩阵（用 scaled exponential kernel）
3. 构建 K 近邻图（KNN graph）
4. 迭代融合：每一层的网络向其他层"扩散"信息，最终收敛到一个融合网络

融合的核心思想是：如果两个样本在多层数据中都相似，它们在融合网络中的连接会被加强；如果只在一层相似，连接会被削弱。

基本流程

library(SNFtool)

# 输入：每层数据矩阵，行是样本，列是特征
rna_t  <- t(rna_final)
meth_t <- t(meth_final)
prot_t <- t(prot_final)

# 参数
K <- 20          # 近邻数
alpha <- 0.5     # 超参数
iter <- 20       # 迭代次数

# 1. 计算距离矩阵
dist_rna  <- dist2(as.matrix(rna_t), as.matrix(rna_t))
dist_meth <- dist2(as.matrix(meth_t), as.matrix(meth_t))
dist_prot <- dist2(as.matrix(prot_t), as.matrix(prot_t))

# 2. 构建相似性网络
W_rna  <- affinityMatrix(dist_rna, K, alpha)
W_meth <- affinityMatrix(dist_meth, K, alpha)
W_prot <- affinityMatrix(dist_prot, K, alpha)

# 3. 融合网络
W_fused <- SNF(list(W_rna, W_meth, W_prot), K, iter)

谱聚类确定亚型

# 用 eigen-gap 方法估计最优聚类数
C_best <- estimateNumberOfClustersGivenGraph(W_fused, NUMC = 2:8)
cat("推荐聚类数:", C_best$`Eigen-gap best`, "\n")

# 谱聚类
clusters <- spectralClustering(W_fused, K = C_best$`Eigen-gap best`)
table(clusters)

# 可视化融合网络
library(pheatmap)
pheatmap(W_fused,
         annotation_row = data.frame(Cluster = factor(clusters)),
         show_rownames = FALSE, show_colnames = FALSE,
         main = "Fused similarity network")

与 Consensus Clustering 比较

SNF 和 consensus clustering（如 ConsensusClusterPlus）都能做亚型发现，但思路不同：

实际项目中可以两种都跑，看结果是否一致。如果一致，结论更可信。

library(ConsensusClusterPlus)

# 拼接矩阵
combined <- rbind(rna_final, meth_final, prot_final)

# Consensus clustering
cc_results <- ConsensusClusterPlus(combined,
                                   maxK = 8,
                                   reps = 500,
                                   pItem = 0.8,
                                   pFeature = 1,
                                   clusterAlg = "hc",
                                   distance = "pearson",
                                   seed = 42,
                                   plot = "pdf")

亚型与临床表型关联

# 把聚类结果和临床信息合并
cluster_df <- data.frame(
  sample  = rownames(rna_t),
  cluster = factor(clusters),
  subtype = col_data$subtype,
  stage   = col_data$stage
)

# Fisher 检验：聚类 vs 已知亚型
fisher.test(table(cluster_df$cluster, cluster_df$subtype))

# 生存分析
library(survival)
library(survminer)

surv_obj <- Surv(col_data$os_time, col_data$os_event)
fit <- survfit(surv_obj ~ cluster_df$cluster)
ggsurvplot(fit, data = cluster_df, pval = TRUE,
           palette = "jco", risk.table = TRUE)

参数敏感性

SNF 对 K（近邻数）比较敏感。建议在 K = 10-30 范围内扫描，看聚类结果是否稳定。alpha 通常固定为 0.5，迭代次数 20 次足够收敛。

常见坑

坑 1：K 值敏感但没做参数扫描

SNF 对近邻数 K 非常敏感。K = 10 和 K = 30 可能给出完全不同的聚类。很多人只跑一次 K = 20 就报告结果。正确做法是在 K = 10-30 范围内扫描，看聚类标签的 NMI（normalized mutual information）是否稳定。如果不同 K 结果差异大，结论不可信。

坑 2：样本量太小时网络融合不可靠

SNF 构建 KNN 图时，如果样本数 < 30 而 K = 20，几乎每个样本都是彼此的邻居，网络失去区分力。小样本场景建议把 K 降到 sample_size / 5 以下，或者换用 MOFA2（对小样本更友好）。

坑 3：某一层数据质量差污染融合结果

SNF 会把噪声层的信息也融合进去。如果蛋白组只有 500 个特征且缺失率高，它的相似性矩阵很嘈杂，融合后会拖累其他层的信号。建议先单独检查每层的 silhouette score，质量太差的层考虑排除。

坑 4：用 eigen-gap 选聚类数时只看一个指标

estimateNumberOfClustersGivenGraph 返回多个指标（eigen-gap、rotation cost 等），它们可能给出不同的最优 K。不要只看一个。如果多个指标一致推荐 K = 3，结果更可信；如果分歧大，说明数据中没有清晰的亚型结构。

下一步

接着深入：

• 07 机器学习预测模型 — 从亚型发现转向临床预测
• 08 整合结果的生物学解读 — 把聚类标签翻译成生物学亚型

横向延伸：

• 05 MOFA2 因子分析实战 — 另一种整合方法，适合小样本和部分缺失
• bulk RNA-seq 07 聚类与亚型 — 单组学聚类方法的对比

参考资源

• SNFtool CRAN
• SNF 原始论文 (Wang et al. 2014)
• ConsensusClusterPlus Bioconductor
• spectral clustering 原理

08 机器学习预测模型

前面几个模块侧重"发现结构"（因子、模块、亚型），本模块转向"预测"：用多组学特征预测临床结局（如药物响应、生存状态、疾病亚型）。核心问题是如何从多层高维数据中选出有预测力的特征子集，以及如何设计合理的交叉验证策略。

特征选择策略

多组学数据的特征数远大于样本数（p >> n），直接建模会过拟合。常见的特征选择路径：

1. 单层筛选后拼接：每层独立做方差过滤或差异分析，保留 top 特征，再拼接
2. 用整合结果做特征：MOFA 因子得分、WGCNA 模块 eigengene 作为输入
3. 嵌入式选择：用 elastic net 或 Random Forest 的内置特征重要性

# 方法 1：每层取 top 500 高变异特征
select_top <- function(mat, n = 500) {
  rv <- apply(mat, 1, var)
  mat[names(sort(rv, decreasing = TRUE))[1:n], ]
}

rna_sel  <- select_top(rna_final, 500)
meth_sel <- select_top(meth_final, 500)
prot_sel <- select_top(prot_final, 200)

# 拼接
X <- t(rbind(rna_sel, meth_sel, prot_sel))
y <- factor(col_data$subtype)

Elastic Net 分类

Elastic net 结合了 L1（特征选择）和 L2（稳定性）正则化，适合高维小样本场景：

library(glmnet)

# alpha = 0.5 是 elastic net；alpha = 1 是 lasso
set.seed(42)
cv_fit <- cv.glmnet(X, y, family = "multinomial",
                    alpha = 0.5, nfolds = 5)

# 最优 lambda
plot(cv_fit)
best_lambda <- cv_fit$lambda.min

# 非零系数（被选中的特征）
coef_list <- coef(cv_fit, s = best_lambda)
selected_features <- lapply(coef_list, function(c) {
  rownames(c)[which(c[, 1] != 0)]
})

Random Forest

Random Forest 不需要特征预筛选，能处理非线性关系，并提供特征重要性排序：

library(randomForest)

set.seed(42)
rf_model <- randomForest(X, y, ntree = 1000, importance = TRUE)

# OOB 错误率
print(rf_model)

# 特征重要性
imp <- importance(rf_model, type = 1)
top_imp <- head(sort(imp[, 1], decreasing = TRUE), 30)
barplot(top_imp, las = 2, main = "Top 30 features by importance",
        cex.names = 0.6)

交叉验证设计

多组学预测模型最容易犯的错误是信息泄露（data leakage）。特征选择必须在交叉验证的内层完成，不能用全部数据选完特征再做 CV。

library(caret)

# 正确做法：用 caret 的嵌套 CV
ctrl <- trainControl(method = "repeatedcv",
                     number = 5,
                     repeats = 10,
                     classProbs = TRUE,
                     summaryFunction = multiClassSummary)

# 在每个 fold 内部做特征选择 + 建模
set.seed(42)
rf_cv <- train(X, y,
               method = "rf",
               trControl = ctrl,
               tuneLength = 5,
               metric = "AUC")
print(rf_cv)

如果样本量很小（< 50），考虑 leave-one-out CV（LOOCV）或 repeated 5-fold CV 来减少方差。

多组学 vs 单组学的预测力比较

一个关键问题是：多组学整合是否真的比单组学预测更好？需要做对照实验：

# 分别用单层数据训练
rf_rna  <- train(t(rna_sel), y, method = "rf", trControl = ctrl, metric = "AUC")
rf_meth <- train(t(meth_sel), y, method = "rf", trControl = ctrl, metric = "AUC")
rf_prot <- train(t(prot_sel), y, method = "rf", trControl = ctrl, metric = "AUC")

# 比较
results <- resamples(list(
  MultiOmics = rf_cv,
  RNA_only   = rf_rna,
  Meth_only  = rf_meth,
  Prot_only  = rf_prot
))
summary(results)
dotplot(results, metric = "AUC")

如果多组学没有显著提升，说明信息冗余或者某一层已经足够。这本身也是有价值的结论。

生存预测

如果目标是生存时间而不是分类，用 Cox 回归配合 elastic net：

library(glmnet)

surv_y <- Surv(col_data$os_time, col_data$os_event)

cv_cox <- cv.glmnet(X, surv_y, family = "cox",
                    alpha = 0.5, nfolds = 5)
# C-index
max(cv_cox$cvm)

常见坑

坑 1：特征选择在交叉验证外面做（信息泄露）

最经典的错误：用全部数据选完 top 500 高变异特征，再做 5-fold CV。这样 CV 的准确率会虚高，因为 test fold 的信息已经泄露到特征选择步骤。正确做法是嵌套 CV：外层评估性能，内层同时做特征选择 + 模型调参。

坑 2：样本量 < 50 时用复杂模型

30 个样本 + 1500 个特征跑 Random Forest 或深度学习，OOB 看着不错但外部验证全崩。p >> n 场景下，elastic net（L1 + L2 正则化）比 tree-based 方法更安全。如果样本实在太少，考虑用 MOFA 因子得分（15 维）替代原始特征。

坑 3：多组学 vs 单组学比较不公平

多组学模型用了 1200 个特征，单组学模型只用了 400 个。特征数不同本身就会影响性能。公平比较应该控制总特征数一致，或者用 AUC/C-index 配合置信区间做统计检验（DeLong test），而不是只看点估计。

坑 4：没有独立验证集

在同一个数据集上调参 + 评估，即使用了 CV 也可能过拟合超参数。如果有条件，保留 20-30% 数据完全不碰，最终只用一次做 blind test。或者用不同来源的数据集（如 TCGA 训练、CPTAC 验证）。

下一步

接着深入：

• 08 整合结果的生物学解读与报告 — 把预测模型的结果写进论文
• 06 SNF 相似性网络与亚型识别 — 从预测转向无监督亚型发现

横向延伸：

• 05 MOFA2 因子分析实战 — 用 MOFA 因子得分作为 ML 模型输入
• bulk RNA-seq 08 报告与可复现 — 单组学预测模型的交付规范

参考资源

• glmnet vignette
• caret 包文档
• randomForest 包
• 交叉验证中的信息泄露 (Hastie et al. ESL Ch.7)
• 多组学预测综述 (Picard et al. 2021)

09 整合结果的生物学解读与报告

跑完 MOFA、SNF 或机器学习模型之后，你手上有一堆因子、模块、聚类标签和特征列表。这些数字本身不是结论——需要把它们翻译成生物学语言，并以可复现的方式写进论文。

从因子到通路

MOFA 因子的权重列表本质上是一个 ranked gene list，可以直接送进 GSEA：

library(MOFA2)
library(fgsea)
library(msigdbr)

# 获取 Factor 1 在 mRNA 层的权重
weights <- get_weights(model, views = "mRNA", factors = 1, as.data.frame = TRUE)
gene_ranks <- setNames(weights$value, weights$feature)
gene_ranks <- sort(gene_ranks, decreasing = TRUE)

# Hallmark 基因集
hallmark <- msigdbr(species = "Homo sapiens", category = "H")
pathways <- split(hallmark$gene_symbol, hallmark$gs_name)

fgsea_res <- fgsea(pathways, gene_ranks, minSize = 15, maxSize = 500)
sig_pathways <- fgsea_res[padj < 0.05][order(NES, decreasing = TRUE)]
head(sig_pathways[, .(pathway, NES, padj)], 10)

如果 Factor 1 富集在 cell cycle 和 DNA repair 通路，而且它在 Mutations 层也有高方差贡献，你可以说"Factor 1 捕捉了与基因组不稳定性相关的跨组学变异"。

从聚类到亚型特征

SNF 或 consensus clustering 给出的聚类标签需要进一步表征：

# 每个亚型的差异特征
library(limma)

design <- model.matrix(~ 0 + factor(clusters))
colnames(design) <- paste0("C", 1:ncol(design))

fit <- lmFit(rna_final, design)
contrast_mat <- makeContrasts(C1 - C2, C1 - C3, C2 - C3, levels = design)
fit2 <- contrasts.fit(fit, contrast_mat)
fit2 <- eBayes(fit2)

# C1 vs C2 的 top 差异基因
topTable(fit2, coef = 1, number = 20)

把每个亚型的 marker 基因列出来，再做 GO/KEGG 富集，就能给亚型一个生物学标签（比如"免疫活跃型""代谢重编程型"）。

多层证据汇总

好的多组学分析不是把每层结果分别报告，而是展示"多层证据指向同一个结论"。一个常见的汇总方式：

# 某个基因在三层的一致性
gene <- "TP53"

# RNA 层：在亚型间差异表达
rna_p <- t.test(rna_final[gene, ] ~ clusters)$p.value

# 甲基化层：promoter 甲基化差异
meth_p <- t.test(meth_final[gene, ] ~ clusters)$p.value

# 蛋白层：蛋白丰度差异
prot_p <- t.test(prot_final[gene, ] ~ clusters)$p.value

cat(gene, "across layers:\n")
cat("  RNA p-value:", format(rna_p, digits = 3), "\n")
cat("  Methylation p-value:", format(meth_p, digits = 3), "\n")
cat("  Protein p-value:", format(prot_p, digits = 3), "\n")

写 Methods 段

多组学论文的 Methods 段需要覆盖以下内容。下面是一个模板：

Multi-omics integration was performed using MOFA2 (v1.8.0).
Input data included transcriptomics (N genes after filtering),
DNA methylation (N CpG sites, promoter-level aggregation),
and proteomics (N proteins). Each layer was independently
preprocessed: RNA-seq counts were variance-stabilized using
DESeq2; methylation beta values were Noob-normalized with minfi;
protein intensities were log2-transformed and batch-corrected
with limma::removeBatchEffect. Features were mapped to gene
symbols using biomaRt. The model was trained with K=15 factors,
slow convergence mode, and random seed 42. Downstream enrichment
analysis used fgsea with MSigDB Hallmark gene sets.

关键原则：

• 写清楚每层的预处理方法和工具版本
• 说明特征映射策略（probe → gene 用了什么注释）
• 报告整合方法的参数（因子数、迭代次数、收敛标准）
• 提供代码仓库链接

可视化建议

多组学论文常见的图：

library(ComplexHeatmap)

# 多层 heatmap：同一批样本，上面是 RNA，中间是甲基化，下面是蛋白
ht1 <- Heatmap(rna_final[top_genes, order(clusters)], name = "RNA",
               show_column_names = FALSE, cluster_columns = FALSE)
ht2 <- Heatmap(meth_final[top_genes, order(clusters)], name = "Meth",
               show_column_names = FALSE, cluster_columns = FALSE)
ht3 <- Heatmap(prot_final[top_genes, order(clusters)], name = "Prot",
               show_column_names = FALSE, cluster_columns = FALSE)

ht_list <- ht1 %v% ht2 %v% ht3
draw(ht_list, column_title = "Multi-omics heatmap by subtype")

可复现性清单

发表前检查：

• 原始数据是否上传到公共数据库（GEO、PRIDE、GDC）
• 分析代码是否在 GitHub/Zenodo 上公开
• R session info 是否记录（sessionInfo()）
• 随机种子是否固定
• 中间结果（模型文件、聚类标签）是否保存

常见坑

坑 1：过度解读单个因子的通路富集

Factor 1 富集到 cell cycle（padj = 0.03），就写"多组学分析揭示了 cell cycle 是关键驱动通路"。但一个 padj 接近阈值的富集可能换个基因集版本就消失了。应该结合多个证据（多个因子、多层数据、外部数据集验证）才能下结论。

坑 2：聚类标签没做外部验证

SNF 给出 3 个亚型，但没有用独立队列验证这些亚型是否可复现。发表前至少应该：(1) 在发现集内做 bootstrap 检验 silhouette score；(2) 有条件时在验证集上重新跑 SNF 看聚类是否一致。

坑 3：可复现性材料不完整

代码上传了但没固定随机种子、没记录包版本、没保存中间模型文件。审稿人或读者无法从代码重现结果。建议用 renv 锁定包版本，sessionInfo() 写进报告，模型 .hdf5 文件存到 Zenodo。

坑 4：Methods 段漏写关键参数

只写"用 MOFA2 做了多组学整合"而没交代因子数、收敛标准、各层的预处理细节和特征数。审稿人第一轮就会要求补充。Methods 模板应覆盖：每层独立预处理方法和版本、特征映射策略、整合方法参数、下游富集所用基因集。

下一步

接着深入：

• 05 MOFA2 因子分析实战 — 回顾因子分析的技术细节
• 07 机器学习预测模型 — 从解读转向预测建模

横向延伸：

• bulk RNA-seq 08 报告与可复现 — 单组学项目的报告规范可直接复用
• 蛋白质组 03 功能富集 — 富集分析的详细操作和解读

参考资源

• fgsea Bioconductor
• msigdbr 基因集
• ComplexHeatmap 完整手册
• ggalluvial 包
• 多组学论文写作指南 (Subramanian et al. 2020)