做单细胞多组学分析这几年最折磨人的环节不是测序而是数据整合。手头既有scRNA-seq又有scATAC-seq或者DNA甲基化数据可这两种数据的特征空间一个在RNA表达量上一个在染色质可及性或甲基化修饰上压根就不是同一个坐标系。硬塞进同一个降维空间里不是被批次效应吞掉就是基因和调控元件对不上号整合出来的结果自己心里都没底。后来用到了北大李程团队开源的GLUEgraph-linked embedding全称《Multi-omics single-cell data integration and regulatory inference with graph-linked embedding》算是把我从“强行凑在一个UMAP里”的泥潭里拉了出来。这篇文章不打算复述论文而是从实际应用的角度把这个方法的核心逻辑、原理细节、实操流程和踩坑记录完整梳理一遍。如果你也在为多组学整合和调控推断头疼这篇应该能帮你省下不少试错时间。1. 为什么多组学整合这么难先搞清楚GLUE在解决什么问题1.1 多组学单细胞数据到底“多”在哪、难在哪要说清楚GLUE的价值得先回到多组学数据本身。单细胞转录组scRNA-seq测的是每个细胞里基因的表达量拿到的是一个细胞×基因的矩阵单细胞染色质可及性scATAC-seq测的是染色质开放区域拿到的是一个细胞×峰peak的矩阵DNA甲基化测的则是每个位点的修饰比例。这三种数据的特征列名完全不同一个是基因名一个是基因组坐标区间一个是CpG位点。问题就出在这你想把所有细胞放到同一个低维空间里比较但每个组学“看到”的特征根本不是同一组变量。拿Seurat那边经典的WNN方法来说它是在每个组学内部先做PCA再用加权最近邻把多个模态的结果揉在一起本质上还是在“同一批细胞”的假设下工作。可实际场景里大部分多组学数据是独立实验产生的细胞群体不完全重合甚至完全没有配对关系这时候WNN那种依赖细胞一一对应的思路就不好使了。更深一层的麻烦在于不同组学的信号类型和信息密度不一样。转录组里能区分出某一群T细胞但到了ATAC的数据里可能是完全不同的一组peak在驱动这群细胞的异质性。你要是硬把它们降维到同一个坐标里很容易出现“转录组里分得开染色质数据里却被拉平”的奇怪现象。GLUE提出的图链接嵌入其实就是针对“特征空间不一致 细胞样本不配对”这两个核心痛点做文章。1.2 已有整合方法为什么不够用从Seurat WNN到LIGER先说Seurat WNN。这个方法在配对多组学数据比如同一个细胞同时测RNA和ATAC上表现还不错因为它天然用了细胞级的对应关系。但一旦面对独立采集的、来自不同批次甚至不同实验室的多组学数据WNN的假设就不成立了你把两组数据硬拼接起来只靠“batch”列去对齐最后要么过度矫正把真实的生物学差异抹掉要么就是那批缺失模态的细胞被硬拉到奇怪的位置。LIGER的思路是用非负矩阵分解把每个组学投影到一个共享的因子空间里但它的对齐强度需要人为调而且对高维稀疏的ATAC矩阵处理效率一般。Harmony适合的就是同一组学不同批次的整合拿到多组学上就得先各自降维再拿降维后的嵌入去做对齐丢失的细节很多。这些方法共同的问题是它们几乎没有利用“基因—调控元件—甲基化位点”之间那些已知的对应关系。而恰恰这些先验知识才是连接不同组学特征的天然桥梁。GLUE敏锐的地方就在于它把这些生物学先验当成一个“图”来用让跨组学的特征之间直接传递信息而不是靠细胞表达量的数值相似度去硬凑。1.3 从ETL视角看数据集成工具泛化与领域差异聊到“data integration”这个词熟悉传统数据工程的人可能第一反应是Pentaho、Kettle那一类的ETL工具。那套体系的思路是把不同来源的数据抽取、清洗、转换最后灌进统一的数据仓库。单细胞多组学整合本质上也是在做类似的事——把基因表达、染色质可及性、甲基化这些异构数据“抽取”出来“清洗”掉批次和技术噪音再“转换”到同一个低维空间里——但生物数据的转换逻辑里多了一层“生物学先验”不能像处理订单表和用户表一样光靠主外键关联。这也解释了为什么传统ETL工具在生物信息领域几乎没有用武之地它们不感知基因组坐标、不识别基因调控关系更不可能帮你推断哪个转录因子调控了哪一群细胞的状态。GLUE的价值在于它把“数据集成”这件事做到了“生物语义”的层面图链接中每一条边都是生物学知识而不是普通的数据表关联字段。理解了这一层后面看它的模型结构就不会觉得抽象了。2. GLUE核心原理图链接嵌入到底在做什么2.1 整体框架双空间变分自编码器加桥接层GLUE的技术底座是变分自编码器VAE这一点和scVI很像但因为面对的是跨组学数据它设计了两个空间原特征空间观测空间和共享的低维潜空间。每个组学数据都配一个专属的编码器和解码器把各自的特征映射到同一个潜空间里去。这一步听起来和“各做各的VAE再向量拼接”差不多但关键在于GLUE在两个组学的特征空间之间加了“图链接”。所谓的图链接graph-linked指的是构建一个跨组学的特征关系图。比如说基因G和它上游某个ATAC峰P存在已知的调控关系那么G和P之间就有一条边。当模型编码RNA数据时基因G的表达信息会通过这条边直接约束ATAC那边峰P的潜空间位置反过来也一样。这样一来两组学虽然在观测特征上完全没有交集但通过这个图它们的特征被从生物学层面“扣”在了一起。把这个思路落到代码层面就是每个组学各自过一套编码器生成均值和方差再从分布里采样得到潜向量。采样之后的潜向量会进入一个“对抗式判别器”用于统一分布同时还有一个“图链接判别器”用来判断潜空间里哪些特征节点应该被连在一起。这两个约束共同作用才让不同组学的细胞嵌入既能对齐又不会丢掉组学特异的信号。2.2 图链接的生物学先验基因-调控元件图谱图链接的具体内容是GLUE从哪儿拿到的论文里用的是基于基因注释和调控元件注释构建的特征关系图。最典型的就是基因转录起始位点TSS上下游的ATAC峰如果一个峰落在某个基因的启动子区域或者已知的增强子区域就把它和那个基因连一条边。对于DNA甲基化数据则是把落在启动子/增强子区域的甲基化探针与对应基因相连。这一步本质上是在把ENCODE、Roadmap这些项目积累的调控元件注释转化为模型输入。这里有一个非常重要的设计细节图链接不是把每个ATAC峰连到它最近的基因就完事而是连到有“调控注释支持”的基因。如果只是按距离连ATAC的peak有一半会落在基因间区连过去根本不是调控关系反而引入了大量噪声。GLUE训练时还会学每条边的权重连错的关系在训练过程中权重会慢慢变低这也是它能做调控推断的基础。实际上图链接的构建质量直接影响整合结果。我一开始偷懒直接用基因上下游2kb范围内的peak连边结果模型训练完ATAC和RNA在UMAP上分布倒是重叠了但细胞类型注释对不上两边的cluster边界明显错位。换成官方推荐的基于peak-genome关联注释构建的图之后整合效果才正常。这个坑后面会细说。2.3 为什么不能“全连”不完全对齐的设计GLUE训练时有个容易让人困惑的地方它刻意不让两个组学在潜空间里“完全重叠”。很多整合方法的目标都是让不同组学、不同批次的细胞在降维图里尽可能混在一起但GLUE反而用了一个“部分对齐”策略。原因是不同组学数据里既包含共享的细胞状态信息也包含组学特异的生物学信息。比如某群细胞的RNA表达差异很大但染色质状态差异相对平缓如果强制对齐负责转录组那部分异质性的信号就会被平均掉最终得到的整合嵌入里这组细胞可能只剩下一团。GLUE的做法是让各组学保留一部分私有的表示只有经过图链接约束的那部分去对齐。这个“不完全对齐”是它优于很多朴素VAE集成方法的地方。在实际操作上控制对齐强度的是模型里的对抗判别器和图链接的权重。判别器太强两个组学会被过度拉近组学特异性信号丢失太弱图链接两边各学各的整合效果聊胜于无。官方默认参数是在大量benchmark上调过的绝大多数情况下不用大改。只有在数据特别复杂、异质性特别强时可以适当把判别器权重调低一些给组学特异信号留更多空间。3. 实操复现用scglue跑通一个标准的多组学整合流程3.1 环境准备与数据格式GLUE官方实现的Python包叫scgluePyPI和conda都能装。我直接在conda里新建了环境Python版本3.9PyTorch 2.0以上基本没问题。保险起见建议按官方GitHub README的install命令装它会自动把依赖的scanpy、bedparse、torch等拉齐。conda create -n scglue python3.9 -y conda activate scglue pip install scglue装好之后数据格式用的是anndata的AnnData对象这一点对用惯了scanpy的人来说非常友好。RNA数据就是常规的细胞×基因矩阵ATAC数据是细胞×peak矩阵peak的名字建议用chrX-起始-终止这种格式后面构建图时省事。另外还需要一个“图”的数据结构scglue推荐用bedparse把基因和peak注释转成图边文件。3.2 图构建与数据预处理图构建这一步官方提供了一个参考基因组注释文件包比如用GENCODE的GTF提取基因坐标用peak文件取交集。实际操作用的是bedparseimport bedparse gene_bed bedparse.read(genes.bed, gtf, tags[gene_id, gene_type]) peak_bed bedparse.read(peaks.bed) # 构建基因-peak图取TSS上下游1kb增强子区域 graph bedparse.build_gene_peak_graph(gene_bed, peak_bed, extend1000) graph.write(gene_peak_graph.tsv)如果你的数据是DNA甲基化就把peak换成甲基化位点文件逻辑一样。预处理方面RNA数据建议用scanpy标准流程过滤低质量细胞、归一化、log1p、筛选高变基因。ATAC数据有个细节值得注意不建议用常规log1p的归一化方式scglue内置了基于TF-IDF的预处理接口对稀疏的ATAC矩阵效果更好。官方pipeline里提供了流程模板我建议大家直接在模板上改比自己从零写要稳得多。另外无论RNA还是ATAC特征名最终必须和图里的节点名完全对应上。我踩过因为基因版本不一致导致图里的边有八成匹配不上的坑训练出来模型效果极差最后只能重新从GTF提取基因坐标统一版本。3.3 模型训练与参数设置数据准备好之后定义一个SCGLUE模型训练起来并不复杂。核心是configure_dataset和fit_SCGLUE这两个接口scglue.models.configure_dataset( rna, NB, use_highly_variableTrue, use_shared_batchesFalse, use_batchFalse ) scglue.models.configure_dataset( atac, NB, use_highly_variableTrue, use_shared_batchesTrue, use_batchTrue ) glue scglue.models.SCGLUE({rna: rna, atac: atac}, graph) glue.fit()几个参数的取舍可以聊一下。use_batch如果设置为True模型会学习每个批次的技术效应但如果数据本身没有明显的批次效应强行设True反而可能把真实生物学差异当批次去矫正。我是把RNA设成False、ATAC设成True因为ATAC数据来自不同构建库技术差异通常更大。模型训练过程中可以隔一段时间看一下训练loss里面几个组成部分的变化趋势重构loss、图链接判别器loss、潜空间判别器loss。如果图链接loss一直震荡不下降八成是图构建时边太少或者连错了特征。训练完成后把RNA和ATAC的细胞嵌入取出来合并拼接直接丢进scanpy的leiden聚类就行。3.4 整合结果的验证与可视化整合后的结果得先从两个维度验证一下一是看两个组学的细胞在UMAP上是否形成一致的细胞类型cluster二是看不同组学的细胞有没有被强行混到一起。前者用常见的marker基因做注释后者要看细胞标签的分布是否合理。还有一招很有用对比整合前后“跨组学的细胞群一致性”。比方说RNA数据里注释出来的CD4T细胞在整合后的ATAC数据里是否也聚在同一个cluster里。有些方法表面上UMAP很漂亮但跨组学邻居一致性很低一到定量评估就露馅。scglue还提供了infer_meta和infer_prior等函数可以算潜空间后验进一步检查嵌合程度。可视化和scanpy无缝衔接emb glue.encode(rna, rna) # 得到RNA细胞的潜空间嵌入 atac_emb glue.encode(atac, atac) combined_emb np.vstack([emb, atac_emb]) import scanpy as sc adata_combined sc.AnnData(combined_emb) sc.pp.neighbors(adata_combined) sc.tl.umap(adata_combined) sc.pl.umap(adata_combined, colorcell_type)如果两个来源的细胞能形成对应的细胞类型团且同一类型的细胞彼此靠近、又不至于糊成一片这个整合基本就成了。4. 调控推断从整合结果里挖cis-regulatory network4.1 解读链接权重GLUE的图链接不是固定的静态关系它在训练时会更新每条边的权重这个权重可以理解为“这条调控关系在当前数据中被支持的程度”。训练完后可以取出学习到的加权邻接矩阵按权重排序就能得到每个基因的主导调控元件。这个过程在论文里对应的是“regulatory inference”的核心输出。我在一个血液数据上跑完GLUE后把权重最高的基因—peak关系按细胞类型拆分发现不少peak实际上落在已知的motif区域里而且和对应细胞类型的谱系特化转录因子高度相关。这比单纯用相关性比如peak和基因的Pearson相关系数推断出的结果要干净得多。因为相关性推断会把大量空间相邻但并无因果调控的peak也算进去GLUE是从调控关系图开始学的天然过滤掉了一部分假阳性。4.2 调控子鉴定有了链接权重可以进一步做调控子regulon分析一个转录因子TF的motif出现在某个peak上而那个peak又通过高权重边连到一个基因那么这条“TF–peak–gene”的关系链就构成了一个候选调控子。实际操作时可以先用HOMER或STREME对高权重peak做motif富集再把富集到的TF跟图链接中的基因对应起来。值得说明的是GLUE本身的输出是“基因—调控元件”的关联强度并不直接给你转录因子结合位点。它和专门的motif分析工具解决的问题不完全一样motif分析回答的是“这个peak上有没有某个TF的可能结合位点”GLUE回答的是“在这个具体数据里这个peak和哪个基因的调控关系被激活了”。把两者结合能得到比单独用其中一个更稳固的结果。4.3 与相关性推断的对比做调控推断的老办法是把同一细胞里的peak可及性和基因表达做相关。这个方法的最大问题在于细胞数有限、数据稀疏相关性的假阳性和假阴性都很高。尤其ATAC数据里零膨胀严重大多数peak在大多数细胞里都是0直接算相关性出来的结果很容易被少数几个高可及性的细胞主导。GLUE用图链接做推断相当于是“先验过滤数据修正”。先验图本身提供了结构性的约束数据再在这个约束里去调整边的权重。所以哪怕两个特征的相关性数值不高只要它们在先验图里有边且权重训练后很高GLUE也会认这个关系反过来相关性很高但图里没有边的“远程关联”则不会被当成直接调控关系大概率是间接效应或者技术噪音。这一点在实际分析里非常有用能显著减少候选调控关系的数量、提升可解释性。5. 常见问题、性能对比与避坑经验5.1 与Seurat、Harmony、scVI的方法学对比为了直观说明GLUE的定位我整理了一个基于实际使用体验的方法对照表方法适用场景是否利用生物学先验跨组学对齐方式调控推断能力Seurat WNN同一细胞配对多组学弱加权最近邻无LIGER多组学、多批次弱非负矩阵分解共享因子无Harmony同组学多批次无模糊聚类迭代校正无scVI/scANVI同组学多批次无深度生成模型无GLUE非配对多组学、调控推断强VAE对抗图链接双约束有从使用体验上讲整合效果层面我在同一套配对数据上对比过Seurat WNN和GLUE两者在RNA主导的细胞类型划分上差异不大但GLUE对ATAC模态相对弱的信号比如某些祖细胞状态保留得更好。而且Seurat WNN拿不到“基因—peak”关系权重也就是说你整合完了还得另找工具做调控推断GLUE则是一套流程全结束。5.2 实操中的坑第1个坑特征名和图节点对不上。这是最隐蔽的坑。RNA用的基因名是SYMBOL图里却用了Ensembl ID或者ATAC的peak坐标是hg19而图是hg38都会导致整合效果大幅下降。建议在构建图之前先用程序检查一下图节点和数据的特征之间的重叠率低于80%就要警惕了。第2个坑高变基因筛选后图边丢失。scRNA预处理时通常筛选高变基因但如果高变基因在调控图里占的比例很低等于把图的一大半连接关系都给切掉了。建议筛选高变基因之后额外保留那些与ATAC peak有连接关系的基因可以显著改善整合效果。第3个坑模型训练显存不足。ATAC的peak数量动辄几万到十几万全量训练时显存容易爆。我通常先对peak做一次覆盖度筛选过滤掉在极少细胞中出现的peak再把高变peak子集用于训练这样显存占用能降低不少速度也快。论文里用的benchmark数据集规模不算大但真实数据动不动几十万细胞这个优化基本是必须的。5.3 什么样的数据适合GLUEGLUE并非万能。如果你的数据只有同一组学、不同批次用Harmony或scVI就够了没必要多此一举引入调控图。如果你有高质量的配对多组学数据同一细胞同时测RNA和ATACSeurat WNN和GLUE都可以前者更轻量后者胜在一套流程还能顺带做调控推断。如果你的数据是非配对的独立多组学或者组学里包含RNAATAC甲基化的三模态数据GLUE几乎是目前最好上手的选择。它能同时整合多模态而且不像某些方法要求各模态特征必须完全对齐。另外GLUE对参考图谱数据也有用——可以在已有注释的转录组图谱基础上把ATAC新数据映射进去做标签迁移这在实际项目里价值很大。从我自己的项目经验看GLUE尤其适合那些“既有公共转录组参考、又有自测染色质数据”的场景。以前想给自测的ATAC数据做细胞注释得先做peak-to-gene关联再用label transfer中间每一步都在丢失信息。现在直接把公共RNA参考图和自测ATAC一起喂给GLUE整合出来的细胞类型标注和marker富集都干净利落。最后的几点实操体会单独说一个训练策略上的细节。GLUE使用了两阶段训练思路先单独对每个组学做预训练让编码器和解码器收敛到合理区间再开启对抗判别器和图链接联合训练。如果一上来就全参数联合训练容易出现生成器还没学好、判别器就已经把梯度带偏的问题最终整合效果会打折扣。官方接口里其实封装了这个过程但如果你自己改模型结构务必保持这个预训练策略。还有关于评价指标我强烈建议别只看UMAP图。UMAP本身会丢失大量局部结构信息有时候看起来分得很开的簇可能只是可视化参数造成的错觉。可以多用跨组学配对分数比如ATAC细胞和RNA细胞之间的kNN重叠率和细胞类型纯度这两个定量指标去做评估。GLUE整合后的数据在这两项指标上通常比Seurat WNN和LIGER更高但前提是你图的构建和数据预处理都做到位了。最后再分享一个小技巧。如果遇到多组学整合后某一群细胞完全被另一组学“吞掉”的情况不要急着调模型参数先去检查这群细胞的marker在另一个组学里是否有对应的特征。比如一群浆细胞的转录组特征非常强但染色质数据里对应的调控元件如果没被纳入图链接那它在ATAC模态里必然缺乏支持的信号整合时自然会被整体压制。这时候合理的做法是在图里补充该细胞群相关的增强子注释而不是去调整损失函数的权重。说到底GLUE最大的价值在于把先验知识引入神经网络的训练过程让“数据驱动”和“知识驱动”不再对立。多组学整合的瓶颈从来不在模型复杂度而在于我们有没有把生物学规律用对。这个思路值得所有做单细胞数据分析的人认真体会。