做单细胞分析这些年我一直有个执念真正能直接跑通的大规模下游代码太少了。大多数论文的代码要么只给定结果要么绕了N道弯跑一次能折腾掉半条命。直到看到Wellcome Sanger研究所Sarah Teichmann团队在Nature发表的肠道细胞图谱研究——超过160万人类肠道单细胞从原始数据处理到细胞注释、到轨迹分析居然把所有代码原原本本公开出来。我当时心想这么全的代码不认真复现一遍真的亏。这篇博文就完整记录我怎么一步步复现这个160万肠道单细胞分析同时把我踩过的坑一并写出来。可能有人会问复现有什么意义意义可太多了。一是这套流程能直接移植到你自己的数据上省掉大量从零搭建的时间二是它暴露了大量在论文里看不到、只有跑代码才能懂的细节——比如某个样本为什么被剔除、某个阈值为什么那样选。这些东西恰恰是做单细胞分析最值钱的经验。1. 项目背景与整体复现思路1.1 这是怎样一项研究肠道不只是一个消化器官它还是人体最大的黏膜免疫界面。肠上皮覆盖面积超过200平方米接触着大量食物抗原和微生物。过去我们只能靠少数组织学标记把肠道分成吸收、分泌、免疫等几大功能块分辨率远远不够。单细胞测序把这件事推进到了逐个细胞点名的水平哪些上皮细胞负责吸收、哪些细胞分泌黏液、哪些免疫细胞驻扎在固有层、它们之间的分子信号是怎么传导的。Sarah Teichmann团队做的这项研究核心是把来自多个队列、覆盖不同发育阶段和多个肠道部位的160万细胞整合到一起建立了一套人类肠道细胞图谱。它回答的问题很具体正常人肠道有多少种细胞状态哪些细胞类型在发育早期就存在、哪些是出生后才上岗的上皮细胞、基质细胞和免疫细胞通过什么分子程序维持组织稳态这些问题的答案直接影响我们对炎症性肠病、肠道肿瘤和感染性疾病的理解。我复现这套分析也是出于同样的考虑先弄清正常肠道里到底有哪些细胞才能讨论疾病状态下哪里出了问题。1.2 复现之前先想清楚三件事第一件事是目标。你是想把论文的图完整复现出来还是想学流程方法用到自己的数据上这两者需要的精力和侧重完全不同。我这次的第一目标是完整跑通代码并理解每个模块的输入输出第二目标才是把关键图复刻出来。第二件事是计算资源。160万细胞的数据量可不是闹着玩的。稀疏表达矩阵大约占用20到40GB内存加上整合、聚类和UMAP建议至少准备64GB内存、8核以上CPU最好有一张支持CUDA的GPU如果跑scVI。没有GPU也不是不行但scVI训练时间会长好几倍跑一个整合模型可能要七八个小时。第三件事是工具栈。这个工作区同时涉及R和Python两套生态。R语言这边是Seurat、harmony、SoupXPython那边是Scanpy、scVI、CellTypist。两手都要硬至少你得能把Seurat对象导出成h5ad或者反过来用Scanpy读h5Seurat。这三件事想明白后再开始动手你会少走很多弯路。2. 环境准备与数据获取2.1 搭建一套不踩坑的conda环境这一步往往是整个复现里最容易出问题的地方。尤其是依赖版本差一个minor版本就可能让某些函数直接报错。我的做法是使用conda/mamba建立两个隔离环境一个给R流程一个给Python流程。R环境里我装的是R 4.2、Seurat 4.3、harmony 1.0、SoupX 1.5。Python环境里是Python 3.9、Scanpy 1.9、scvi-tools 1.0、CellTypist 1.5。版本不需要和官方仓库一模一样但也不能差太远尤其是numpy、scipy、torch这些大件最好按官方lock文件来。我建议用mamba代替conda因为多平台依赖解析速度差异太大。创建环境时直接用mamba env create -f environment.yml mamba activate gut_atlas如果官方仓库没有给environment.yml你就自己列清单但一定要把关键包的版本固定下来。实际中最大的坑往往是Python环境里的torch与CUDA版本不匹配导致scVI训练时直接报libcudart.so not found之类的错误。这种情况建议先用nvidia-smi确认驱动支持的CUDA版本再用pip install torch对应版本去对齐。2.2 数据下载与元数据对齐最费时的部分160万细胞不是一次性采集的而是来自多个公开数据集。常见的数据源包括GEO、ArrayExpress和人类细胞图谱HCA数据门户。文件格式也比较杂有的是10x Genomics的Cell Ranger输出barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz有的是已经处理好的h5ad文件甚至还有loom和RDS。下载数据时紧凑做法是先只下载表达矩阵和配套元数据不要一上来就把几百GB的BAM拉下来。BAM文件之后需要时可以再补但表达矩阵才是分析的主战场。元数据对齐是整个流程中最容易漏掉的一步。论文中的样本编号、供体信息、组织部位通常分散在几个单独的表格里。你需要确认每个样本来自哪个供体、什么年龄段组织部位是回肠、结肠还是直肠样本是否经过筛选、为什么被筛选。我一开始只取了10x的标准输出结果后面做按组织分面绘图时发现一大堆样本的metadata缺失。重新翻README和补充材料一个个对回来白白浪费了两天时间。所以这里强烈建议开工之前先建一张样本映射表里面包含sample_id、文件路径、供体、组织部位、发育阶段、排除标志等字段。后续分析全程围绕这张表来能减少大量重复劳动。2.3 需要从FASTQ重新跑Cell Ranger吗看具体情况。如果官方代码提供的是已经过滤好的表达矩阵你完全可以跳过比对直接从矩阵开始。但如果你下载的是原始FASTQ或者你希望复现比对环节那就要跑Cell Ranger count。以10x v3试剂为例命令大致是cellranger count --idsample_001 \ --fastqs/data/fastq/sample_001 \ --samplesample_001 \ --transcriptome/ref/refdata-gex-GRCh38-2020-A \ --include-introns \ --localcores16 \ --localmem64有几个地方要特别注意。一是参考基因组必须下载匹配10x官方的refdata别自己从UCSC乱抓注释文件否则基因名对不上二是--include-introns要不要加取决于论文代码里用的是什么设置加了intron会把pre-mRNA也算进去会导致部分基因表达量系统性偏高三是肠道样本里经常混入植物、微生物来源的reads比对后要检查一下人类reads比例如果明显低这类样本要标记出来单独讨论。如果你发现某个样本的比对率特别低还要警惕是不是barcode泄漏或者建库污染。这个在肠道样本里不算罕见毕竟取样过程容易混入肠道内容物。处理方式是把样本单独放在一边先不参与整合等整体数据跑顺了再回头确认它是不是可以保留。3. 核心流程与原理拆解3.1 质量控制别用固定阈值偷懒质量控制是单细胞分析里最考验判断力的环节。很多教程就是两个阈值跑到底min_genes200、pct_counts_mt20。放到这个160万细胞项目里这种一刀切的做法很容易把好细胞误杀。不同肠道组织、不同发育阶段线粒体比例的正常范围差异非常大。肠道上皮细胞本身代谢活跃线粒体读数占比经常偏高而某些内分泌细胞因为转录本总量低算出来线粒体比例也很高。你直接设一个20%的硬阈值很可能把整类少见细胞全部过滤掉。我复现时采取的策略是分层QC先按样本为单位画出n_genes、total_counts、pct_mito的分布先过滤掉明显异常的低质量细胞比如基因数小于200、总UMI小于500再对每个样本单独决定线粒体比例阈值参考中位数和MAD中位数绝对偏差对个别偏离很大的样本看看是否存在建库质量问题。双细胞检测也不能少。我是在每个10x样本内部跑Scrublet把predicted_doublet列存进metadata整合后再按需过滤。不建议在全量160万细胞上跑DoubletFinder计算量太大而且模拟双细胞的策略在全量数据上反而容易失真。3.2 整合方法选型Harmony还是scVI样本一多批次效应就躲不掉。不同实验室、不同建库批次、不同测序深度这些技术差异如果不在分析里校正聚类结果就会被批次绑架。在整合环节这个项目里同时存在Harmony和scVI两条路径我实际体验如下特性HarmonyscVI原理迭代软聚类线性嵌入校正变分自编码器非线性生成模型输入PCA嵌入原始稀疏表达矩阵速度快CPU即可慢强烈建议GPU内存较低较高建模潜力可解释性更强能吸收更复杂的技术噪声参数batch_key指定批次变量batch_key指定批次变量可调n_latent/n_layers我的经验是先用Harmony快速做一轮探索性聚类看看主要细胞类型是否浮现如果发现批次簇依然明显再用scVI重新整合。特别要注意的是scVI的输入建议先做Log1p归一化但不建议只保留高变基因后再训练——scVI本身能从全基因中学到更多信息只喂高变基因反而会损失信号。还有一个容易忽略的点整合时要不要保留组织部位和发育阶段这些生物学变量论文代码里的做法是只把样本ID作为批次变量其他生物学因素年龄、组织部位保留在数据中参与聚类。千万不要图省事把所有供体ID也塞进batch_key那会把真实的生物学差异也一并抹掉。3.3 细胞类型注释自动工具与人工验证要配合Sarah团队的工具CellTypist是这次注释环节的核心。它本质上是基于逻辑回归的细胞类型分类器预训练模型已经包含了大量人类免疫、肠道等组织的数据。使用起来很简单celltypist --indata gut.h5ad --model Gut_HCA.pkl --outdir celltypist_output如果模型是针对肠道训练过的那么大部分细胞的注释准确率会非常高。但自动注释不等于完事。我强烈建议保留置信度输出而不是只看最佳匹配标签。当某个细胞在多个类型之间的概率接近时这类细胞往往就是注释灰区需要重点检查。我每次跑完自动注释一定会做以下三件事把我关心的关键marker基因的表达画到UMAP上检查对每个细胞类型做一个平均表达谱和已知的marker列表做对比把低置信度注释的细胞单独拎出来重新聚类分析。比如肠道中的潘氏细胞Paneth cells它们表达LYZ和DEFA5但如果只看LYZ很容易和巨噬细胞混在一起。必须同时看DEFA5、DEFA6这些潘氏细胞特异基因才能区分开。这种组合marker验证是注释环节里必不可少的一步。4. 完整实操流水线与核心代码4.1 用Scanpy完成QC、整合与聚类以下是我在复现中经常会用到的Scanpy核心流水线代码结构可以直接套用到自己的数据上。第一步读取10x输出并进行基础过滤import scanpy as sc import pandas as pd adata sc.read_10x_h5(filtered_feature_bc_matrix.h5) adata.var_names_make_unique() adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.filter_cells(adata, min_genes200) sc.pp.filter_genes(adata, min_cells3) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, inplaceTrue)pct_mito阈值按样本分组灵活处理。如果你用阈值列表太僵硬也可以改用mito比例高于样本中位数三倍MAD的细胞剔除策略这在大量样本项目中更好用。接着做标准化、高变基因选择和PCAsc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, n_top_genes3000, flavorseurat_v3) adata.raw adata.copy() sc.pp.pca(adata, n_comps50, use_highly_variableTrue)然后跑scVI整合。注意在setup_anndata时指定batch_key为样本ID列import scvi scvi.model.SCVI.setup_anndata(adata, batch_keysample_id) model scvi.model.SCVI(adata, n_latent30, n_layers2) model.train(max_epochs100, early_stoppingTrue) adata.obsm[X_scVI] model.get_latent_representation()紧接着用潜在表示建图并聚类sc.pp.neighbors(adata, use_repX_scVI, n_neighbors15, n_pcs30) sc.tl.leiden(adata, resolution1.0) sc.tl.umap(adata)如果复现时发现聚类结果偏碎或者偏粗优先调resolution。leiden聚类的resolution从0.5往上逐步试直到主要细胞类型都能被分离出来。4.2 Seurat与Scanpy互转的几个注意事项有些环节官方代码是在Seurat里完成的比如SoupX去污染、某些轨迹分析。那我需要把数据在R和Python之间来回搬运。最常用的转换链是SeuratDisklibrary(Seurat) library(SeuratDisk) obj - readRDS(gut_integrated.rds) SaveH5Seurat(obj, filename gut_integrated.h5Seurat) Convert(gut_integrated.h5Seurat, dest h5ad)在Python侧用Scanpy读回adata sc.read_h5ad(gut_integrated.h5ad)这里有几个我每次都会犯的错提前提醒你Seurat对象里可能有多个assayRNA、SCT转换时会默认保存第一个assay你确定自己在用哪个原来的factor列转成字符串后整数型metadata会被读成object类型后续做数值操作前要先转换.uns里如果有特殊对象比如graph、DBI连接保存h5ad时可能直接报错建议先清空或简化objmisc。4.3 内存优化与并行调优160万细胞的矩阵不是拿来就完事儿。我一开始直接跑UMAP32GB内存瞬间告警进程被杀连原始对象都没存住。后来学乖了只保留需要的列去掉大段的原始字符型注释把不需要的矩阵转成稀疏存储不要轻易转稠密在跑UMAP前先用PCA降到50维再喂给UMAP的backend如果机器核数多可以适当调高邻居搜索的并行度但注意别把内存吃满。Scanpy有部分函数会隐式地将稀疏矩阵densify比如某些sc.tl.score_genes的实现。如果你发现内存迅速上涨留意警告信息最好手动检查数据类型。细胞数多、基因数多保持adata.X为csr_matrix尤为关键。5. 常见问题与排查技巧实录5.1 为什么UMAP图是一团糊我复现过程中第一次看到UMAP图时简直想砸电脑160万细胞挤成一团没有任何可分的子结构。排查下来问题出在那次scVI训练时把batch_key设成了donor_id而不是sample_id。供体这个变量和生物学差异年龄、性别纠缠不清模型把真实差异当成批次成分抹平了。改回样本ID再训练UMAP结构立刻变得清晰。这里也给一个通用排查顺序先确认批次变量是否选对再检查高变基因是否选对然后试着调低n_neighbors最后再调节分辨率。如果排查完还是一团糊还有一种可能是你的数据里混入了过多低质量细胞回QC环节再清理一遍。5.2 注释结果和论文对不上我复现时的细胞类型比例和论文图左右差了一个数量级排查了很久最后发现是数据筛选不一致。论文把一批低质量样本提前排除了而我为了省事没有下载那份样本筛选清单直接把所有样本都跑进去了。比例差异的直接原因是样本筛选策略不一致跟注释算法没太大关系。核对样本清单的正确姿势是把论文补充材料里的excluded samples表拿过来逐行比对。同样如果把数据源搞混了比如把不同的GEO accession当成同一个数据那结果就更没法对了。所以建议从一开始就建立一份sample_id → 数据文件 → 排除与否的映射表全程用它管理数据。5.3 批次矫正失效的几个信号如果整合后仍然看到按样本聚类的现象有几个常见信号UMAP上出现一个样本独享的独立小簇同一细胞类型表达谱高度依赖技术批次聚类树中上层分支和样本ID对应得过于整齐。遇到这种情况大概率是batch_key设置不对也可能是你用的高变基因包含了太多批次敏感的干扰基因比如线粒体或核糖体基因。只保留高变基因时可以先排除线粒体、核糖体和免疫球蛋白基因减少技术噪声带入整合。另外Harmony对大规模数据的lambda参数也可以适当调低默认值是1有时调成0.5能改善生物学信号保留不过要小心过拟合。5.4 降采样调试一个真香策略复现160万细胞最怕的就是每一步都要等很久。我的做法是先做降采样调试每个初步clusters随机抽一部分细胞比如总共2万细胞先把整个流程跑通、参数调顺、可视化修好再刷全量数据。这一步能把调试周期从小时级压到分钟级极大提升效率。不过降采样也有讲究不能只从某个数据集抽要保证抽出来的子集涵盖所有样本和主要细胞类型否则你看到的调试结果和全量结果会差别很大。我习惯在leiden注释之后先按样本做分层随机抽样这样既保留了批次多样性又保证了细胞类型覆盖。5.5 其他小坑内存爆炸、进程被杀、临时文件把磁盘塞满最后再补几个环境层面的坑。第一临时目录别放/tmp160万细胞的中间文件可能上百GB最好指定到一个空间充裕的目录。第二重复运行同一流程时覆盖文件会引起磁盘碎片建议每个实验单独开文件夹统一命名。第三scvi模型在CPU上训练虽然可行但别开太多进程不然内存会被线程池吃空。写到这里我顺便把一些高频问题整理成了速查表方便你在复现卡壳时快速定位问题。现象可能原因排查手段整合后UMAP无结构batch_key选错改回sample_id重训细胞类型比例与论文不符样本筛选清单不一致逐行核对excluded样本表注释标错类型marker基因重叠组合marker验证不看单一基因内存崩溃矩阵被densify检查adata.X是否为稀疏类型模型训练过慢无GPU或epoch太大用GPU或先降采样调试基因名大量未映射注释版本不一致用模型特征列表做基因对齐6. 复现之外还能拿这套代码做什么6.1 给自己的新数据做注释这是最直接的收益。论文开源的CellTypist模型训练好了肠道数据可以直接用现成的模型跑注释不用自己再从头训练。即使是别的组织肺、肝、皮肤只要你能找到一个训练好的参考模型这套逻辑一样适用。我自己就试过把新收集的肠道样本直接喂给模型大部分细胞注释结果和手工判断一致速度简直可以用秒级形容。需要做的只是确保输入矩阵的基因名和模型特征名对齐。6.2 与空间转录组数据联合分析在160万细胞的单细胞参考图谱之上你可以做空间转录组deconvolution。比如把10x Visium的空间斑点数据映射到参考细胞类型上估算每个物理位置上细胞类型的组成比例。这样单细胞图谱和空间信息就结合起来了能回答哪些细胞靠近哪些细胞的问题。6.3 跨数据集、跨疾病的比较有了统一的参考图谱你可以把疾病样本如溃疡性结肠炎、克罗恩病的细胞向参考图谱映射计算细胞类型比例偏移、基因表达差异。这套思路在细胞图谱研究中非常流行实际操作时只需要把参考数据作为anchor用标准映射算法比如Seurat的MapQuery或者Harmony把你自己的数据映射上去。这篇文章写到这里其实已经把复现一套160万细胞单细胞分析的完整链路拆开了。最后再分享一点个人体会复现不是目的理解才是。那些看似简单的QC阈值、批次变量选择、注释验证才是真正决定分析成败的地方。希望你跑完这套流程后再面对自己的单细胞数据时心里能更有底气。