最近帮一个课题组处理多组学数据样本本身很珍贵同一批组织既跑了WES又做了单细胞核RNA测序和ATAC测序还送了蛋白质组、空间转录组和空间代谢组。当时最大的困境不是单个组学不会分析而是五个组学各自的报告堆在一起之后没办法形成一条连贯的故事线。这篇文章就想把这些环节重新捋一遍——从WES到snRNA-seq、snATAC-seq再到蛋白质组、空间转录组和空间代谢组每个模块我都会给出真正能跑的分析代码和出图逻辑同时讲清楚为什么这一步要这样做、和下一步怎么接上。这篇文章更适合两类人看一类是刚接触多组学项目、手里的数据已经到齐但不知道从哪下手的分析同学另一类是只做过单一组学、想往整合方向扩展的研究者。整篇文章我尽量用“拿到数据后我会怎么操作”的方式来讲代码都是可以直接替换成自己数据路径去跑的版本而不是那种只放个函数名就结束的伪代码。1. 多条数据放一起前先把科学问题的分析主线定下来1.1 多组学不是多个单组学报告的拼接很多多组学项目最后失败不是测序质量差而是每个分析环节各自为政。WES找到几个驱动突变单细胞分出十几个亚群蛋白组筛出几十个差异蛋白空间转录组图上看起来有区域结构代谢组图像也漂亮——但放到一起谁和谁有关系谁也说不清。我习惯在正式动代码之前先画一张整合草图。这张草图上只需要回答三个问题样本来自什么样的队列、每个组学分别能回答哪一层问题、层与层之间在哪个尺度上能对上。比如同一批组织切片用于空间转录组和空间代谢组时具体是哪一张切片、什么方向包埋、切片序号是否一致这些如果不在分析前记录清楚后面做共定位分析时只能靠猜。我遇到过一个实际场景基因表达层面看到某个代谢酶的mRNA在特定区域显著升高但蛋白组和代谢组结果里这个通路并没有激活。后来排查发现mRNA升高和蛋白表达之间存在时间差和翻译调控单看转录组就会得出误导性结论。这就是多组学的价值——用不同维度的数据互相约束、互相验证而不是各自讲一个孤立的故事。1.2 一张分析路线表让每组数据都服务于最终故事组学层次数据形态核心分析主要输出图回答的科学问题WESFASTQ → BAM → VCF/MAF突变检测与注释oncoplot、lollipop图哪些基因携带驱动突变snRNA-seq表达矩阵聚类、细胞类型注释UMAP、DotPlot、feature plot组织中有哪些细胞类型和状态snATAC-seq片段文件/fragmentspeak calling、motif富集峰可及性图、motif富集图哪些调控元件在这些细胞中开放蛋白质组蛋白定量矩阵差异分析、富集分析热图、火山图、富集气泡图哪些蛋白真正发生了变化空间转录组空间表达矩阵空间聚类、区域marker空间特征图、空间聚类图基因表达在组织中的空间分布如何空间代谢组质谱成像数据/imzML单离子成像、空间相关代谢物空间分布图、相关热图代谢物在哪里富集与转录空间是否一致这张表会贯穿整个分析周期。每做完一个组学都要回头问一句这个结果对整体的科学假设是支持、反驳还是新增视角。我在给学生的建议里经常说多组学分析最大的风险不是缺少工具而是缺少一个把数据串起来的主线。主线一旦清晰每个分析模块的优先级也就清楚了。2. WES先把体细胞突变的“候选名单”筛干净2.1 WES 标准分析管线中的核心步骤WES的核心是捕获外显子区域并测序目的是找到和疾病相关的体细胞突变。拿到FASTQ后常规流程是质控过滤、比对到参考基因组、去除PCR重复、突变检测、突变注释。每一步的软件选择已经相当成熟我这边给一套能直接跑的短流程重点不是把命令铺全而是标注出每个步骤里最容易出问题的地方。# 1. 质控与接头修剪fastp 一步完成 fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \ -o sample_trim_R1.fq.gz -O sample_trim_R2.fq.gz \ --html sample.fastp.html \ --thread 16 # 2. BWA-MEM 比对 排序 建索引 bwa mem -t 16 hg38.fa sample_trim_R1.fq.gz sample_trim_R2.fq.gz | \ samtools sort - 16 -o sample.sorted.bam samtools index sample.sorted.bam # 3. 标记并去除PCR重复最好用 GATK MarkDuplicates也可以直接用 picard gatk MarkDuplicates \ -I sample.sorted.bam \ -O sample.dedup.bam \ -M sample.dedup.metrics.txt # 4. 突变检测肿瘤/正常配对用 Mutect2单样本则用 HaplotypeCaller gatk Mutect2 \ -R hg38.fa \ -I tumor.bam \ -I normal.bam \ -tumor TUMOR_SAMPLE \ -normal NORMAL_SAMPLE \ -O sample.somatic.vcf.gz # 5. 用 Funcotator 或 ANNOVAR/VEP 做突变注释 gatk Funcotator \ -R hg38.fa \ -V sample.somatic.vcf.gz \ -O sample.somatic.maf \ --ref-version hg38 \ --data-sources-path /path/to/funcotator_dataSources.v1.8.hg38.20230218/一个非常常见的坑是建索引不完整。包括参考基因组的.dict、.fai、bam 的.bai缺一个都会在中间步骤报错而且报错信息经常让人看不懂。建议在拿到参考基因组后第一时间把索引都建好。另外Mutect2 的肿瘤/正常配对中如果 normal 样本本身也有杂质过滤 threshold 可以适当放宽但如果是公开数据分析缺少正常配对时可以用 panel of normalsPON来过滤种系突变否则结果里会混入大量来自个体本身的 SNP。2.2 突变全景图和驱动基因的几种可视化解法突变注释完成后MAF 文件就是后续可视化的基础。R 里我常用maftools它可以直接读 MAF画 oncoplot瀑布图、lollipop 图、TMB 比较等。这些图在文章里几乎是标配代码也很简单library(maftools) # 读入 MAF 文件同时可以带上临床信息表 maf - read.maf( maf sample.somatic.maf, clinicalData clinical.tsv ) # 瀑布图展示每个样本的突变谱 oncoplot( maf, top 20, drawRowBar TRUE, drawColBar TRUE, showTumorSampleBarcodes FALSE ) # 单个基因的突变位点分布图适合展示驱动基因 lollipopPlot( maf, gene TP53, AACol Protein_Change, showDomain TRUE )很多人直接拿top 20画图但这里有个细节top默认按突变频率排序如果你的队列里存在一个突变极多的样本瀑布图会被这个样本拉得很长其他样本的突变信息被压缩。我通常会在画图前先按样本突变总数排序再考虑是否要对 top 基因做二次筛选。如果队列样本数超过 50还可以叠加tmb()函数计算肿瘤突变负荷结合临床信息做分组比较。WES 层面的分析目标并不复杂——找到可信的驱动事件为后续转录组、蛋白组层面的验证提供“上游候选基因”清单。只有在 WES 阶段把突变名单筛干净后面做突变-表达关联时才不会把时间浪费在噪音上。3. snRNA-seq细胞聚类和单细胞手动注释是最花时间的环节3.1 Seurat 标准流程从表达矩阵到 UMAPsnRNA-seq 的分析主力目前还是 Seurat。核 RNA 数据和普通单细胞数据的区别在于核内捕获到的转录本更多偏向未成熟 mRNA 和长链非编码RNA胞浆 mRNA 比例下降所以某些细胞类型比如中性粒细胞在核 RNA 数据里会比较难注释。分析流程上差异不大但基因过滤阈值和 marker 选择需要针对核数据做调整。library(Seurat) library(dplyr) # 读入 10X 输出目录或直接读表达矩阵 obj - CreateSeuratObject( counts Read10X(filtered_feature_bc_matrix/), project multiomics, min.cells 3, min.features 200 ) # 质控线粒体比例和基因数过滤 obj[[percent.mt]] - PercentageFeatureSet(obj, pattern ^MT-) obj - subset(obj, subset nFeature_RNA 200 nFeature_RNA 6000 percent.mt 10) # 标准化、高变基因、PCA、UMAP、聚类 obj - NormalizeData(obj) obj - FindVariableFeatures(obj, selection.method vst, nfeatures 3000) obj - ScaleData(obj) obj - RunPCA(obj, npcs 30) obj - RunUMAP(obj, dims 1:20) obj - FindNeighbors(obj, dims 1:20) obj - FindClusters(obj, resolution 0.8) # 保存聚类结果后续手动注释会反复用到 DimPlot(obj, reduction umap, label TRUE)关于聚类分辩率的设置0.8 在大多数组织里是一个比较常用的起点但并没有一个绝对最优值。实际操作中我会用 0.4 到 1.2 的范围各跑一遍然后检查 cluster 的稳定性一个合理的聚类结果应该是分辨率小幅变化时大部分 cluster 能保持稳定而不是被切碎或完全合并。单细胞分析的“最佳分辨率”本质上是由下游注释需求决定的想把巨噬细胞M1/M2分开分辨率就要高一些只想分大类0.4 足够。3.2 单细胞手动注释自动化注释工具只是参考替代不了人工确认现在很多人喜欢用 SingleR、CellTypist 这类自动注释工具结果一跑出来就直接用这是比较危险的习惯。自动注释的核心依赖参考数据集一旦你的组织和参考数据来源差异大比如疾病状态、物种、组织部位不同错的概率会明显上升。我个人的流程是先做自动注释得到一个初步方向然后进入“单细胞手动注释”流程——直接用已知的 canonical markers 在所有 cluster 上打分逐个检查。# 手动注释的核心用已知 marker 列表做模块打分 markers - list( CD4 T c(CD3D, CD4, IL7R, LEF1), CD8 T c(CD3D, CD8A, CD8B, GZMB, NKG7), NK c(NCAM1, NKG7, KLRD1, KLRF1), B c(CD79A, MS4A1, CD19), Plasma c(MZB1, SDC1, JCHAIN, XBP1), Myeloid c(LYZ, CD68, CD163, AIF1), Endothelial c(PECAM1, VWF, CLDN5, FLT1), Fibroblast c(COL1A1, COL1A2, DCN, PDGFRA), Epithelial c(EPCAM, KRT8, KRT18, KRT19) ) # 直接为每个 cluster 计算这些 marker 集的模块得分 obj - AddModuleScore(obj, features markers, name ManualScore) # 用热图/小提琴图看每个 cluster 对应哪个 marker 集 VlnPlot(obj, features ManualScore1, pt.size 0) DotPlot(obj, features unique(unlist(markers)), group.by seurat_clusters) RotatedAxis()建议不要完全依赖打分结果做自动分配。正确的手动注释姿势是先看模块得分高的细胞群分布再回到 UMAP 上调几个典型 marker 的 FeaturePlot确认空间位置是否符合生物学直觉同时对边界不清的 cluster单独拉出来重新聚类看亚群。这个流程听起来费时间但一个 cluster 注释错了后面所有差异分析和组成分析都会受影响。尤其是髓系细胞这个“垃圾筐”几乎所有组织里都会存在一个marker表达不高不低、难以归类的 cluster仔细观察后通常会发现是双细胞或过渡态细胞需要按实际表达特征单独命名而不能简单标记为“unknown”就放弃。4. snATAC-seq从染色质开放区域看基因调控线索4.1 ArchR 处理 snATAC-seq 的标准分析姿势snATAC-seq 和 snRNA-seq 经常来自同一个核悬液也就是常说的 Multiome 实验同一个细胞核里同时捕获 RNA 和开放染色质信息。分析染色质可及性的主流工具是 ArchR 和 Signac两者各有优势。ArchR 在处理大样本时速度更快内置了 topic modelIterative LSI来做降维Signac 的优势是能和 Seurat 生态无缝衔接。我的习惯是如果数据来自 10X Multiome优先用 ArchR 做完整的 snATAC 流程再把结果降到 Seurat 里做多组学整合。library(ArchR) addArchRGenome(hg38) addArchRThreads(threads 16) # 输入是 10X 输出的 fragments.tsv.gz 文件 ArrowFiles - createArrowFiles( inputFiles c(sample1_fragments.tsv.gz, sample2_fragments.tsv.gz), sampleNames c(sample1, sample2), minTSS 4, minFrags 1000, maxFrags 5e5 ) proj - ArchRProject(ArrowFiles) proj - addIterativeLSI(proj, useMatrix TileMatrix, name IterativeLSI) proj - addClusters(proj, reducedDims IterativeLSI, resolution 1.0) proj - addUMAP(proj, reducedDims IterativeLSI) # 检查 TSS 富集分数评估数据质量 getTSSEnrichment(proj, fast TRUE)TSS 富集分数是判断 snATAC 数据质量的重要指标低于 4 的样本建议直接过滤掉。实际操作中我遇到过一批样本碎片数非常高但 TSS 分数很低后来发现是细胞核裂解不充分导致背景噪音太高换一批核重新建库后数据就好了。所以拿到数据第一步先看 TSS 分数再决定是否进入下游分析。调用 peak 之后最有价值的信息是不同细胞类型的特化调控元件。我会在 ArchR 里用getMarkerFeatures()找到每个 cluster 特异的 peak再做 motif 富集看哪些转录因子结合位点在这些开放区域中显著富集。这一步能把“哪个基因在哪个细胞类型里可能被调控”这个问题的线索找出来。# 找每个 cluster 特异开放的 peak markersPeaks - getMarkerFeatures( ArchRProj proj, useMatrix PeakMatrix, groupBy Clusters, bias c(TSSEnrichment, log10(nFrags)), testMethod wilcoxon ) # 对特异 peak 做 motif 富集找出关键的转录因子 proj - addMotifAnnotations(proj, motifSet cisbp) proj - addBGD(proj, type Motif) enrichMotifs - peakAnnoEnrichment( ArchRProj proj, peakSet markersPeaks, background proj, motifSet cisbp )4.2 把 snATAC 和 snRNA 放进同一个坐标系当同一个细胞核同时有 RNA 和 ATAC 数据时整合分析通常有两种策略一种是把 snATAC 的 peak 矩阵映射到 snRNA 的细胞嵌入空间用 Seurat 的WNN加权最近邻方法进行联合降维另一种是用 MOFA 这类多因子分析工具同时分解两个数据矩阵的共享结构。选择的依据取决于你想回答的问题——如果要定义细胞类型WNN 更合适如果想看调控共变模块MOFA 更合适。WNN 在 Seurat 里的用法很直接前提是细胞 barcode 能一一对应library(Seurat) # 前提rna_obj 和 atac_obj 使用同一批细胞 barcode rna_obj - RunPCA(rna_obj, npcs 30) atac_obj - RunTFIDF(atac_obj) atac_obj - FindTopFeatures(atac_obj, min.cutoff q75) atac_obj - RunSVD(atac_obj) # 设置两个 assay 后直接合并为 multi-assay 对象 obj - merge(rna_obj, y atac_obj) obj - FindMultiModalNeighbors(obj, reduction.list list(pca, lsi), dims.list list(1:30, 2:30)) obj - RunUMAP(obj, nn.name weighted.nn, reduction.name wnn.umap, reduction.key wnnUMAP_) obj - FindClusters(obj, graph.name wsnn, algorithm 3, resolution 0.8)整合之后我会重点检查两种数据在 WNN 空间里的一致性是否存在明显偏置——比如某些细胞类型只靠 RNA 数据驱动、另一些只靠 ATAC 数据驱动。出现这种情况时通常意味着其中一组数据对特定细胞类型的捕获效率有问题。WNN 的优势在于它会根据每个细胞的数据质量自动调整权重但如果你发现某个大类完全由 ATAC 主导而 RNA 贡献极低需要回到原始数据检查而不是直接接受整合结果。5. 蛋白质组把转录层面的推测变成蛋白层面的证据5.1 蛋白定量矩阵的前期处理蛋白质组学数据通常来自 DIA/DDA 质谱输出是一张“蛋白 × 样本”的定量矩阵常见列名是Protein.ID、Gene.Symbol以及各样本的强度值。拿到矩阵后我一般不急着做差异分析而是先做三件事缺失值过滤、数据标准化、批次效应检查。如果某个蛋白在超过 50% 的样本中缺失直接过滤对于剩下的缺失值如果是 DIA 数据可以查看原始定量是否真的低于检测限而不是被错误过滤处理完缺失值后再做 log2 转换让数据分布接近正态最后用中位数标准化或分位数标准化消除样本间的总强度差异。遇到明显批次效应的数据我还会用 limma 的removeBatchEffect()加一重校正。library(tidyverse) # 读取蛋白定量矩阵行为蛋白列为样本 prot - read.csv(protein_matrix.csv, row.names 1) # 缺失值过滤保留在至少 50% 样本中检出的蛋白 keep - rowSums(is.na(prot)) ncol(prot) * 0.5 prot - prot[keep, ] # log2 转换 prot_log - log2(prot 1) # 中位数标准化 prot_norm - sweep(prot_log, 2, apply(prot_log, 2, median, na.rm TRUE), -)5.2 差异蛋白、热图和功能富集可视化蛋白质组的差异分析我常用 limma 或普通的 t 检验。limma 的优势是它能借用所有蛋白的方差信息来稳定小样本下的方差估计对样本量较小的队列更友好。差异分析结果出来后用热图展示关键差异蛋白的表达模式再用 clusterProfiler 做 GO/KEGG 富集看这些蛋白富集在哪些生物学通路。library(limma) library(pheatmap) library(clusterProfiler) library(org.Hs.eg.db) # design 矩阵假设有两组样本group 是分组向量 design - model.matrix(~ group, data meta) fit - lmFit(prot_norm, design) fit - eBayes(fit) # 提取差异蛋白 res - topTable(fit, coef 2, number Inf, sort.by P) deg_proteins - res[res$adj.P.Val 0.05 abs(res$logFC) 1, ] # 关键差异蛋白热图 pheatmap( prot_norm[rownames(deg_proteins)[1:50], ], scale row, annotation_col meta %% select(group), show_rownames TRUE, cluster_cols TRUE ) # GO 富集分析 ego - enrichGO( gene rownames(deg_proteins), OrgDb org.Hs.eg.db, keyType SYMBOL, ont BP, pAdjustMethod BH, qvalueCutoff 0.05 ) dotplot(ego, showCategory 15)蛋白质组数据在整合里的作用是给转录组结果提供一个“执行层面”的验证。mRNA 水平和蛋白水平并不总是线性相关这个过程里存在转录后调控、翻译效率差异、蛋白降解速率等因素。如果转录组显示某个通路激活而蛋白组完全没有相应变化你需要意识到这里面可能存在重要的调控层。反过来如果转录组没有变化但蛋白组变化显著这类蛋白往往是翻译调控或蛋白稳定性层面的候选是一类容易被单组学忽略的信号。6. 空间转录组把基因表达钉回组织坐标6.1 Seurat 处理空间转录组的标准流程空间转录组最常用的数据是 10X Visium它的输出包括一个表达矩阵和对应的组织图像、spot 坐标。分析时可以直接用 Seurat 的 Spatial 流程本质上仍然是标准 scRNA-seq 流程只是多了一个空间坐标信息和图像信息。空间数据的质控要额外注意 spot 是否落在组织区域、总 UMI 数和基因数是否过低以及线粒体比例是否异常升高通常与组织坏死区域相关。library(Seurat) # 读入空间转录组数据 obj - Load10X_Spatial( data.dir filtered_feature_bc_matrix/, filename filtered_feature_bc_matrix.h5, assay Spatial, image spatial ) # 标准化使用 SCTransform效果比 NormalizeData 更稳定 obj - SCTransform(obj, assay Spatial, verbose FALSE) # 降维和聚类 obj - RunPCA(obj, assay SCT, verbose FALSE) obj - FindNeighbors(obj, reduction pca, dims 1:30) obj - FindClusters(obj, resolution 0.6) obj - RunUMAP(obj, reduction pca, dims 1:30) # 在组织图像上看聚类和基因表达 SpatialDimPlot(obj, label TRUE, label.size 3) SpatialFeaturePlot(obj, features c(EPCAM, CD3D, LYZ))6.2 空间可变基因和区域 marker空间转录组区别于单细胞的最大价值在于它的空间位置信息你可以看到特定基因是在组织某个边界、某个结构内、还是在弥散分布。寻找空间可变基因是一个关键步骤Seurat 提供了FindSpatiallyVariableFeatures()其中markvariogram模式用变异函数来识别具有空间自相关特征的基因比简单比较 spot 间差异要更合理。# 找空间可变基因模式选 markvariogram 更稳健 obj - FindSpatiallyVariableFeatures( obj, assay Spatial, selection.method markvariogram, nfeatures 50 ) top_spatial - head(SpatiallyVariableFeatures(obj, selection.method markvariogram), 6) SpatialFeaturePlot( obj, features top_spatial, ncol 3, alpha c(0.1, 1) )用空间可变基因做后续分析时我一般会额外加一步把这些基因映射到单细胞数据里已经注释的 cell type marker 上看哪些细胞类型的 marker 具有强空间结构。这样做的原因很实际——空间转录组的 spot 通常是多个细胞的混合信号纯从 spot 表达定义出来的“结构 markder”可能来自多种细胞类型只有和单细胞的细胞类型 marker 对应起来才能解释这个空间结构背后的细胞学本质。空间转录组在整合中还有一个不可替代的用途它是唯一能把“细胞类型”和“组织微环境区域”关联起来的组学层次。7. 空间代谢组以组织为尺度的代谢物成像7.1 质谱成像数据形态和 Python 处理空间代谢组对很多测序出身的人来说比较陌生但它其实和空间转录组非常互补。主流技术包括 DESI-MSI 和 MALDI-MSI它们通过质谱对组织切片上的每个像素进行扫描得到每个像素位点上一个完整的质谱。数据导出的形式通常是 imzML 格式或者被厂商软件转换为“像素 × m/z 特征强度”的二维矩阵。处理空间代谢组数据的常用工具包括 SCiLS Lab、METASPACE以及 Python 生态里的 pyimzML 和 scikit-image。我讲一下最通用的 Python 处理方式以 imzML 文件为例import numpy as np import matplotlib.pyplot as plt import pyimzML.ImzMLParser as ImzMLParser # 读取 imzML 文件 p ImzMLParser(dataset.imzML) mzs, intensities p.getspectrum(0) # 将逐像素谱转换为空间图像 # 假设已知图像的宽度和高度可以从 imzML 元数据中读取 width, height 100, 80 mz_list np.array([mz for mz, _ in [p.getspectrum(i) for i in range(len(p))]]) # 提取目标 m/z 的强度图m/z 容差通常设为 0.01 Da target_mz 524.36 tol 0.01 ion_image np.zeros((height, width)) for x, y, (mz_array, int_array) in p: mask (mz_array target_mz - tol) (mz_array target_mz tol) if mask.sum() 0: ion_image[y-1, x-1] int_array[mask].sum() plt.imshow(ion_image, cmapviridis) plt.colorbar(labelintensity) plt.title(fm/z {target_mz}) plt.axis(off) plt.show()这段代码的价值在于帮你理解空间代谢组的数据结构——每个像素是一整条质谱你可以把它类比成空间转录组里的每个 spotspot 里的维度是基因这里的维度是 m/z 特征。质谱成像分辨率通常在几十到一百微米左右和 Visium 的 spot 间距在同一量级这让两种数据在组织切片层面的对齐成为可能。7.2 代谢物空间图和转录组联动的思路拿到代谢物空间图之后最有价值的分析不是单独看某个代谢物的绝对分布而是把它和空间转录组做共定位分析。实际操作中先把空间转录组的 spot 坐标和质谱成像的像素坐标统一到同一张组织图像坐标系里然后对每个空间转录组 spot取其覆盖范围内的代谢物像素均值生成一个“spot × 代谢物强度”的关联矩阵最后计算代谢物强度和空间表达基因之间的相关性。from scipy.stats import spearmanr # gene_expression: spot x gene 矩阵来自空间转录组 # metabolite_level: spot x metabolite 矩阵由质谱成像重采样得到 # 计算每个基因和每个代谢物之间的 Spearman 相关 corr_matrix np.zeros((gene_expression.shape[1], metabolite_level.shape[1])) pval_matrix np.zeros_like(corr_matrix) for i in range(gene_expression.shape[1]): for j in range(metabolite_level.shape[1]): rho, pval spearmanr(gene_expression[:, i], metabolite_level[:, j]) corr_matrix[i, j] rho pval_matrix[i, j] pval # 取显著相关的基因-代谢物对画热图 plt.imshow(corr_matrix, cmapRdBu_r, vmin-1, vmax1) plt.colorbar()做这一步时有一个常见的坑空间相关分析容易产生“共线性假象”。在组织里很多基因和代谢物的分布都受同一个组织结构主导比如都富集在上皮层即使二者没有直接的生物学调控关系也会表现出显著相关。我通常加入一个基于空间距离的置换检验把样本空间坐标打乱后重复计算相关性得到一个经验 p 值来校正空间自相关带来的误差。如果置换后显著性基本消失说明观察到的“共定位”可能只是共享空间分布模式而不是真正的代谢-转录关联。8. 跨组学整合从模块拼接走向闭环验证8.1 三个层次的整合框架多组学整合没有唯一的正确方法但有一个比较实用的分析思路分三个尺度做整合。第一个是样本尺度。比如 WES 得到的突变状态某个基因是否突变和蛋白质组中的蛋白丰度或单细胞中的细胞类型比例做相关性分析看突变是否在系统层面带来表型差异。这个尺度适合回答“这个驱动基因到底影响什么”的问题。第二个是细胞尺度。snRNA-seq 和 snATAC-seq 天然共享同一个细胞做 WNN 或 MOFA 联合分析定义出既有转录组特征又有调控特征的细胞状态。这个尺度能回答“我定义的这个细胞类型在调控层面是否真的站得住脚”的问题。第三个是组织尺度。空间转录组和空间代谢组都发生在组织上通过空间坐标对齐后做共定位和相关性分析。这个尺度能回答“基因表达和代谢物分布是否按同一空间模式组织”的问题。三个尺度之间不是独立的而是互相印证WES 找到了某个代谢基因的激活突变那么你在单细胞层面应该能看到这个基因所在的通路在特定细胞类型里上调空间层面应该能看到对应代谢产物在相应区域富集。整条链走通故事就闭环了。8.2 MOFA 多因子分析的代码示例如果需要做更系统的多组学整合我推荐 MOFA。它能同时处理多个数据矩阵提取共享的和各组学特异的因子每个因子相当于多个组学共变的“趋势模块”。MOFA 的使用流程比较清爽library(MOFA2) # 构造多组学数据列表不同组学同一批样本 data_list - list( RNA rna_matrix, # 样本 x 基因 ATAC peak_matrix, # 样本 x peak Protein protein_matrix # 样本 x 蛋白 ) # 创建 MOFA 对象并运行 mofa - create_mofa(data_list) data_opts - get_default_data_options(mofa) model_opts - get_default_model_options(mofa) model_opts$num_factors - 10 mofa - prepare_mofa(mofa, data_options data_opts, model_options model_opts ) mofa - run_mofa(mofa, outfile mofa_model.hdf5) # 查看因子和样本的关系、因子与特征权重的对应 plot_factor(mofa, factor 1, color_by group) plot_weights(mofa, factor 1, nfeatures 20)MOFA 的结果解读需要一点经验。因子排名不一定对应重要性排序要结合每个因子解释了多少方差、和哪个临床特征或细胞状态相关来判断。因子的生物意义需要回到原始的 feature 权重去解读而不能只看因子得分。8.3 整合过程中最常见的坑第一个坑是样本不对齐。同一份组织样本不同组学的样本名可能来自不同命名规则如果不仔细核对批次信息和样本ID后续所有关联分析都会错位。这个问题在单细胞转录组和蛋白组之间特别常见因为前者用采样管编号后者用上机顺序编号看似对应实则错乱。拿到数据的第一步先做一张样本对应表逐项核对。第二个坑是批次效应。不同批次处理的样本其技术差异很容易被误判为生物学差异。在整合前至少要查看每个组学的 PCA 或者 UMAP 是否按批次聚类。如果有明显的批次分离需要通过 Harmony、ComBat 等方法在校正后再做下游整合。第三个坑是过度整合。有些人把所有组学的所有成分全部放进一个模型里跑结果因子非常复杂解释起来无比痛苦。我更推荐的做法是先用最简单的相关性分析建立初步假设再用 MOFA 这类工具做系统性的共变结构分解。分析到最后你会发现真正能支撑结论的往往不是几十个成分而是几个清晰、可重复的模式。我自己这些年每次做多组学项目最深的体会是整合不是靠 fancy 的算法而是靠清晰的实验设计和对每一层数据质量的严格把关。分析流程中多花两小时做元数据核对和质控可以避免在后续关联分析中浪费两周去排查问题。先把每一层数据做成可信的、可解释的结果再讲整合的故事这才是多组学分析真正可靠的工作方式。