实际做精准医学和基因组学项目的时候,我经常被问到一个问题:单细胞和空间转录组的数据都拿到了,然后呢?这个问题背后是很多团队的真实困境——测序费用花了不少,fastq堆了一硬盘,但分析管线怎么搭、两组学数据怎么互相印证、最后怎么落到AI模型上,很多人是模糊的。这篇文章是“多组学数据整合与AI建模”章节的开篇,我聚焦在单细胞与空间转录组技术栈上,把从实验原理、数据形态、分析管线、质控逻辑,到多组学整合和AI建模落地的完整路径梳理一遍。我会直接讲我在真实项目里怎么选工具、怎么定参数、在哪里踩过坑,尽量让不同基础的读者都能从中找到自己能用的东西。
1. 单细胞和空间转录组为什么总是绑在一起谈
1.1 单细胞测序解决的核心问题
传统转录组测序(bulk RNA-seq)拿到的是一个组织块里所有细胞的平均信号。问题在于:一个肿瘤组织里可能有癌细胞、T细胞、巨噬细胞、成纤维细胞,它们的基因表达模式完全不同,平均之后很多关键信号会被稀释甚至掩盖。单细胞转录组测序(scRNA-seq)的出现就是为了解决这个问题——把组织解离成单个细胞,给每个细胞单独测序,从而得到“细胞×基因”的表达矩阵。常规10x Genomics平台一次实验可以测几千到上万个细胞,稀疏矩阵里每行是一个细胞,每列是一个基因,矩阵中的数值代表该基因在这个细胞中的转录本捕获量。
但单细胞测序有一个天然代价:组织解离的过程把细胞从原本的“地理位置”上扯下来了。细胞和细胞之间原本挨着还是隔着、肿瘤细胞是浸润在免疫细胞里还是形成一个孤立团块、某个基因的高表达到底发生在组织哪个区域——这些信息在单细胞数据里完全丢失。正因为这个缺点,空间转录组测序被越来越多地纳入实验设计。
1.2 空间转录组:把坐标信息加回来
空间转录组的技术路线与单细胞不同,它不做细胞解离,而是在组织切片上直接进行转录组捕获。目前应用最广的是10x Visium平台:组织切片被贴在带有条形码探针的捕获区域上,每个捕获点位(spot)的直径约55微米,会捕获到该点下方几个到十几个细胞的混合转录组信号,同时记录每个spot在组织上的二维坐标。简单说,空间转录组数据的表达矩阵里,每一行是一个空间spot,每一列是一个基因,另外还有每个spot的x、y坐标,以及对应的组织切片HE染色图像。
除了Visium,还有Slide-seq(分辨率接近单细胞但基因检出量偏低)和Stereo-seq(华大智造,纳米级分辨率)等平台。不同平台的分辨率、通量和成本差异很大,选型时需要根据研究问题权衡。Visium的优势是流程成熟、与10x单细胞数据兼容性好,是目前大多数实验室入门空间转录组的首选。
1.3 整合的必要性:单细胞给分辨率,空间给位置
单细胞数据有细胞类型分辨率但缺位置,空间转录组有位置但每个spot是混合信号。把两者放在一起看,才是完整的组织图景。比如我们研究肿瘤免疫微环境,单细胞数据能告诉我们这个组织里有哪些免疫细胞亚群,而空间转录组能告诉我们CD8+T细胞到底聚集在肿瘤边缘还是浸润进了肿瘤核心;再比如细胞通讯分析,配体-受体互作不仅需要知道哪些细胞表达了配体和受体,更需要知道它们是否在物理上邻近。解决这些问题,本质上是一道“把单细胞分辨率映射回空间位置”的整合题,这也是后面所有分析流程的出发点。
2. 技术栈全景:从fastq到表达矩阵的那套标准管线
2.1 单细胞定量工具:Cell Ranger是标尺,但不是唯一选择
拿到单细胞测序的原始fastq数据后,第一步是把reads比对到参考基因组并完成定量。绝大多数人接触到的第一个工具是10x官方提供的Cell Ranger。它的命令行设计直观,一个count命令就能完成比对、过滤、barcode计数和表达定量:
cellranger count --id=PBMC_sample \ --transcriptome=/data/ref/refdata-gex-GRCh38-2020-A \ --fastqs=/data/fastq \ --sample=pbmc_donor1 \ --localcores=16Cell Ranger最大的优点是“无脑且标准”,输出的filtered_feature_bc_matrix.h5文件可以直接被Seurat或Scanpy读取。它的缺点是资源消耗比较大,并且是一个相对封闭的黑盒逻辑。如果样本量很大或需要更灵活的定制,很多团队会改用STARsolo或者kallisto|bustools。STARsolo用STAR比对器直接生成稀疏表达矩阵,内存占用比Cell Ranger低不少,速度更快;kallisto|bustools走的是伪比对路线,不需要完整比对,速度最快,适合大规模筛选场景。我们在实际项目中,小样本量直接上Cell Ranger,几百个样本以上的大队列会用STARsolo做并行化处理。
2.2 空间转录组的Space Ranger与图像配准
空间转录组的原始数据处理沿用10x的生态,对应工具是Space Ranger。和单细胞流程相比,它多了一个关键步骤:把测序得到的barcode信号与组织切片HE图像的像素位置对应起来。运行count时需要指定slide编号和捕获区域信息,再把组织切片的HE图像一并传入:
space ranger count --id=Visium_sample \ --transcriptome=/data/ref/refdata-gex-GRCh38-2020-A \ --fastqs=/data/fastq \ --image=/data/images/tissue.tif \ --slide=V19B01-073 \ --area=Capture_Area_A1 \ --localcores=16Space Ranger会输出每个spot的坐标、组织下的spot覆盖率,以及和单细胞格式兼容的feature-barcode矩阵。这一步完成之后,空间转录组数据就已经具备了和单细胞数据“对话”的基础语言——同一个参考基因组、同一种稀疏矩阵格式。
2.3 表达矩阵之后:Seurat和Scanpy两个生态的分工
拿到表达矩阵之后,分析的主战场就转移到了R或Python生态。R语言生态里Seurat是绝对主力,它的图形可视化、差异表达分析、富集分析(比如GO富集)体验非常好;Python生态里对应的是Scanpy,处理百万级细胞时内存更友好,而且和后续的机器深度学习工具链衔接更顺。很多单细胞组间GO富集分析的实际操作,目前仍然是R生态最顺手——跑完Seurat的FindMarkers,把差异基因列表喂给clusterProfiler做富集,整套链路成熟可靠。
两个生态之间可以通过sceasy或zellkonverter包做对象互转,AnnData和Seurat对象互相转换,满足混合使用的需求。我的习惯是:小规模探索性分析用R+Seurat,速度直观;数据量大到几十万细胞以上,或者后续要上深度学习模型,就切到Python+Scanpy。
3. 质量控制:真正决定分析成败但最容易被忽略的一环
3.1 单细胞QC的三个标准指标
很多人在拿到表达矩阵后急着跑UMAP,这是比较常见的误区。质量控制如果没做到位,后面的聚类和注释都会失真。单细胞数据质控主要看三个指标:每个细胞捕获到的总UMI数、检测到的基因数,以及线粒体基因的reads占比。线粒体reads占比高的细胞通常意味着细胞膜破损、胞质内容物泄漏,只剩下线粒体还在活跃,属于典型的低质量细胞。
关键的一点是,阈值不应该照搬教材里的固定值。不同组织来源、不同建库条件产生的数据分布差异很大,合理的做法是先画出这三个指标的分布图,观察是否呈现双峰,在低谷处切分。比如外周血单个核细胞(PBMC)数据,基因数一般在500~2500之间,线粒体比例阈值为10%~15%;而某些组织解离后压力较大的样本,线粒体比例阈值可能得放宽到20%以上。
3.2 Doublet问题:被忽视的干扰源
质控里另一个经常被忽略的问题是doublet,也就是两个或多个细胞被同时包裹进同一个液滴,导致一个barcode下面混着多套转录组信号。10x平台的doublet率大约每捕获1000个细胞增加0.8%,样本越大doublet绝对数越多。如果不处理,doublet会形成一批“中间态”聚类,严重干扰细胞类型注释。
检测doublet的成熟工具有DoubletFinder(R包)和Scrublet(Python包),原理都是基于每个barcode的基因表达是否像两张表达谱的混合。建议在聚类之前先跑一遍doublet检测,把高概率doublet的barcode剔除掉。这里有个小技巧:doublet检测的结果通常作为“先验概率”使用,不要一刀切删完,保留概率处于中间地带的细胞,等聚类后再根据marker基因表达情况做人工判断,效果更好。
3.3 空间转录组的QC:维度更多
空间转录组的质控除了UMI数和基因数之外,还要额外关注spot的组织覆盖率、测序饱和度和透化条件的质量。Visium切片上只有组织覆盖区域才有生物学信号,所以第一步是去掉组织边缘外的空spot;测序饱和度低说明当前测序深度还不够把每个spot的转录本都捕获全,继续增加深度有意义,而饱和度已经很高情况下再测就是浪费。组织透化处理时间过长或过短,都会让RNA捕获效率出现偏差,这种实验层面的问题往往会表现为整个切片所有spot的UMI普遍偏低或空间梯度异常。如果遇到这种情况,check一下Space Ranger输出的图像比对结果,往往能看到问题所在。
空间转录组数据经过质控后,还需要做一道和单细胞不一样的流程——对spot的UMI进行归一化和方差稳定化,同时结合组织切片图像,把spot聚类结果映射回空间坐标上,肉眼判断聚类边界是否与解剖结构一致。
4. 空间坐标与邻域结构:空间转录组独有的数据形态
4.1 坐标系统与图结构
处理空间转录组数据,最难适应的是它多了一套坐标系统。每个spot不仅有基因表达向量,还有在组织切片上的row、col坐标。Visium的spot排列在一个六边形栅格上,相邻spot之间的距离是一个固定值。这种结构天然适合用图(graph)来表达:以spot为节点,以物理距离在一定阈值内的spot为边,就构成了一张空间邻接图。邻接图的价值在于,它把“空间邻近”转化成数学上可以计算的结构,后续的空间可变基因识别、细胞生态位分析、甚至图神经网络建模,都建立在这张图上。
Scanpy里可以用squidpy或者GraphTools库来构建邻接矩阵。比如squidpy中用sq.gr.spatial_neighbors(adata)一句就能完成,然后用这个邻接结构去做空间的Moran’s I自相关分析,识别哪些基因的表达存在空间聚集趋势。
4.2 反卷积:从spot到细胞类型
Visium每个spot捕获的是多个细胞的混合信号,直接做细胞类型注释并不现实。通常会采用反卷积(deconvolution)策略:利用同一组织或类似组织的单细胞数据作为参考,推断每个spot中各种细胞类型所占的比例。主流工具包括cell2location、RCTD、SPOTlight等。
cell2location是目前比较稳的选择,基于变分贝叶斯推断,可以把参考单细胞数据的细胞类型特征映射到spot空间,输出每个spot的细胞类型组成比例矩阵。完成后我一般会做一个“可信度检查”:把推断出的免疫细胞丰度与同一个切片HE染色里相应区域的病理结构做肉眼比对。如果肿瘤核心区域在反卷积结果里显示巨噬细胞富集,这就和病理知识互相印证,分析结果才敢往下用。
4.3 空间可变基因与细胞生态位识别
有了邻接结构,就可以回答“哪些基因的表达在空间上不是均匀分布”这类问题。空间可变基因识别工具,如SpatialDE(高斯过程模型)、SPARK(空间自相关模型),可以从几千个基因里筛选出具有显著空间模式的基因。这些基因往往是组织结构形成的分子基础,比如区域特异性标记物、肿瘤边界的分子信号等。
再进一步,细胞生态位指的是由邻近多种细胞类型构成的局部微环境。通过MISTy或Squidpy的细胞邻域分析,可以量化每个spot周围的细胞类型组成,识别出“高T细胞浸润区”“富成纤维细胞区”等不同的生态位类型。有了生态位信息,细胞通讯分析才更有说服力——CellChat或CellPhoneDB这类配体受体分析工具,结合空间距离做过滤,得到的互作关系才是真正有空间意义的。
5. 多组学整合与AI建模:核心不在“拼数据”而在“对齐”
5.1 多组学整合的两种模式
所谓多组学数据整合,现实中通常指两类模式。一类是横向整合:同一批样本做了单细胞转录组、单细胞ATAC-seq(染色质开放性)、蛋白质组等多模态数据,分析目标是把同一细胞亚群的不同模态信号对应起来;另一类是纵向整合:不同批次的单细胞数据来自不同患者或不同实验室,需要把数据合并后进行统一的细胞注释和比较分析。后者在队列研究中极其常见,比如收集100例结直肠癌患者的单细胞数据,首先要做的是消除批次效应,然后再谈疾病亚型差异。
5.2 批次效应:整合绕不开的坎
批次效应来自实验处理、建库试剂批次、测序深度差异,甚至操作人员的不同,是一种技术噪声。如果不加纠正直接合并数据,聚类结果往往先是按批次分开而不是按生物学差异分开。当前主流的批次整合方法有两类:一类是嵌入式方法,代表性的有Harmony和Seurat的CCA整合;另一类是生成模型方法,代表性的有scVI等深度隐变量模型。
Harmony的原理是把数据映射到低维空间后,通过迭代聚类方式最大化生物学差异同时最小化不同批次间的差异,速度快、容易上手,是批量处理时的默认选择。scVI基于变分自编码器,可以学习每个细胞的隐表示,对复杂非线性批次效应的建模能力更强,但训练耗时也更长。这里要特别提醒:批次矫正不能过度。把不同患者的数据矫正到完全重叠,可能会抹掉真实的个体差异和疾病亚型信号。我的做法是矫正后先跑一次聚类,再分批次检查每个聚类的细胞组成比例。如果发现某个聚类里所有细胞都来自同一个样本,那大概率是过矫正了。
5.3 AI建模的几种落地场景
多组学整合的最后一步是AI建模,实际项目中我遇到过比较高频的需求有三种。
第一种是细胞类型自动注释。传统注释依赖人工检查标记基因,费时且主观性强。现在可以用参考数据集训练分类器,或者用Seurat的映射函数把新数据映射到注释好的参考数据上;更进一步则可以直接用单细胞基础模型做零样本注释,比如scGPT和Geneformer这类预训练模型,输入基因表达谱就能输出细胞类型预测。实测下来,预训练模型在常见组织类型上表现不错,但遇到稀有细胞类型或非模式物种时效果打折扣,需要做少量数据微调。
第二种是基因调控网络推断。输入是单细胞表达矩阵和ATAC-seq的染色质开放性数据,输出是转录因子到靶基因的调控关系网络。经典工具如SCENIC利用转录因子结合motif信息推测调控子活性,还有基于互信息的算法如GRNBoost2。这类网络的输出可以帮助锁定疾病驱动转录因子,但网络推测本质上是一种相关性分析,后续需要实验验证。
第三种是药物响应和扰动预测。比如输入某类细胞在药物处理前后的转录组变化,训练模型预测新的药物组合对细胞状态的扰动结果。这个方向现在还在早期阶段,但已经有工作室在用生成模型模拟“如果给这个细胞敲除某个基因,它下一时刻的转录组会变成什么样子”。这类预测一旦稳定,对靶点发现和药物筛选的加速作用会非常明显。
5.4 基础模型时代的单细胞分析
近两年最热的趋势是把大语言模型的思路搬到基因表达数据上。Geneformer、scFoundation、scGPT等模型,用海量单细胞数据做自监督预训练,再用少量标注数据做微调,就能完成细胞类型注释、基因表达预测、扰动效应预测等任务。这个方向的想象空间很大,但目前仍有明显局限:预训练数据的细胞组成偏倚会影响模型在新组织类型上的表现;tokenization策略(把基因怎么映射成离散token)还没有统一标准;模型推理结果的可解释性也远不如传统marker基因验证路径。我的建议是可以用这些模型做初步的注释和预测,但关键结论一定要回到常规分析流程里做交叉验证。
6. 我在真实项目里踩过的坑和现在的推荐路径
6.1 内存和算力的现实约束
单细胞数据进入百万细胞级别后,内存管理是第一个现实问题。Seurat在默认设置下处理20万细胞以上时速度明显下降,有时指令执行后觉得卡死了,其实是在等内存。Scanpy的AnnData基于h5py的稀疏存储结构,对大规模数据友好得多。处理超大样本时我建议使用Scanpy配合Zarr格式做分块读取,或者先用Harmony完成批次整合后用Leiden聚类,这两步都做了并行化优化,可以撑住主流的队列规模。GPU不是刚需但能提速,有一个可见的好处是跑scVI或scGPT这类模型时训练时间从几天缩短到几小时。
6.2 数据格式和版本依赖的坑
R和Python两套生态互转是绕不开的坑。我遇到过一个具体问题:用sceasy把Seurat对象转为AnnData时,由于基因名存在重复(比如线粒体基因MT-与核基因组重复符号),转换后AnnData的var_names不唯一,导致后续Scanpy的many操作报错。后来统一走zellkonverter转换、并先做var_names_make_unique()才解决。另一个常见问题是工具版本不一致:Cell Ranger的参考基因组版本、Seurat版本、Scanpy版本之间偶发兼容性问题,建议在项目开始时就把环境版本冻结,用conda或docker做环境隔离,避免跑完分析后因为版本问题无法复现。
6.3 我推荐的入门路线
如果拿到一套单细胞+空间转录组的匹配数据,我建议按这样推进:先用Space Ranger和Cell Ranger把两家数据分别处理到表达矩阵;用Seurat或Scanpy做单细胞的质控、聚类和注释,得到相对可靠的细胞类型图谱;然后把单细胞数据作为参考,用cell2location对空间转录组做反卷积,得到每个spot的细胞构成;接着构建空间邻接图,分析空间可变基因和生态位;最后把单细胞和空间数据整合到同一个AnnData对象里,结合病理区域标注做AI建模,比如用图神经网络预测某个组织区域属于哪种病理状态。每走一步,都回到原始切片图像上做一次视觉验证,这是我认为整个流程中最重要的工作习惯。
最后分享一点个人体会:单细胞和空间转录组的技术栈已经不再高不可攀,但真正的门槛不是跑通流程,而是每次分析都能守住生物学直觉。多组学整合和AI建模只是放大镜,放大的是数据里的真实信号,同样也会放大数据里的噪声和批次假象。如果你在做类似项目时遇到拿不准的质控阈值或者整合结果和新数据对不上,多半是某个环节的信号失真了,回头从QC开始逐层排查,往往比继续加模型更管用。