GLUE多组学单细胞整合：从特征链接图到调控推断的完整指南

我第一次拿到 10x Multiome 的数据时想得很简单把同一批细胞的 RNA 表达矩阵和染色质可及性矩阵拼在一起然后跑一个常规聚类。结果显而易见地翻车了——两个矩阵的特征维度差了将近两个数量级稀疏度完全不同单纯拼接之后除了少数高表达基因整个低维表示几乎被 ATAC 的千万个 peak 淹没转录组信息完全没法体现。后来读到高歌老师团队提出的 GLUEGraph-Linked Unified Embedding核心就是那篇《Multi-omics single-cell data integration and regulatory inference with graph-linked embedding》我才意识到多组学整合的关键根本不在于“拼数据”而在于“找锚点”。GLUE 解决的核心问题是用一个先验的跨组学特征关系图谱把原本没有共同坐标的不同组学数据映射到同一个语义空间。它不但能做多组学单细胞数据整合还能在整合之后顺势做调控推断——比如基因和远端 enhancer、基因启动子和甲基化位点之间的调控关系推断。这套思路对正在做 multiome、CITE-seq、sci-CAR 或者任何“多组学 单细胞”场景的人都有参考价值哪怕你只是想把两个批次的数据对齐也能从中理解“基于图约束的嵌入”到底比暴力拼接强在哪。1. 先别急着拼接矩阵多组学数据整合卡在哪一步多组学整合听起来像是一个存储问题实际上是一个生物语义对齐问题。不同组学的特征不是一个单位的RNA 的特征是基因ATAC 的特征是基因组上的开放区域甲基化的特征是 CpG 位点。基因表达量和染色质可及性虽然都反映细胞状态但它们之间没有天然的“同维度对应关系”。你不可能把基因 A 的表达量直接和 peak B 的可及性做减法也不可能把二者当成同一组特征的两次重复测量。1.1 不同组学的“特征”从来就不对等基因表达矩阵大约两万多个基因ATAC 矩阵少则数万、多则上百万个 peak甲基化矩阵又是几百万个 CpG 位点。即便做了特征筛选降到一个量级也仍然存在根本差异基因是人转录本的单位peak 是基因组区间的单位CpG 是单个甲基化位点的单位。它们背后虽然共享同一套基因组坐标但坐标本身并不等于“同一个东西”。更麻烦的是数据分布完全不同。RNA 表达适合负数二项分布建模ATAC 的 peak 可及性通常被处理成二值或计数DNA 甲基化则是 beta 值。如果你强行把它们放到一个矩阵里等于逼着同一个模型同时解释三种完全不同的统计机制结果往往是低方差、高维度的特征主导全局结构而那些真正有调控意义的弱信号全被淹没了。1.2 细胞对齐与特征对齐其实是同一件事很多人以为多组学整合的重点是让“同一类细胞”在不同组学数据里聚到一起所以注意力全部放在细胞 embedding 上。但实际做下来你会发现细胞能不能对齐取决于特征能不能对齐。如果两个组学里的基因表达和染色质可及性在共享空间里没有对应关系细胞即便被强行拉近也是“假整合”——下游一算 peak-to-gene 关联就露馅因为特征层面的对应关系本身就是乱的。所以真正有效的多组学整合必须同时完成两件事第一把细胞映射到共享低维空间第二让不同组学的生物学相关特征在这个空间里彼此靠近。GLUE 的思路就是先把第二个问题显式建模出来用它作为“路标”去引导第一个问题。这也是 graph-linked embedding 与早期 CCA、MOFA 这类方法的本质区别——它不是在找两个矩阵的最大协方差方向而是假设特征之间的调控连接是已知的、可构造的图结构然后让整个嵌入过程服从这张图。2. GLUE 如何用“锚链图”把不同组学牵到一起graph-linked embedding 原理拆解GRAPH-LINKED EMBEDDING 这个名字值得拆开看。“Linked”指的是预先定义的跨组学特征连接也就是链接图。以转录组和染色质可及性为例我们可以构造基因启动子区域和开放染色质 peak 之间的连接如果一个 peak 落在某个基因启动子附近或通过 loop 结构锚定在某个基因上就给这两个特征加一条边。这张图就是整个模型的“锚链”。2.1 链接图不是从数据里现学出来的而是来自先验知识GLUE 和很多端到端深度学习模型的不同点在于它不指望模型自己发现基因和 peak 之间有什么关系。它要求你在训练前显式提供一张特征链接图比如“基因 ENSG00000141510 的启动子与 chr12:12345-12400 这个 peak 有连接”。这些连接可以来自基因组距离、ENCODE 的 ChIP-seq 数据、Roadmap 的增强子注释、外显子捕获到的 loop 数据等。为什么一定要先验图原因很实际单细胞数据本身的噪声太大如果让模型自己从零学习基因和 peak 的关联看到的更多是测序深度和 dropout 带来的假信号而参考基因组注释和染色质构象数据是经过几十个样本、多种技术反复验证的用它们约束模型可以让模型把精力放在“如何用已有调控知识去解释新数据”上而不是从头猜测生物学关系。2.2 图自编码器把特征映射到共享空间在技术上GLUE 对每个组学都构建一个组学特异性编码器把细胞的特征向量映射到低维共享空间。但真正核心的是对特征 embedding 的约束它把上一步构造的链接图当作一个图网络用一个图自编码器去学习每个组学特征在低维空间里的表示。学过图神经网络的人都知道图自编码器天然适合处理这种“节点是特征、边是调控关系”的图结构。节点嵌入需要保持图上邻近结构如果基因 A 和 peak B 之间存在链接边那么它们的特征 embedding 在共享空间里就应该靠近。这只是第一步。更重要的是细胞 embedding 在训练时会通过某种相似度函数与特征 embedding 对齐。粗略理解就是一个细胞如果在 RNA 组学里高表达基因 A那么它的细胞 embedding 应该靠近“基因 A”的特征 embedding如果同一个细胞在 ATAC 组学里对应 peak B 可及性高而“peak B”的特征 embedding 又恰好离“基因 A”很近那么这个细胞的跨组学表示就被链接图“捏”到了一起。类比一下基因 A 和 peak B 是两份不同语言写的说明书里的同一关键词链接图是双语词典GLUE 做的就是同时学习“英文版说明书”和“中文版说明书”里句子的语义向量要求“关键词处含义一致”。2.3 对抗策略拉平不同组学的分布差异即使有了特征链接图细胞 embedding 在不同组学之间的分布仍然可能有系统偏移。比如同一个细胞类型在 scRNA-seq 里的转录组测序深度普遍比在 sci-CAR 里高如果只看表达水平模型很容易把这些“技术差异”学进 embedding。GLUE 的做法是用一个对抗判别器去判断某个细胞 embedding 来自哪个组学与此同时细胞编码器会拼命让判别器分不出来。这个对抗过程最终让不同组学来源的相同细胞类型在 shared space 中落到同一片区域。这里要强调“特征链接”和“对抗对齐”是互补的链接图管的是特征层面的生物学关系对抗判别器管的是组学层面的技术偏移。只有对抗没有图结果就是单纯的数据分布对齐组学内真正的功能关系丢失只有图没有对抗又可能被技术差异带偏。GLUE 把两者放进同一个目标函数里联合优化所以它做了“整合 调控推断”两件事而不是先整合再单独做调控预测。3. 构建输入feature link 是整条链的灵魂很多人第一次跑 GLUE 时最不习惯的就是它不只吃一个表达矩阵还要你额外提供一个特征链接表。如果这一步没做好后面模型再精巧也白搭。我见过最典型的错误是拿 GRCh38 版本的基因注释去匹配 hg19 版本的 peak 坐标结果是基因和 peak 的链接表基本是乱连训练出来的 embedding 自然经不起下游验证。3.1 数据格式与基本清洗GLUE 基于 scanpy/anndata 生态RNA 和 ATAC 各自是独立的 AnnData 对象。每个 AnnData 的 obs 代表细胞var 代表特征细胞名、特征名必须是统一标准。RNA 的 var_names 建议直接用 Ensembl gene ID 或官方 symbol 都可以但一定要和链接表里的基因 ID 类型保持一致。ATAC 的 var_names 我习惯用“chr1:1000-2000”这种经典 BED 格式表示因为后续和链接表做坐标匹配最方便。清洗方面RNA 要过滤低质量细胞和低表达基因ATAC 也要过滤峰覆盖度过低的细胞。这里没有统一阈值我通常按每个数据集的实际分布取拐点。有一点要特别注意ATAC 数据的稀疏程度远高于 RNA在过滤时千万不要只保留在大多数细胞里都出现的 peak——那样会把大量有意义的远端调控元件丢掉而远端调控元件恰恰是 GLUE 做调控推断最喜欢的东西。3.2 链接表的构造坐标、距离和注释来源链接图的核心就是一张三列甚至两列的表左边是 RNA 特征的 ID右边是 ATAC 特征的 ID中间是权重可选。以 gene-peak 链接为例常见的构造方式是取每个基因的启动子区域一般是转录起始位点 TSS 前后各 2.5kb 或 5kb然后检查哪些 peak 与之有交集。想纳入远端调控关系可以把范围扩到 TSS 上下游 50kb或者引入 loop 数据比如 Hi-C 和 HiChIP 的锚点互作。下面是我常用的构造流程# 伪代码实际操作根据基因组版本调整 gene_bed read_gtf(gencode.v40.annotation.gtf) # 提取gene区间 peak_bed read_bed(peaks.bed) # scATAC peaks links overlap_by_coords( gene_bed[[chr, start, end, gene_id]], peak_bed[[chr, start, end, peak_id]], promoter_flank5000, distal_flank50000 ) # 输出 links.tsv包含 gene_id, peak_id, weight除了 gene-peak 链接GLUE 还支持更复杂的多组学链接比如基因启动子与 CpG 甲基化位点的链接、基因启动子与 TF motif 所在 peak 的链接。链接不是越多越好而是要“有生物学依据”。我宁可只用 ENCODE 和 SCREEN 的高置信注释也不在第一步就把所有距离小于 10kb 的组合全塞进去因为错误链接会像一个坏锚点把无关特征拉在一起对 embedding 的污染是全局性的。3.3 配对不同组学时的特殊处理如果你拿到的是 10x Multiome 这种同一个细胞上同时测了 RNA 和 ATAC 的数据那么细胞 level 是配对的。GLUE 完全支持这种配对数据但它并不依赖配对关系做整合真正决定整合质量的仍然是特征链接图。实际上有个经验配对数据里如果只把两个矩阵拼在一起跑还会引入大量技术噪声而 GLUE 是先分别编码两个组学再通过链接图在 latent space 对齐所以即使配对也应该用“非配对”的心态去设计实验——不要把“同一细胞”当作整合的唯一依据要留出一部分配对信息做验证。4. 训练 GLUE 时我反复调整的几个关键点模型本身是复杂的但真正决定成败的往往是那些容易忽略的训练细节。我在不同数据集上跑 GLUE 时遇到过 latent 空间完全不分离、对抗 loss 归零、训练到一半直接 OOM 等各种情况。下面这些点是我每次新数据集开工前都会确认一遍的。4.1 超参数latent dim 与权重系数latent dimension 的设置不能拍脑袋一般 50 到 100 之间比较常用。太小的话不同生物信号会挤在一起太大的话后续调控推断的 feature embedding 反而很难解释。ATAC 特征数量特别大的时候我会适当提高 latent dim但不会超过 128。另一个关键超参数是链接图重构 loss 与对抗 loss 的权重论文里默认的平衡点做得不错但实际数据噪声更大时权重需要适当调整。判断标准很简单如果整合后不同组学来源的细胞明显分成两团说明对抗权重太低如果所有细胞不分三七二十一都压到一团说明对抗过强把真实的生物学差异也抹掉了。4.2 训练过程监控别只看总 loss我通常同时盯着三个曲线总 loss、链接图重构 loss、对抗判别器 accuracy。判别器 accuracy 如果长时间保持在 0.5 附近说明编码器成功骗过了判别器但如果从第一个 epoch 就维持在 0.5也要警惕判别器可能根本没训练起来这时整合其实没有发生。相反如果判别器 accuracy 一直升到 0.9 以上说明编码器还没跟上来此时要让对抗权重加大或者调整学习率。链接图重构 loss 是我最看重的指标。它如果没降下去说明模型并没有学到基因和 peak 之间的对应关系后面做出来的调控关系就不必看了。我会在每轮结束保存一次 feature embedding随机抽几十对 gene-peak 链接手动看一下它们的余弦相似度是否真的比非链接对更高。这个方法很土但比任何漂亮曲线都直观。4.3 数据规模与显存管理GLUE 在大规模数据集上确实比较吃显存。没有 32G 以上显存的话我建议先用代表性细胞做子采样比如每种注释细胞类型各取 500-1000 个细胞先跑通流程、确认参数合理再放全量数据。另外特征链接图的大小也会影响显存如果 peak 数量超过几十万构建完整邻接矩阵是不现实的需要做 batch 式的图采样。GLUE 本身对大规模图做了优化但你在传数据时一定要检查自己的 feature 数量不要带着两三百万个 peak 直接跑除非你对自己的显存非常有信心。4.4 随机种子与可复现性GLUE 这类基于深度学习的工具对随机种子敏感不同 seed 跑出来的 embedding 在小概率上会有差异尤其是聚类边界上的细胞。我每次正式训练前会固定 seed并且保留最终模型的 checkpoint方便下游反复抽查。做团队项目时我会把数据的版本、链接表版本、seed 和所有超参数全部记录在配置文件里。这看起来是常识但我在实际合作里见过太多“上周跑出来的结果这周复现不了”的尴尬多半都是因为这些细节没记录。5. 整合之后调控推断是怎么“顺水推舟”出来的GLUE 这个名字从广义上说就是“整合 调控推断”。它把特征链接图作为监督信号让基因 embedding 和 peak embedding 在共享空间里靠近。训练完成之后这些 feature embedding 就是天然的调控关系打分器。5.1 从 cosine 相似度到候选调控对给定任意一对基因和 peak我直接取出它们在 latent space 里的 embedding算一个余弦相似度就能得到这个基因-peak 对的调控关联分数。这个分数等价于模型在观察了所有细胞表达和可及性之后对先验链接的一个“重新加权”。先验链接表里可能有一万条边其中不少是基因组距离近但实际不调控的“假阳性”整合后只有那些在真实数据里被持续共表达的 peak 和基因才会获得高相似度。我常用的做法是从链接表里筛出 top 10%-20% 的高分基因-peak 对再叠加 motif 信息看该 peak 里是否富集了对应基因启动子的转录因子结合位点。这样能进一步把直接调控和间接关联分开。拿 BENGI 基准数据做验证时GLUE 的调控推断精度整体上是优于单纯用距离或共开放程度的尤其是在远端 enhancer-gene 关系上。5.2 细胞类型特异的调控网络很多人的需求不是拿到一个全局 gene-peak 列表而是想知道“某个细胞类型里到底谁调控谁”。这时可以把训练好的细胞 embedding 做聚类得到细胞类型注释然后对每一类细胞单独构造一个“细胞子集”重新计算 feature embedding 之间的局部关联。因为细胞 embedding 已经对齐子集里的基因表达和 peak 可及性天然对应之后再做相关性分析会比直接拿原始矩阵算准得多。举个例子我处理过一个 CD4 T 细胞分化的数据全局链接表里某个 FOXP3 上游 enhancer 和启动子的分数并不高但把细胞 embedding 按 Treg 亚群单独拿出来后这个对分数瞬间冲到前 5%。这种细胞类型特异的信号恰恰是单细胞多组学数据最值钱的信息。5.3 注意推断的是“调控可能性”不是因果机制GLUE 的调控推断本质上是“用高维共变信息对先验调控网络重排序”它不能告诉你某个转录因子到底是怎么结合到那个 peak 的也不能告诉我们敲除之后会不会真的改变下游基因表达。如果你想进一步验证需要结合扰动实验、CRISPR screen 或者至少用 eQTL/QTL 或等位基因特异表达数据做外部验证。我见过有人把 GLUE 的一堆高分组直接当成因果结论写进文章这是很危险的审稿人大概率会追问你实验验证。6. 验证整合质量只看 UMAP 是不够的我承认 UMAP 和 t-SNE 是看整合效果最方便的工具但它们也最容易骗人。不同组学的细胞被拉得越近UMAP 看起来越“融合”但融合不一定是好事过度融合会把所有细胞类型都糊成一团。所以我在正式分析里会同时报告两类指标一类衡量批次/组学混合程度另一类衡量生物类别保真度。6.1 组学混合度与生物保守性两个方向都要测组学混合度可以用 kBET、LISI 或 silhouette score 来评估理想情况是来自不同组学但属于同一类型的细胞在空间上充分混合也就是“batch effect”小。但别忘了如果所有细胞都混在一起这类指标反而很高。所以还要同时看生物保守性比如 ASWcell-type silhouette width、ARIAdjusted Rand Index或者直接做无监督聚类的纯度。最稳妥的做法是拿一个独立训练的 scRNA-seq 注释结果去和 GLUE 整合后的聚类比对如果 ARI 掉得厉害说明整合牺牲了真实生物学结构。6.2 我常用的评估清单下面这张表是我每次做完多组学整合后必跑的检查项供你参考评估维度常用指标/方法理想表现注意点组学混合kBET, LISI, 判别器准确率同类型细胞难以区分来源要按细胞类型分层计算生物保守ARI, NMI, cell-type ASW注释标签可分性高不要用过大的聚类数特征对齐链接图中 gene-peak cosine 相似度显著高于非链接对抽样检查避免全量计算下游可迁移label transfer 准确率跨组学预测准确率稳定注意类别不均衡稳定性不同 seed 重复训练结果聚类/标签基本一致至少在主要类型上稳定这个清单看起来复杂但实际跑起来很快。我一般会先跑第 1 和第 2 项如果两者都不差再去抽第 3 项如果第 1、2 项出现明显矛盾比如组学混合度很高但 ARI 很低那大概率是过度整合这时候就需要回头调低对抗权重而不是继续堆 epoch。6.3 别忽略了“翻译一致性”这个特有指标GLUE 还允许你做一件很有意思的验证把某个细胞从 RNA 编码器得到 embedding然后用另一个组学的解码器看它“翻译”过去之后能不能还原出该细胞在 ATAC 里预期的 peak 可及性模式。简单说就是做 cross-omics 重建。如果整合质量好RNA embedding 应该足以预测同一细胞在 ATAC 空间里的大致开放状态。这个指标传统整合工具里没有但对 GLUE 这种“基于链接图做嵌入”的方法来说特别敏感强烈推荐在项目里用上。7. GLUE 与其他多组学整合工具的分工划界现在单细胞多组学整合工具很多新手很容易选错。如果是配对数据很多人会用 Seurat 的 WNN如果是非配对多组学数据有人用 MOFA、Harmony、scVI 或 totalVI。GLUE 在这些工具里的位置其实很特殊它的核心是特征链接图。7.1 工具对比工具输入特点是否利用先验调控知识是否支持调控推断适用场景Seurat WNN配对多组学同一细胞否主要靠细胞间相似性有限Multiome/CITE-seq 常规分析MOFA多组学配对或非配对否线性因子分解弱想找共享/特异因子解释性好Harmony同质性特征有效整合批次否否单组学批次效应校正scVI/totalVI单组学或蛋白RNA否弱表达 / 蛋白表达建模GLUE多组学配对或非配对是强需要调控关系推断的复杂多组学整合这张表不是一个非此即彼的排名而是告诉你怎么按场景选。如果你手里是 10x Multiome且只想看细胞亚群图谱Seurat WNN 确实够用但如果你要把 scRNA、scATAC、DNA 甲基化几个来源的数据整合到同一空间还希望顺带看看“开放染色质变化可能影响哪些基因”那么 GLUE 几乎是目前最直接的工具。MOFA 的解释性更好但它是线性模型对于基因调控这种强非线性过程通常不够Harmony 只处理同构特征根本没法直接用于跨组学特征。7.2 GLUE 的边界什么时候不该用它GLUE 不是万能的。我在以下三种情况会主动建议别人别用它第一你对所用物种的注释比较差比如基因组组装不完整、没有可靠的基因注释和调控注释那链接图本身就建不起来效果还不如用简单的共嵌入方法第二你想解释每个维度对应的生物学过程GLUE 的 latent space 解释性明显不如 MOFA 这种因子模型因为它学出来的是一个整体语义空间不是几个可解释的因子第三你只有单组学数据只想做批次整合——这种时候 Harmony/scVI 更轻、更快没必要为一个不需要调控推断的问题引入整套先验图体系。7.3 我的选择策略做实际项目时我倾向于“先用简单方法再上 GLUE”。先跑一个 Seurat 或 Harmony 的基线整合看看细胞类型分群是否合理如果分群没问题但下游的调控网络一塌糊涂再引入 GLUE 做一轮增强。这样既能快速交付也能在需要时用 GLUE 得到一个更可信的调控推断结果。很多团队把 GLUE 当成“最后一步的裁判”——不是再画一张 UMAP而是用它把 peak 到基因的真实关系理清楚。这个角色定位在我做过的多组学项目里确实没有其他工具能完全替代。如果你正准备开始跑 GLUE我建议不要一上来就堆全部数据。先取一个中等规模的子集把链接表构造好跑一轮训练抽出几个已知的调控关系检查 feature embedding 是否合理。这个过程通常只要半天却能帮你避开后面 90% 的坑。多组学整合本身不难难的是你想清楚“我到底要整合什么、用什么知识去连接它们”。GLUE 提供了一套非常优雅的答案剩下的事就看你怎么用好手里那根锚链了。