处理单细胞多组学的朋友应该都经历过同样的痛苦手上有同一批样本的转录组和染色质可及性数据甚至在同一个细胞上同时拿到了RNA和ATAC可一打开软件就发现两个数据矩阵的维度对不上几千个细胞、两万个基因对上几十万个峰根本没法直接丢进同一个聚类模型。GLUEgraph-linked embedding这篇论文的做法当时让我眼前一亮它绕开了“强行把特征空间变成一样”的思路而是先把基因、峰、motif之间的已知关系画成一张特征图再让两组数据在图的约束下共同嵌入到同一个低维空间。这篇文章就把GLUE的来龙去脉、原理和上手细节完整聊一遍适合刚接触多组学整合、想弄明白工具背后逻辑的实操型选手。1. 先想清楚多组学整合到底在整什么1.1 组学之间的“语言不通”问题我们做多组学整合本质上想回答一个问题一个细胞在转录组、染色质可及性、甲基化等不同层面上分别是什么样的状态这些状态能不能放在同一个坐标系里比较。道理很简单做起来很麻烦。RNA-seq的feature是大约两万个基因scATAC-seq的feature是几十万甚至上百万个峰甲基化芯片/测序的feature又是几十万个CpG位点。每个feature承载的信息类型也不同表达量是counts可及性是二元的open/closed或counts甲基化率是0到1之间的连续值。更要命的是同一个生物学事件在不同组学里长得完全不一样。比如一个T细胞转录组里看到CD3D高表达ATAC里看到的是CD3D附近的启动子区域开放甲基化数据里则是某些CpG位点的甲基化水平下降。三者本质上都在说同一件事但数值形态、特征名字都对应不上。如果没有一个跨组学的“翻译层”两个矩阵直接拼在一起数学上可行生物学上却没法解释。1.2 为什么不能像传统ETL那样直接拼接做过传统数据集成的人可能会想这不就是个多表关联问题吗给两张表一个共同主键JOIN一下不就好了Pentaho Data Integration这类ETL工具的玩法确实是这样业务数据靠统一主键或者统一ID就能拼接。可组学数据最大的问题在于基因、峰、CpG位点之间根本不存在一一对应的主键关系。一个基因对应若干个启动子近端峰、若干远端增强子峰一个CpG位点可能落在启动子区也可能落在gene body甚至基因间区。如果硬要把ATAC峰注释到基因上再和RNA合并必然会丢掉远端调控信息而远端增强子恰恰是细胞类型特异的染色质调控核心区域。还有一类尝试是把所有组学统一投影到“基因活性分数”Seurat早期的做法就是基于这个思路简单但损失信息。GLUE并不是在细胞层面硬拼而是把“特征之间的关系”本身当成桥梁让特征先对齐细胞再跟着特征走。这个思路上的差异是GLUE和很多传统方案的根本区别。1.3 GLUE给出的答案用图把特征桥接起来GLUE全称是Graph-Linked Unified Embedding作者Cao和Gao在2022年发表在Nature Biotechnology上。名字里的“graph-linked”强调得非常到位它构建了一张跨组学的特征图节点是不同组学的特征边是特征之间已知的生物学关系比如peak落在基因启动子附近、转录因子的motif出现在某peak里、TF和靶基因有调控关系等。这张特征图相当于一个“语义桥梁”。有了桥之后每个组学内部先做各自的降维编码编码过程同时受特征图中的边约束。图里有关联的特征在低维空间里就会被拉近没有关联的特征自然就分开。这样一来不同组学不是被强行塞进同一套坐标而是通过共享的特征空间间接对齐。这个思路最大的好处是即使两个组学的feature名字完全对不上只要它们在图里有边就能在嵌入空间里产生联系之后细胞嵌入也随之对齐。2. GLUE核心原理graph-linked embedding的三步棋2.1 第一步特征图feature graph怎么建GLUE输入的特征图我建议理解成“特征版的知识图谱”。节点是每个组学的feature边是跨组学的生物学关系。常见的关系有几类peak-gene关系某个ATAC峰落在某基因的启动子区域或远端调控区域这是最常用的一类边可以通过基因组坐标直接算也可以从ENCODE等数据库拿。TF-motif-peak关系某个转录因子的结合motif在某个峰序列里出现过说明这个TF可能结合在这个peak上。TF-gene关系转录因子和靶基因之间的调控关系可以从DoRothEA、TRRUST、以及文献积累的GRN里拿。同一组学内部的特征关系比如基因-基因相互作用这种边可选不影响主流程。这张图的构建质量直接决定整合效果。我在实际项目中遇到最典型的失败案例是ATAC峰没有注释到正确的基因组版本基因名又是另一套命名导致peak和gene之间几乎没有边模型虽然能跑但那本质上成了两个独立VAE的拼接完全没有“link”的效果。所以动手前一定要统一参考基因组版本基因注释用同一套GTFpeak坐标也要对应同一版本。2.2 第二步图变分自编码器如何嵌入同一空间特征图只是先验结构真正把数据嵌进去依赖变分自编码器VAE。GLUE的主体框架是图变分自编码器它的思路粗略理解是这样每个细胞的原始特征向量经过编码器网络映射到一个低维隐变量z再通过解码器还原出原始特征空间的分布。关键在解码器这边。普通VAE的解码器是从隐变量z直接映射回原始特征而GLUE在解码之前先用图神经网络GCN/GAT对特征图中的每个特征节点做编码得到每个基因、每个峰的低维表示。之后解码器把细胞隐变量z和特征表示拿去做某种相似度/内积运算重构出细胞在这个特征上的数值。这样一来细胞嵌入和特征嵌入被放在了同一个低维空间里。我个人的理解是特征嵌入描述了“这个基因/峰在这个空间里的位置”细胞嵌入描述了“这个细胞在这个空间里的位置”。如果特征图里基因A和峰B有条边那么它们的特征嵌入会被约束在某个距离内细胞在表达基因A和开放峰B时自然就会被引导到相近的位置。这才是graph-linked embedding的核心逻辑。2.3 第三步链式约束与对抗批次校正如何协同GLUE的目标函数大致可以拆成三块实际训练时这三块不是简单的相加而是加权协同重构损失reconstruction loss保证细胞隐变量能还原原始数据。具体分布会根据组学类型选择RNA-seq通常用负二项分布来建模countsATAC峰可以考虑二项或负二项甲基化用Beta分布。这一点比简单Z-score标准化靠谱得多因为不同组学的数据生成机制不一样。图约束损失graph constraint对特征图中每一条边要求两端特征在嵌入空间中的表示尽量接近。这个约束可以用距离的平方项实现也可以用更复杂的拉普拉斯约束。它决定了“link”的强度也是GLUE防止不同组学被错误对齐的关键。对抗性批次损失adversarial loss跟scVI接近加一个判别器区分细胞来自哪个组学或哪个批次编码器则被训练成让判别器分不出来。这样做的目的是消除组学和批次来源带来的技术差异。这三块需要权衡。图约束太强跨组学确实能合得很彻底但可能把真实的生物学差异也抹掉对抗太强会把不同细胞类型的真实差异也当作“技术变异”给拉平。论文里有一套默认权重但换数据集后我基本都会重新调。2.4 GLUE与Seurat/Harmony/scVI的定位差异很多同学问我有Harmony就能整合为什么还要上GLUEHarmony解决的是“同一组学内多个批次的校正”它没有跨特征空间的建模能力。Seurat WNN虽然能处理多组学但前提是几乎一定是同一细胞测了RNAATAC的paired数据而且本质上还是在特征或细胞维度做加权融合。scVI的VAE框架和GLUE有共同点但scVI通常也是把不同组学的feature全部拼到一个输入空间里组学feature维度悬殊时低维表示很容易被大维度的组学主导。用一张表总结更直接工具核心思路跨组学特征对齐是否需要同一细胞批次校正适用场景Seurat CCA/WNN典型相关分析或加权近邻融合弱通常需要特征同名或先转基因活性最好paired有同一细胞RNAATAC流程成熟Harmony迭代聚类校正嵌入无只处理同一特征空间不需要强单组学多批次整合scVI条件VAE拼接输入特征弱直接拼特征不强制强单组学多批次、大数据量GLUE图变分自编码器特征图桥接强通过feature graph跨特征对齐不强制支持非paired强RNAATAC/甲基化等多组学联合整合、调控推断GLUE不可替代的点在于它允许非配对样本参与跨组学整合。比如你有100个病人的RNA数据和其中50个病人的ATAC数据不需要同一批细胞/同一个人只要有足够的先验特征关系也能放进同一个嵌入空间。这在真实项目里非常实用毕竟不是每个课题都有条件做10x Multiome。3. 实操过程从两个h5ad到整合后的一张大表3.1 环境与安装GLUE的官方实现是Python包scglue依赖PyTorch、scanpy、networkx等。我的建议是单独建一个conda环境不要和日常分析环境混在一起因为scglue对scanpy、numpy的版本比较敏感混装容易把环境搞坏。conda create -n glue python3.9 -y conda activate glue conda install pytorch cudatoolkit11.3 -c pytorch -y pip install scglue scanpyGPU不是硬性要求但我强烈建议用GPU。ATAC峰矩阵随随便便几万到几十万维CPU跑起来太煎熬。装完可以跑一下python -c import scglue; print(scglue.__version__)确认安装成功。我遇到过的情况是PyTorch和CUDA版本不匹配导致模型初始化后原地卡死排查半天最后是把cudatoolkit版本对齐才解决。如果不想折腾CUDA可以退而求其次装CPU版小数据量也能跑。3.2 构建输入Anndata、基因注释、motif注释GLUE的输入是标准的AnnData对象RNA和ATAC各一个都已经做过基本QC。数据怎么预处理有讲究我的经验是RNA数据保留高变基因可以峰值矩阵也建议先做一次标准化和特征筛选但千万不要做传统意义上的批次校正。GLUE自己会处理批次如果你提前用Harmony拉平过一次后面GLUE的对抗学习反而容易过度。关键步骤是给每个AnnData补上特征注释。RNA数据需要有基因注释信息ATAC数据需要峰值注释和motif注释。以人类为例大致流程是这样import scglue import scanpy as sc rna sc.read_h5ad(rna.h5ad) atac sc.read_h5ad(atac.h5ad) # 1. 基因注释 scglue.data.get_gene_annotation( rna, gtfgencode.v38.annotation.gtf.gz, bysymbol ) scglue.data.get_gene_annotation( atac, gtfgencode.v38.annotation.gtf.gz, bysymbol ) # 2. motif注释ATAC需要 scglue.data.get_motif_annotation(atac, specieshuman, motifjaspar)需要说明的是不同版本的scglue接口有差异具体函数签名打开帮助文档看一眼最稳妥。motif注释这一步很关键因为TF-motif-peak这条边是GLUE建立跨组学调控联系的重要来源。如果你用的物种不是人类或小鼠公共数据库覆盖不全特征图会稀疏很多整合效果会打折扣。3.3 模型训练与结果读取特征注释完成后构建feature graph并训练模型。核心流程可以用下面的伪代码表达具体类名和参数以你安装版本为准# 3. 构建跨组学特征图 graph scglue.data.merge( adatas[rna, atac], keys[rna, atac], ongenes, howinner # 取共有基因作为桥接关系 ) # 4. 建立模型并训练 model scglue.models.SCGLUEModel(graphgraph, latent_dim64) model.fit( adatas[rna, atac], graphgraph, max_epochs200, batch_size128, seed0 ) # 5. 导出整合后的嵌入 cells model.encode([rna, atac]) combined scglue.data.merge(cells, keys[rna, atac]) combined.X combined.obsm[X_glue] sc.pp.neighbors(combined, use_repX_glue) sc.tl.umap(combined)这段代码只是主线示意千万不要盲复制。我见过很多人在版本更新后拿着旧教程的类名硬跑结果各种AttributeError。我的习惯是跑之前先help(scglue.models)、help(scglue.data)把接口过一遍花五分钟比报错后查一小时的效率高。训练完成后combined.obsm[X_glue]就是整合后的细胞低维表示可以接scanpy的聚类、UMAP、差异分析流程。也可以把combined导出h5ad作为下游GRN分析、细胞类型注释的输入。3.4 我的训练参数心得GLUE的超参数里我最常调的是三个latent_dim、图约束强度、max_epochs。latent_dim默认值一般是50或100数据量大、组学数量多的时候建议稍微调大一点给模型足够的容量去容纳复杂的特征关系但如果只有RNAATAC两个组学且样本量不大太大会过拟合太小又表达不开关键差异。图约束强度是GLUE的命门。约束太弱跨组学相当于“各跑各的”ATAC和RNA在UMAP上会形成两个“大陆”约束太强则会把不该拉近的细胞类型也搅在一起看起来整合了生物学也毁了。我自己的做法是先用默认参数跑一版看轮廓如果ATAC和RNA明显分群就把图约束权重往上调如果分群很干净但已知marker在两种组学里对不上说明约束过强了往回退。论文里的默认权重是一个起点不是终点。max_epochs建议先跑一个短版比如10个epoch看一眼loss曲线稳定程度再决定。我踩过的最典型的坑是一上来跑200 epoch跑到一半发现学习率太大loss发散白白浪费几小时。先短训调试再全量训练这是我最想强调的经验。4. 整合之后能做什么调控推断与应用落地4.1 跨组学协同聚类与标记基因GLUE整合完成后最直接的价值是跨组学细胞类型对齐。以前RNA和ATAC各自聚类然后靠人工比对marker去猜两组结果的对应关系整合之后两种组学的细胞进的是同一个UMAP同一个细胞类型不管来自RNA还是ATAC都聚在一起。这一步对多组学细胞图谱项目特别有用。实战中我一般会做几个检查第一确认整合后的cluster里RNA细胞和ATAC细胞的比例没有极端偏斜第二抽几个已知的谱系marker基因看它在RNA模态的细胞群中和ATAC模态的peak活性是否一致第三用差异分析找出每个cluster的marker基因和差异可及性peak从两个层面解释这个cluster的生物学身份。GLUE整合做得好的话这两套结果通常是自洽的比如一个T cell clusterRNA里CD3D高ATAC里CD3D附近的调控元件开放度也高。4.2 TF-增强子-靶基因的调控推断GLUE不能直接给你一张GRN但它能把调控推断这件事做得顺畅很多。原因是整合后的数据天然包含了TF-motif-peak-gene这条链式信息特征图里有TF和peak的关系有peak和gene的关系整合后细胞嵌入和特征嵌入又处在同一个空间于是我们可以把“某细胞类型中TF活性高、某peak开放、某靶基因表达上调”这几个事件关联起来。我的流程通常是在整合后的对象上分别跑RNA模态差异表达和ATAC模态差异可及性然后在显著差异peak里做motif富集找到候选TF。接下来回到GLUE的特征图中看这些peak是否和候选TF的靶基因有边连接。如果有就形成了一个很有说服力的调控链路TF在细胞类型A中结合增强子区域该增强子开放度上调靶基因表达随之上升。这比单纯在RNA数据里算TF和靶基因表达相关要可靠因为多了一个染色质层面的证据。4.3 参考映射与新数据投放GLUE训练好的模型可以当成一个参考图谱来用。以后来了一批新的ATAC数据不用重新和RNA一起训练直接用训好的模型把新数据的细胞映射到原有嵌入空间里实现“参考转录组查询染色质”的映射模式。这个功能对临床样本或大规模队列特别有价值。比如你已经有了一套完整的健康人免疫系统RNAATAC图谱后续拿到病人的ATAC数据可以快速映射到图谱中看到病人细胞在哪些细胞类型上偏离了正常状态。需要注意的一点是查询数据的峰注释、基因注释必须和训练时一致否则特征对不上映射结果基本不可信。5. 常见问题与排查技巧实录5.1 训练报错与版本坑scglue和scanpy、numpy、pytorch的版本兼容是个老大难。我自己遇到最多的报错无非几类AttributeError是接口变了KeyError是特征名对不上CUDA out of memory是显存不够。这里分享一个排查思路先在CPU上用最小数据集跑一遍验证代码逻辑无错再切到GPU跑全量。这样能把“代码问题”和“资源问题”分开省很多时间。报错信息常见原因处理办法AttributeError: module scglue has no attribute xxx版本API变动查看官方文档和help按新版写法KeyError: gene_annotation特征注释缺失检查是否运行了get_gene_annotationCUDA out of memory峰值矩阵太大或batch_size过大减小batch_size、latent_dim或先用top可变峰RuntimeError: all elements of input should be between 0 and 1ATAC数据不是二元/计数检查是否做了不合适的标准化5.2 整合效果差过度混淆还是清高过头整合效果差通常分两种。第一种是UMAP上RNA和ATAC各占一坨完全没有融合。这种情况大概率是feature graph太稀疏peak和gene之间的边太少或者motif注释缺失导致TF-peak边根本没建起来。还有一种可能是图约束权重太低。我的排查路径是先看feature graph里的边数够不够再检查图约束损失是否在下降最后一步才是调权重。第二种是整合得“过头”细胞类型完全糊在一起连已知的T细胞、B细胞marker都分辨不出来。这往往是图约束太强把不同细胞类型的真实差异也当作技术差异抹平了。遇到这种情况我会先降图约束再检查批次校正强度。对抗性损失太强也会导致过度混淆因为判别器把不同细胞类型也识别成了需要消除的“批次”。5.3 显存不足与运行时间优化峰值矩阵动辄几十万列显存压力非常大。我的几个实用招数一是事先把低质量峰过滤掉保留全基因组范围内top 5万个可变峰就够大多数分析用了二是把batch_size调小显存不够时从128改到64往往立竿见影三是在模型支持的情况下用混合精度训练四是实在不行就换CPU跑但把max_epochs压到200以内小数据量也能接受。时间方面RNAATAC各一万个细胞50000个peak在单卡V100上大约跑2到3小时这个量级完全可以接受。如果数据量到十万细胞级别建议先在几千个细胞上做参数调试再全量训练避免参数还没调好就白烧一晚上GPU。6. 选型建议什么场景值得上GLUE6.1 场景化对比工具选型不需要赶时髦关键是看你的数据和问题是否匹配。我把常见场景拆一下数据情况核心目标推荐方案只有RNA多批次去除批次效应后集群注释Harmony或scVI同一批细胞的RNAATAC且主要是聚类单细胞多模态注释Seurat WNN流程成熟非配对RNAATAC想把两组装的细胞整合到一起跨组学细胞图谱构建GLUE优先RNAATAC甲基化三组学多组学联合整合调控推断GLUE特征图优势明显已有参考图谱新来一批ATAC/甲基化新数据映射到参考图谱GLUE参考映射或spatial/其他映射工具GLUE尤其适合“两个组学没有天然的配对关系”的项目。比如公共数据库里RNA数据是一批样本ATAC数据是另一批独立样本这种情况下CCA和WNN都很难派上用场而GLUE通过特征图桥接可以正常工作。6.2 我的个人判断如果你的项目只是“同一个multiome数据做个聚类注释”GLUE的收益没那么大用Seurat WNN可能更快更省事。但如果你的目标是构建跨样本、跨组学的整合图谱或者想从ATACRNA联合数据里挖调控机制GLUE的feature graph思路是绕不开的核心资产。我甚至觉得GLUE最大的价值不是那个整合嵌入而是它逼着你去认真整理基因注释、motif注释、peak-gene关系这些先验知识。这堆信息整理清楚之后下游无论跑什么分析都是加分项。最后再分享一点我的操作习惯GLUE跑完之后不管UMAP看起来多干净我都会做一个固定动作抽5到10个已知的谱系marker逐个看它们在RNA模态的表达和ATAC模态的可及性确认两者在对应细胞类型上真的“对齐”了。这一步听起来简单却是很多人偷懒跳过的地方。技术指标再好生物学对不上也是白搭。另外feature graph的构建质量几乎决定了GLUE的上限基因命名不一致、基因组版本混乱、motif数据库选错都会让整合效果大打折扣。我这几年的体会是花半天时间把数据和注释整理干净比调一周参数都管用。这大概是GLUE教给我最重要的一课。