多模态数据拿到手,最容易踩的坑是先merge再谈"整合"。这两个动作在 Seurat 里根本不是一回事:merge只是把两个对象的矩阵按列拼到一起,细胞还是各过各的;而跨模态整合要做的是让两个原本互不相识的特征空间产生对应关系,把一种模态上已经积累的生物学结论"翻译"到另一种模态上去。我见过太多项目卡在这一步——单细胞转录组和单细胞开放染色质明明是同一批细胞测的,画出来的 UMAP 却是两条互不相交的轨道,聚类结果各说各话。这篇内容就聊 Seurat 里跨模态整合这条链路:它到底包含哪几个性质完全不同的问题、每种问题该走哪条函数路径、参数背后的取舍逻辑是什么、以及我在实际项目里反复踩到的那些坑。适合已经跑通单模态流程、正准备把两种或多种模态拼到一起的分析者,也适合那些跑完了整合但不确定结果可不可信的人。
1. 先把"跨模态整合"拆开:Seurat 里其实是三件性质不同的事
很多人把"整合"当成一个按钮,按下去就期待两个模态自动对齐。实际在 Seurat 的语境里,这个词至少覆盖三个独立问题,它们的目标、输入形态和输出结果都不一样。分不清这三件事,后面的参数调整基本等于瞎调。
1.1 第一件事:同模态内部的批次效应消除
这是最常被误称为"跨模态整合"的场景。你手上有三批转录组数据,来自不同时间点、不同操作者,同一群细胞在 UMAP 上按批次分成了三坨。这时候要做的是批次校正,输入是同一个模态(都是基因表达),输出是一个消除了技术差异的共享嵌入。CCA、RPCA、Harmony、scVI 这些都是干这个的,跟"跨模态"没有关系。
我之所以强调要先把这件事单独拎出来,是因为它经常和后面两件事混在一起做,导致问题定位困难。正确的顺序是:每个模态内部先各自完成批次校正,再考虑模态之间的对齐。反过来做——先把不同模态拼起来再做批次校正——你会发现校正算法把"模态差异"也当成批次差异给抹掉了,结果就是两种细胞类型被强行捏成一团,看起来漂亮但完全错。
1.2 第二件事:跨模态的"翻译",也就是标签迁移
真正的跨模态问题之一:我有一个模态标注得很清楚(比如有蛋白标记、有明确的门控结果、或者有大量已发表注释),另一个模态只有裸矩阵,没有标签。我要把标签从前者搬到后者,或者把后者的细胞投影到前者的坐标系里。
这个任务的关键在于必须存在一个"双向都测过"的中间层。桥接整合(bridge integration)就是为这件事设计的:它要求你提供一个较小的、同时具备两种模态的数据集当字典,用这个字典建立两套特征空间之间的映射关系,再把大数据集映射过去。没有这个中间层,跨模态映射在原理上就无从谈起——两个特征空间之间没有任何共享锚点。
1.3 第三件事:多模态联合嵌入,同一批细胞的两个视角互相投票
同一批细胞既测了转录组又测了染色质开放性(或者蛋白丰度),两种测量描述的是同一群细胞的不同侧面。这时候不需要"翻译",需要的是融合:让每个细胞的最终邻居关系由两个模态共同决定,而不是简单取交集或取并集。加权最近邻(WNN)解决的就是这个问题。
它的价值在于鲁棒性:某个细胞在转录组上看不出明显特征(比如表达谱很稀疏),但在染色质开放性上有清晰的调控信号,WNN 会自动把权重倾斜到信息量更大的那个模态上。反过来,如果某个细胞的开放染色质数据质量很差(peaks 数极低),权重就会偏向转录组。这种逐细胞的自适应是简单拼接降维做不到的。
1.4 三条路径的对照与选型依据
把上面三件事整理成一张表,实际项目里对着数据形态查一下,基本不会走错路:
| 你手上的数据形态 | 想达成的目标 | 该走的路径 | 主要输出 |
|---|---|---|---|
| 单一模态,多个批次 | 消除批次差异 | IntegrateLayers+ CCA/RPCA/Harmony | 新的降维结果 |
| 模态 A 有标签 + 模态 B 无标签 + 一个双向都测的小数据集 | 标签迁移 / 投影 | CreateBridgeReferenceSet+BridgeIntegration | 映射后的嵌入与预测标签 |
| 同一批细胞同时有两种模态 | 联合聚类与降维 | FindMultiModalNeighbors+ WNN | weighted.nn图与联合聚类 |
| 两种模态测在不同的细胞上,也没有桥接集 | 无法直接做 | 必须先补桥接数据或改用伪配对策略 | 无 |
表格最后一行是我特意加的。有些项目在立项阶段就没规划桥接样本,后期想补只能重新实验,代价极高。如果你正在设计实验,哪怕多花一点成本,也建议留一小批细胞做双模态共测——这批数据的价值会在整个项目周期里反复兑现。
1.5 为什么 v5 要把这套东西重构成 Assay5
Seurat v5 引入Assay5和layers的概念,不是纯粹的重构洁癖,而是跨模态场景逼出来的。老版本的Assay对象里,一个基因矩阵只能存一份 counts,多批次或多模态要么拆成多个对象放在 list 里,要么用split挂在一起但操作起来很别扭。桥接整合需要在同一个对象里同时持有"参考嵌入"和"查询嵌入",还要能按数据集切片计算,老结构撑不住。
Assay5把矩阵按 layer 组织:counts.1、counts.2、data.1……这种命名方式让IntegrateLayers可以只对指定 layer 做计算,而不动其他数据。代价是很多老代码会失效——object@assays$RNA@counts这种直接访问槽位的写法在 v5 里拿不到东西,得改成LayerData(object, assay = "RNA", layer = "counts")。这个改动看起来琐碎,但它是后面所有跨模态操作能跑通的前提,建议在动手做整合之前,先花半小时把项目里所有直接访问槽位的代码清理干净。
2. 塞进同一个对象之前的四道检查
整合跑出来的结果诡异,八成问题出在数据准备阶段,而不是整合算法本身。下面这四项检查我现在的习惯是写成脚本固定执行,每次新数据进来先跑一遍,能省掉大量返工。
2.1 barcode 对齐:字符串拼接带来的隐性错位
两种模态分开建库、分开测序,barcode 天然会重复。常见的处理是给每个模态的 barcode 加后缀,比如AAACCTGAGAAACCAT-1变成AAACCTGAGAAACCAT-1_rna和AAACCTGAGAAACCAT-1_atac。这一步本身没问题,问题出在后缀风格不统一:有的流程用_1、_2,有的用-rna、-atac,还有的直接在合并时重命名成cell_00001这种序号。
后果很隐蔽:FindMultiModalNeighbors不报错,但它按细胞名去匹配两个模态的降维矩阵,名字对不上就会静默地用错行——你看到的联合聚类结果是两个模态的随机配对。我遇到过一次,UMAP 上出现了一条奇怪的"中间带",查了两天才发现是后缀不一致导致部分细胞错配。
检查方法:
head(colnames(obj[["RNA"]]))和head(colnames(obj[["ATAC"]]))逐字符对比,再用length(intersect(a, b))确认交集数量等于预期细胞数。交集数量对不上就不要往下走。
2.2 特征命名的三套体系:基因符号、Ensembl ID、peak 坐标
转录组这边,同一个基因可能有GENE1、Gene1、gene1三种写法,取决于参考基因组版本和注释文件来源;也可能一部分是基因符号、一部分是 Ensembl ID。开放染色质这边,peak 坐标的格式差异更大:chr1-1000-2000、chr1:1000-2000、chr1_1000_2000,甚至有的工具输出 0-based 有的输出 1-based。
单模态分析时这些问题往往不致命,因为内部一致就行。但跨模态场景下,如果你的桥接字典和查询数据用了不同的注释版本,映射关系会整体偏移——不是完全错,而是错一点点,检测不出来但结果就是不够干净。
| 不一致类型 | 典型现象 | 处理方式 |
|---|---|---|
| 基因符号大小写混用 | 部分基因在两个对象里被视为不同特征 | 统一转大写后重建对象 |
| 混用符号与 Ensembl ID | 桥接时锚点数量异常少 | 全部转成同一种 ID 体系 |
| peak 坐标分隔符不同 | 交集为空,桥接直接失败 | 统一成同一格式的字符串 |
| 参考基因组版本不同 | 映射结果系统性偏移 | 必须回到原始数据用同一版本重新定量 |
最后一行是最麻烦的:如果两个模态用了不同版本的参考基因组,靠后期字符串处理是修不回来的,只能重跑定量流程。做多模态项目时,把所有模态的参考版本统一写进项目文档,是我认为性价比最高的一条规范。
2.3 Assay5 的分层语义:counts、data、scale.data 谁该是谁
跨模态整合对数据处于哪一层非常敏感。以开放染色质为例,进入 LSI 降维之前必须做 TF-IDF 变换,变换后的结果应该落到data层,而原始的 peak 计数留在counts层。如果你把 TF-IDF 的结果写回了counts,后续任何依赖原始计数的操作(比如差异可及性分析)都会得到荒谬的结果。
转录组这边同理:NormalizeData之后的对数归一化值在data,ScaleData之后的在scale.data。做整合时,IntegrateLayers默认基于data层计算锚点,如果你只跑了ScaleData而没跑NormalizeData,锚点质量会明显下降。
我现在的习惯是在整合前打印一份数据状态清单:
check_layers <- function(obj, assays) { for (a in assays) { cat("=== Assay:", a, "===\n") cat("cells:", ncol(obj[[a]]), " features:", nrow(obj[[a]]), "\n") for (lyr in Layers(obj[[a]])) { m <- LayerData(obj, assay = a, layer = lyr) cat(sprintf(" %-12s dim = %d x %d, range = [%.3f, %.3f], nnz = %.2f%%\n", lyr, nrow(m), ncol(m), min(m), max(m), length(m@x) / (nrow(m) * ncol(m)) * 100)) } } } check_layers(obj, c("RNA", "ATAC"))输出的稀疏度(nnz)尤其值得盯一眼。转录组的data层稀疏度通常在 5%–20%,开放染色质的counts层可能只有 1%–5%。如果某个层稀疏度接近 100%,说明你误把稠密矩阵(比如 scale 之后的结果)当成 counts 用了,这种数据进整合流程基本必然出问题。
2.4 元数据列的类型陷阱
元数据看起来人畜无害,但它在跨模态流程里承担了关键角色:数据集来源、模态标识、桥接标记都要靠它来区分。常见问题有两个。
一是列名冲突。合并两个对象时,如果两边都有叫orig.ident的列,merge会把它们拼在一起或者覆盖,结果就是你无法区分哪个细胞来自哪个模态。解决办法是合并前先把关键列重命名成带模态前缀的名字,比如rna_batch、atac_batch。
二是类型不一致。数值型的批次标签和字符型的混用,会导致IntegrateLayers里的分组逻辑失效——它可能把1和"1"当成两个不同的组。这种事不会报错,只是结果会莫名其妙地多出一组。用sapply(obj@meta.data, class)快速扫一遍,把所有分类变量统一成factor或character,别偷懒。
3. 桥接整合:用"双语细胞"当字典的完整落地链路
桥接整合是 Seurat v5 里我认为最有价值也最容易被低估的功能。它的思路很像翻译:你不懂两种语言之间的对应关系,但手上有一本双语词典(那个同时测了两种模态的小数据集),就能借助词典把一种语言的句子翻成另一种。
3.1 桥接集的选择标准:为什么不能随便挑
桥接集质量直接决定映射质量,但在实际项目里它经常是最将就的一环——可能是实验时顺手多测的几个样本,也可能是从公开数据里扒的。我总结了三条筛选标准,按优先级排。
第一,模态覆盖度。桥接集必须在两种模态上都有足够深度的数据。如果桥接集的开放染色质部分每个细胞只有几百个 peaks,那它提供的映射信号还不如随机猜。我一般要求桥接集在两个模态上的中位特征数不低于主数据集的 60%。
第二,细胞类型覆盖度。桥接集里出现的细胞类型,必须是查询数据集里也存在或至少相关的。如果你的桥接集只有外周血免疫细胞,而查询数据是实体组织的上皮细胞,映射结果会把这些上皮细胞强行塞进免疫细胞类型里,而且不会报错。这一点极其隐蔽,因为映射后的 UMAP 看起来完全正常。
第三,与目标数据的批次接近性。桥接集和查询数据之间的技术差异越小,映射越准。如果桥接集来自完全不同的建库方式或不同的操作时间,映射的噪声会显著上升。理想情况下,桥接样本和主样本应该是同批次制备的。
3.2 构建 bridge reference set 时的两个关键决定
CreateBridgeReferenceSet的参数不算多,但有两个必须想清楚。
第一个是bridge.assay.name。这指的是桥接集里那个"跨模态"的 assay——如果你要桥接的是转录组和蛋白,那桥接集里的蛋白 assay 就是关键;如果是转录组和开放染色质,那通常选开放染色质那个 assay(或者它的活性打分矩阵)。选哪个取决于你想让哪一侧的信息作为映射的"锚"。
我倾向于选数据更稀疏、更难直接比较的那一侧。理由是:信息密集的那一侧(比如转录组)本身就能提供比较稳定的降维,映射的主要不确定性来自稀疏侧,把它作为桥接锚能让映射过程更聚焦在真正需要解决的问题上。
第二个是bridge.reduction。这需要提前在桥接集上算好一个降维结果,通常是联合降维(比如 APCA 或者直接对两种模态的特征拼接后做 PCA)。这个降维结果的质量决定了整个映射的骨架——它一旦歪了,后面所有的映射都会跟着歪。建议在做这一步之前,先用桥接集自己的数据验证一下:在这个降维空间里,已知的细胞类型是否已经分开了?如果分不开,回头调整降维参数,别急着往下走。
# 1) 在双模态桥接集上先算好联合降维 bridge <- FindMultiModalNeighbors( bridge, reduction.list = list("pca", "lsi"), dims.list = list(1:30, 2:30), k.nn = 20, l2.nn = 200 ) bridge <- RunUMAP(bridge, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_") # 2) 构建桥接字典 bridge.set <- CreateBridgeReferenceSet( object = bridge, bridge.reduction = "wnn.umap", # 或 "apca",取决于你的降维方案 bridge.assay.name = "ATAC", # 稀疏侧作为锚 bridge.assay.key = "bridge_atac", dataset.name = "bridge", verbose = TRUE )提示:
bridge.reduction这一项在官方示例里常写成一个专为桥接准备的联合降维名字。不同版本对这个参数的取值要求略有差异,跑之前务必用?CreateBridgeReferenceSet确认当前版本接受什么,别照抄旧教程。
3.3 IntegrateLayers 走 BridgeIntegration 的参数取舍
真正执行映射的是IntegrateLayers的BridgeIntegration方法。这里最需要权衡的是k.anchor和k.filter这类锚点相关参数。
k.anchor控制每个细胞在另一侧找多少个候选锚点。默认值偏保守,在细胞类型差异明显的数据上表现不错;但如果你的查询数据里有稀有细胞类型,默认值可能直接把它们过滤掉——因为稀有细胞在桥接集里找不到足够多的相似邻居,锚点数量不达标。这种失败模式不报错,只是稀有亚群会在映射后的 UMAP 里凭空消失,或者被硬塞进最近的大类。
我的经验做法是做一次参数扫描:k.anchor分别取 5、10、20,各跑一遍,然后统计每种设置下映射后保留的细胞类型数量(用已知 marker 打分快速判断)。如果某个设置下稀有亚群消失,就说明这个值太大,得往小调。
k.filter则控制锚点的可靠性过滤强度,值越大过滤越狠,锚点越可靠但数量越少。数据质量好的时候可以调大,质量差的时候调大反而会让锚点不够,映射变得不稳定。
obj <- IntegrateLayers( object = obj, method = BridgeIntegration, orig.reduction = "pca", new.reduction = "integrated.bridge", bridge.reference = bridge.set, k.anchor = 10, k.filter = 200, verbose = TRUE )跑完之后立刻做一件事:看obj@commands里记录的锚点数量。锚点数太少(比如低于查询细胞的 5%),说明桥接集和查询数据的匹配度不够,后面的结果基本不用看了,老老实实回去检查数据准备工作。
3.4 映射完先别急着出图:三个质量指标
映射结果的可信度检查不能靠"图看着顺眼"。我给自己的固定流程是三项。
第一项,映射前后细胞类型的分布一致性。如果你的查询数据集来自几个不同的样本或条件,映射之后每个条件下的细胞类型比例应该大体稳定。如果某个条件下突然冒出一类在原始分析里完全没见过的细胞,那多半是映射错了,不是发现新类型。
第二项,标签预测的置信度分布。桥接整合通常会给出预测标签和对应的分数,看这个分数的分布形态:如果大部分细胞集中在高位、只有少数低分,说明映射整体可靠;如果分数呈双峰或者整体偏低,说明桥接集覆盖度不足。
第三项,也是最容易被忽略的一项:反向映射的合理性。从查询数据里挑几个映射置信度最高的细胞,看它们在桥接集里最近邻的细胞类型是否和预测标签一致。这一步等于人工复核,我一般抽 20 个细胞看一眼,几次下来就能积累出对这套流程可靠性的直觉。
4. WNN 加权最近邻:两个模态互相投票的机制与调参
WNN 解决的是另一类问题:同一批细胞被测了两次,两个视角该怎么融合。它不属于"翻译",属于"共同决策"。
4.1 权重是算出来的,不是调出来的
先纠正一个常见误解:WNN 里的模态权重不是用户设定的超参数,而是算法逐细胞计算的结果。这个区别很重要,因为它意味着你没法"手动调高转录组的权重",只能通过调整近邻结构间接影响权重分布。
计算的思路大致是这样:对每个细胞,先在每个模态里各找 k 个近邻(k.nn控制)。然后对每个细胞和每个模态,计算该细胞到其近邻里其他细胞的平均距离,这个平均距离作为该细胞在该模态下的"带宽"。带宽的作用类似于一个自适应尺度——稀疏数据带宽大,稠密数据带宽小。
接下来是关键一步:用模态 A 的邻居集合去模态 B 的距离空间里衡量距离,得到跨模态亲和度。亲和度用的是指数衰减形式,距离越远衰减越快,带宽控制衰减速度。最后每个细胞的模态权重,是两个方向亲和度归一化的结果,权重和为 1。
这个机制解释了一个现象:为什么 WNN 在稀疏数据上特别有价值。转录组数据里那些只测到几十个基因的细胞,其转录组带宽会非常大,跨模态亲和度相应很低,权重自然被压到很低的位置,最终这些细胞的邻居关系几乎完全由开放染色质决定。反过来也成立。
4.2 两个近邻参数的实际影响
k.nn和l2.nn这两个参数,名字看起来差不多,作用完全不同。
k.nn是每个模态内部找近邻的数量,默认 20。这个值影响的是模态内部结构的精细程度:调小会让局部结构更敏感,适合细胞类型边界清晰的数据;调大会让结构更平滑,适合噪声大或者细胞类型连续过渡的场景。我的习惯是先按默认值跑一遍,然后取k.nn为 10 和 30 各跑一遍,看聚类结果的稳定性——如果三种设置下大类划分一致,只是亚群有差异,说明数据结构本身是清楚的,按需选细粒度;如果三种设置下大类都不一样,那说明数据本身有问题,先别调参数,回去查质量。
l2.nn是计算权重时用的第二层近邻范围,默认 200。它不直接影响聚类,而是影响权重计算的平滑程度。调大这个值会让权重变化更平缓,代价是丧失一部分逐细胞的适应性。这个值在数据规模差异很大时需要相应调整——如果你的数据只有几百个细胞,200 这个默认值就显得太大了,可以考虑缩小到细胞数的 10%–20%。
obj <- FindMultiModalNeighbors( obj, reduction.list = list("pca", "lsi"), dims.list = list(1:30, 2:30), k.nn = 20, l2.nn = 200, prune.SNN = 1/15 ) obj <- FindClusters(obj, graph.name = "wsnn", algorithm = 3, resolution = 0.5) obj <- RunUMAP(obj, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_")注意FindClusters在 v5 里必须显式指定graph.name,不写会落到默认的 RNA_snn 图上,得到的结果跟 WNN 完全无关——这个坑我踩过,聚完类发现和单模态结果一模一样,还以为是 WNN 没起作用。
4.3 用权重分布做诊断:哪些细胞在"靠一个模态撑场面"
权重本身是极好的诊断工具。取出来画一张按细胞类型的权重分布图,信息量很大。
如果某个细胞类型在两个模态上的权重接近 50/50,说明两种测量对该类型都有清晰信号,这类结果通常最可靠。如果某个细胞类型在转录组上的权重长期低于 20%,说明它的转录组信号很弱(可能是低 RNA 含量的细胞,或者测序深度不够),这时候该类型的聚类结果要谨慎对待——它几乎完全由开放染色质决定,任何染色质侧的偏倚都会直接反映到聚类上。
# 较新版本可直接取权重表;老版本可从 weighted.nn 图中反推 w <- ModalityWeights(obj, reduction = "weighted.nn") colnames(w) <- c("rna_weight", "atac_weight") obj$rna_weight <- w$rna_weight # 按细胞类型汇总,看是否存在严重倾斜 agg <- aggregate(cbind(rna_weight, atac_weight) ~ seurat_clusters, data = obj@meta.data, FUN = mean) print(agg[order(agg$rna_weight), ])如果输出里出现某个簇的rna_weight平均值低于 0.15,我通常会在后续分析里给这个簇单独标注,任何基于转录组的差异表达结论都要先存疑。
4.4 WNN 聚类和单模态结果打架时的处理顺序
大概率你会遇到这种情况:WNN 分出来 12 个簇,转录组单独分出来 10 个,开放染色质单独分出来 15 个,三套结果对不上。这时候的处理顺序很关键。
第一步,先看差异簇的权重构成。如果某个簇在 WNN 里被合并了,但两个单模态都把它分开,那说明这个簇在两个模态里的信号方向一致,是 WNN 的l2.nn平滑过度导致的,考虑调小k.nn。
第二步,看是否是分辨率问题。WNN 的图结构和单模态不同,同样的 resolution 在不同图上对应的簇数不一样。别用同一个 resolution 横向比较,而是以"簇数"为锚点,把 WNN 的 resolution 调到和单模态簇数相当,再比较簇的组成。
第三步,如果前两步都排除了,那就是真实的生物学差异——某个细胞群在转录组上和邻居很像,但染色质状态明显不同(比如同一类型细胞的不同激活状态)。这种时候不要强行让三套结果统一,WNN 分得更细恰恰是它的价值所在。我个人的做法是把这类簇在命名时保留层级信息,比如标成T_cell_resting和T_cell_activated,而不是硬塞成一个T_cell。
5. 结果不信怎么办:一套可执行的交叉验证清单
跨模态整合最大的风险不是跑不出来,而是跑出来了但结果是错的,而且错得很像对的。下面这套检查我基本每次都过一遍,覆盖了绝大多数失败模式。
5.1 marker 层面的正向验证
最直接的办法:拿一批你确信无疑的 marker,在两个模态上分别看它们的分布,再在联合嵌入上看它们是否落在预期的位置上。
转录组侧的 marker 好办,常规的差异表达就能做。开放染色质侧需要多做一步——把基因附近的 peak 集合起来算一个可及性打分,再看这个打分是否与该基因的表达在细胞水平上正相关。这个相关性如果很低(比如低于 0.1),说明两个模态的细胞对齐可能有问题,回头查 barcode。
联合嵌入上的验证,我一般挑 5–8 个跨类型差异明显的 marker,在 WNN 的 UMAP 上标出表达量,看它们是否集中在预期的簇里。如果某个 marker 散落在多个相距很远的簇上,而这两个簇在单模态分析里本来就是一回事,那就是整合把不该拆的东西拆了。
5.2 标签转移置信度的正确读法
桥接整合给出的预测分数,很多人直接拿它当"准确率"看,这是个误解。这个分数反映的是"该细胞在桥接字典里的近邻有多集中",而不是"这个标签有多正确"。
具体来说,一个细胞如果有 100 个近邻全部集中在同一种细胞类型,分数会很高;如果近邻均匀分布在三种类型里,分数会很低。但低分不一定是错的——处于分化过渡态的细胞,本来就应该在多种类型之间。真正需要警惕的是"高分但错":某群细胞拿了很高的分数,但你看它们的 marker 明显属于另一个类型。
这种情况通常意味着桥接集本身缺少对应类型,算法只能把细胞塞给最近的那个。所以看分数的同时,一定要结合 marker 做独立判断,两者冲突时以 marker 为准。
5.3 四种典型失败模式与对应现象
| 失败模式 | 表面现象 | 根因 | 排查手段 |
|---|---|---|---|
| 模态不对齐 | UMAP 上两个模态分层或呈两条平行带 | barcode 或特征名不一致 | 交集数量检查、特征名逐条比对 |
| 桥接集覆盖不足 | 稀有细胞类型在映射后消失 | 桥接集里没有对应类型 | 映射前后类型数量对比 |
| 模态主导 | 结果几乎等同于单模态聚类 | 某一模态权重普遍接近 1 | 权重分布直方图 |
| 过度平滑 | 该分开的亚群被合并 | k.nn或l2.nn过大 | 参数扫描对比簇数 |
这四种里,第二种最难发现,因为它的表现是"少了一类细胞",而人往往会下意识地认为那类细胞本来就不存在。养成习惯:映射前后各存一份细胞类型计数表,做差。
5.4 一个容易被忽略的检查:模态间的相关性随时间漂移
如果你的数据是分批采集的,值得做一件事:按批次分别计算两个模态在联合嵌入上的局部一致性,看有没有随时间漂移。
具体做法是:对每个批次,随机抽若干细胞,对每个细胞找它在联合嵌入里的近邻,再分别算这些近邻在转录组和开放染色质空间里的平均距离。如果某个批次的染色质侧平均距离明显偏高,说明这批数据的染色质质量下降,它的聚类结果里染色质的贡献被系统性削弱了。这类漂移在合并分析里完全看不出来,但会影响对批次间生物学差异的判断。
6. 规模上去之后的工程取舍
前面聊的都是方法论,最后说说实操层面的现实问题。跨模态整合的计算开销远高于单模态,一个几万细胞的双模态数据集,全套流程跑下来内存和时间都可能超出预期。
6.1 稀疏矩阵在 Assay5 里的真实开销
直觉上,稀疏矩阵的内存应该和存储的非零元素数量成正比。但实际用起来不是这样:Assay5在多个 layer 并存的情况下,会有相当一部分内存消耗在 layer 的索引结构和对象的元信息上。
我实测过一个十万细胞左右的双模态数据集:转录组counts层非零元素约 3 亿,占用几个 GB;开放染色质counts层非零元素少一半,但内存占用并没有少一半,因为 peak 数量多、列索引开销大。再叠加 TF-IDF 变换后的data层,总占用就上去了。
我的做法是及时清理中间层。比如 TF-IDF 算完之后,如果不打算重跑,可以只在需要的时候保留data层,把临时的中间矩阵用rm()加gc()清掉。另外DietSeurat在跨模态场景下要慎用,它默认会删掉非主 assay 的一些信息,容易把刚要用的层删掉。
6.2 先降维还是先整合
这是个纯粹的效率问题。两种顺序在数学上不等价,结果会有差异。
先降维后整合:每个模态各自算降维(转录组走 PCA,开放染色质走 LSI),然后用高维的降维结果去计算近邻和权重。优点是计算量小,高维降维天然起到降噪作用;缺点是降维过程中可能丢掉模态间的一致信息。
先整合后降维:在原始特征空间上做对齐,然后再统一降维。理论上更"完整",但计算量巨大,而且稀疏数据在原始空间上的对齐效果未必比降维后更好。
我的实测结论是:对于细胞数超过五万的数据集,先降维后整合的效果和效率综合更好,而且差异通常在选择最优降维数(dims.list)的过程中被吸收掉了。只有在数据量小、且你怀疑某个模态的关键信号被降维抹掉时,才值得尝试先整合的路线。
6.3 中间产物落盘策略
跨模态流程的每一步都很贵,落盘策略直接影响迭代效率。
务必把桥接字典对象单独存一份。它和主数据无关,可以跨项目复用,重新算一遍的成本往往是最高的。其次是把每个模态的降维结果存下来——参数调整时最常改的就是整合和聚类这几步,降维通常不动。
最后是weighted.nn图的存储。这个图对象体积不小,但重算成本高,而且如果你要做参数敏感性分析,存下来能省很多时间。用saveRDS压缩存,别用默认设置。
6.4 一个我觉得值得的取舍:宁可慢一点,也别跳过诊断
跨模态整合有个特点:错误的结果和正确的结果长得一样漂亮。单模态分析里,聚类结果诡异你一眼就能看出来;跨模态场景下,权重失衡、桥接覆盖不足这些问题都会产生看起来很正常的 UMAP。
所以我在项目里给自己定了条规矩:任何一次跨模态整合,必须先出三张图再往下走——权重分布直方图、映射前后类型计数对比、marker 在联合嵌入上的分布。这三张图花不了十分钟,但能拦住大部分后知后觉的返工。踩过几次坑之后我越来越确信,这部分"额外"的工作才是整个流程里最省时间的一环——真正昂贵的从来不是算力,而是在错误结果上继续做两周的下游分析。
关于参数扫描,我最后分享一个小技巧:把k.nn、l2.nn、k.anchor这几个关键参数写成一个列表,用循环批量跑,每次把簇数、权重均值、稀有类型保留情况记进一张表。跑完之后这张表本身就是判断数据质量的依据——如果参数在合理范围内变化时结果剧烈波动,说明数据本身不够扎实,先回去查数据,别指望调参数能救。