我很少会为一篇Nature论文的代码仓库专门写长文但这次是个例外。Sarah Teichmann团队发表在Nature上的160万肠道单细胞分析GitHub仓库里几乎是“超全代码”式地公开了从QC、批次整合、聚类到细胞注释的全流程管线完整到你根本不用猜步骤。我照着这套流程完整复现了一遍过程中踩了不少坑也摸清了这套基因组学工程的真正设计逻辑。如果你手头有单细胞数据或者想学习大型细胞图谱是怎么搭建的这篇从环境到代码的拆解应该能帮你少走很多弯路。1. 为什么值得死磕这套肠道图谱代码——它的地位和复现门槛1.1 这篇工作解决的是什么问题肠道是一个极度复杂的器官上皮、间充质、免疫、内皮细胞交织在一起不同肠段的功能差异又很大。过去大部分单细胞研究只盯着几万到几十万个细胞样本量小、个体少很难覆盖稀有细胞类型也容易把技术批次当成生物学差异。Teichmann团队这篇工作把规模拉到了160万细胞级别目的是构建一张覆盖多个供体、多个肠段的人类肠道细胞图谱。Sarah Teichmann本人是Wellcome Sanger研究所的细胞图谱计划核心推动者之一她团队做的东西有一个明显特点图谱不是画出来就完事而是尽量把方法论做成可复用、可迁移的流程。所以这篇论文的附加价值比很多同类工作高一大截——论文背后那套代码几乎可以被当成一份“大型单细胞图谱分析模板”来读而不是单纯的复现材料。1.2 代码公开度与可复现性的真实水平单细胞领域的论文代码经常两极分化。有些论文的仓库里只有一个环境配置文件和半截notebook数据动不动就是受控访问复现基本靠猜。但这套肠道图谱不一样分析流程用的是Python生态scanpy为核心兼顾了不同水平的分析者数据发布渠道比较友好原始测序数据和矩阵数据都做了公开步骤被拆成了模块化脚本不是一个几千行的“上帝notebook”每一步都有独立入口细胞注释结合了marker基因和CellTypist模型后者是Teichmann团队自己的工具用起来顺手得反常。我经常说判断一套代码值不值得跟不看它star有多少先看两点一是步骤能不能从原始数据一路跑到最终图二是中间产物有没有被缓存。这套仓库两点都满足。1.3 复现它需要什么基础和时间直说这不是给纯小白准备的教程。你最好已经会用scanpy处理过至少一份10x数据对h5ad对象的操作不陌生知道Leiden聚类和UMAP/tsne到底在干什么。如果这些还没概念建议先拿一个10万细胞以下的公开数据集练一遍再来看它。时间成本方面如果你只用公开矩阵跑核心分析流程普通工作站建议64GB内存以上大约需要两天到一周主要时间耗在环境配置和数据拷贝上。如果你想从原始FASTQ开始重现CellRanger定量那就要准备磁盘阵列和更长的计算时间我个人的建议是除非你想专门学定量流程否则直接从count矩阵开始性价比最高。2. 复现前必须做的三件事仓库地图、环境锁定与数据落盘2.1 先读懂仓库结构再动手很多人在复现时犯的第一个错误是直接双击notebook从头跑到尾。这套流程的仓库目录是按阶段划分的我建议你花一刻钟先把它摸清楚再动手指。我用表格整理一份常见结构不同分支略有差异但思路基本一致目录/文件职责我的理解scripts/核心分析脚本按01、02、03编号严格按编号顺序执行这是管线的主干notebooks/交互式探索脚本用来画图、看中间结果一般不是必须跑通config/路径、样本ID、参数配置文件最容易被忽略但改路径基本都在这data/原始矩阵和中间产物缓存有些体积太大不会进仓库需要单独下载environment.ymlconda依赖清单第一优先要看的文件别自己瞎装包这个结构本身就在传达一个工程化理念可复现的流程不是“一份代码跑完”而是“每个阶段都有独立产物、独立验证节点”。复现的时候不要改动原始脚本的逻辑只改配置里的路径和参数。如果你想在这个流程上做自己的分析不要在原仓库里改直接fork一份出来按自己的数据单独开分支。2.2 环境配置锁版本而不是追最新单细胞分析的环境问题是复现时最大的隐形杀手。scanpy、anndata、scikit-learn、leidenalg这些包版本互通性很差装成最新版反而经常跑不起来。我的经验是以仓库里environment.yml为主如果需要调整记住一个原则锁一套互相兼容的版本组合不要贸然升级某个核心包。我自己跑通的组合大概是这样conda create -n gut_atlas python3.9 conda activate gut_atlas pip install scanpy1.8.2 anndata0.8.0 harmonypy0.0.6 leidenalg0.9.1 pip install celltypist scrublet scikit-misc为什么不追最新版scanpy因为scanpy从1.9开始对部分API有过调整而CellTypist和harmony这些第三方依赖对旧版scanpy的适配更成熟。用python 3.9我个人觉得比3.11更稳不会遇到一些编译型依赖装不上的问题。提示如果你在服务器上装leidenalg遇到编译错误先试试conda install -c conda-forge leidenalg python-igraph别急着用pip硬刚。2.3 数据获取与校验这块卡住了大多数人160万细胞的图谱不是一份文件而是由几十个样本的矩阵组成。Sanger团队一般会把处理后的数据放出来常见形式是h5ad或10x格式的文件夹。下载的时候有两个大坑第一个坑是文件太大。几十GB的h5ad文件用浏览器下载几乎必断。建议用命令行下载工具比如wget加断点续传wget -c https://.../sample_A.h5ad-c参数很重要断了之后能接着下。有条件的话用aspera这类高速传输工具速度能快一个量级。第二个坑是校验。下完文件一定要核对md5或sha256不要觉得没提示就是成功。我遇到过文件显示下载完了但读h5ad时一直报“file is corrupted or truncated”最后发现就是下载阶段出了静默错误。这一步看起来浪费时间实际上省的是你后面排错的几天时间。3. 核心流程逐段走读从原始矩阵到160万细胞的注释图谱3.1 细胞QC与doublet过滤肠道组织的坑比想象中多拿到每个样本的count矩阵之后第一步不是急着合并而是逐样本做QC。这一点是大型图谱流程和小型分析最重要的区别之一小数据集你可以跑完再统一过滤但160万细胞的数据如果直接合并再做QC内存和效率都会出问题而且个别质量差的样本会污染整体聚类结果。QC的核心维度还是那三个每个细胞检测到的基因数n_genes_by_counts、总UMI数total_counts、线粒体基因比例pct_counts_mt。代码骨架长这样import scanpy as sc adata sc.read_10x_h5(sample_A.h5) sc.pp.filter_cells(adata, min_genes200) sc.pp.filter_genes(adata, min_cells3) adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone) # 肠道黏膜组织线粒体比例普遍偏高阈值要比外周血样本宽松 adata adata[adata.obs[pct_counts_mt] 25, :] adata adata[adata.obs[n_genes_by_counts] 500, :]这里有个肠道特有的坑肠道组织本身存在大量菌群和管腔内容物ambient RNA比例比培养细胞高很多线粒体基因比例也普遍偏高。如果你照搬外周血或细胞系的20%阈值会误杀一大批真实的上皮细胞。我在复现时尝试过把阈值放到25%保留下的细胞数量和后续marker表达特异度都更合理。doublet双细胞问题也是矩阵过滤阶段必须处理的。大图谱中双细胞率会随着上机载量升高而上升。建议用scrublet或类似方法逐样本运行不要在全量数据上跑。单个样本的doublet分数分布更干净容易找到阈值。import scrublet as scr scrub scr.Scrublet(adata.X, expected_doublet_rate0.08) doublet_scores, predicted_doublets scrub.scrub_doublets() adata.obs[doublet_score] doublet_scores adata adata[~predicted_doublets, :]提示expected_doublet_rate按10x官方说明上机回收率不同取值也不同1万细胞回收率建议按0.08~0.1估算而不是默认的0.06。3.2 批次整合为什么非做不可Harmony怎么用最稳160万细胞来自不同供体、不同肠段、不同测序批次如果不做整合聚类结果会按样本成团而不是按细胞类型聚集。这一步是整个流程的技术核心也是代码仓库里最值得细读的部分。仓库里用的批次整合方案以Harmony为主。它的核心思路是先把数据做PCA降维然后在主成分空间里迭代识别并矫正“批次相关的子空间”可以在保留真实生物学差异的前提下去除技术差异。相比BBKNN这类基于图的快速整合方法Harmony在保留细胞类型特异表达信号方面更稳相比Scanorama它的内存效率在百万细胞规模下明显占优。标准用法是跑完PCA再接Harmonyimport scanpy as sc adata sc.read_h5ad(gut_atlas_combined.h5ad) # 先对数归一化和高变基因选择再PCA sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, n_top_genes5000, flavorseurat_v3) adata adata[:, adata.var[highly_variable]].copy() sc.tl.pca(adata, n_comps50, svd_solverarpack) sc.external.pp.harmony_integrate(adata, keysample_id)这里有两个参数需要说透。一个是n_comps我用50个主成分作为Harmony输入比默认的20~30个多因为细胞类型多的时候稀有群体的差异往往落在靠后的PC里PC太少容易把稀有细胞类型当噪声抹掉。另一个是harmony_integrate的key必须填你数据里标记样本/批次的obs列名这里是sample_id。整合完成后后续邻居计算和聚类全部要用X_pca_harmony这个嵌入而不是原来的X_pca。很多复现者在这里翻车聚类时忘记改use_rep参数导致整合完全没起作用。后续Leiden聚类时一定记得传入use_repX_pca_harmony。3.3 聚类与分辨率选择图谱不是一次聚类就能成型的聚类阶段放在整合之后图和聚类算法的选择直接影响最终的细胞类群划分。代码仓库里用的是Leiden算法而不是更早期的Louvain。Leiden相比Louvain能避免产生断连的社区生成的聚类结果更稳定在百万细胞级别上跑起来效率也可接受。但我的体会是单次聚类根本达不到图谱级分析的需求。160万细胞里高丰度细胞类型比如T细胞可以轻松聚出一大群而LGR5阳性的肠道干细胞、BEST4阳性上皮细胞这类稀有群体如果分辨率设置不对要么被并入其他群体要么被切成碎片。实际操作时推荐做分辨率扫描而不是拍脑袋定一个数sc.pp.neighbors(adata, use_repX_pca_harmony, n_neighbors30) for res in [0.1, 0.3, 0.5, 0.8, 1.2]: sc.tl.leiden(adata, resolutionres, key_addedfleiden_{res})跑完后对比不同分辨率下的cluster数量同时抽查几个已知marker基因的表达是否在特定cluster里保持一致。判断标准是新提高分辨率能分出有独特marker表达的子群而不是把同一个细胞群随机劈成几份。我自己最后在复现时用的是resolution0.8在主要细胞类型和稀有群体之间取得了比较好的平衡。聚类之后做UMAP可视化时也建议调参但不要盲目追求“好看”的图。n_neighbors30比默认的15更平滑适合大型图谱。min_dist0.5时全局结构清晰但局部细节差我最终选0.3既能看到细胞大群又被分开了又不会把相互关联的过渡态切成碎片。3.4 细胞注释marker为主CellTypist辅助拿到cluster后最核心的问题是这些数字编号的群体到底是什么细胞。仓库的方案是marker基因验证和CellTypist辅助注释双轨并行。先看marker基因。肠道图谱涉及的上皮、免疫、间充质、内皮四大类群每个都有相对特异的标记。我这里列几个我当时反复确认的关键marker细胞类型关键marker基因备注肠上皮干细胞LGR5、OLFM4、SMOC2位置在隐窝底部杯状细胞MUC2、TFF3几乎只存在上皮类群中潘氏细胞DEFA5、LYZ、REG1A小肠特有结肠几乎没有簇细胞tuftPOU2F3、TRPM5、GNAT3极稀有高分辨率下才分得清BEST4上皮细胞BEST4、CA7近年图谱分析新定义的亚型T细胞CD3D、CD3E区分辅助/杀伤看CD4、CD8A浆细胞MZB1、CD79A肠道固有层大量存在用scanpy计算cluster的marker表达时我建议用Wilcoxon检验而不是t检验单细胞表达谱的分布远非正态t检验的p值参考意义不大。sc.tl.rank_genes_groups(adata, groupbyleiden_0.8, methodwilcoxon) sc.pl.dotplot(adata, var_namesmarker_dict, groupbyleiden_0.8)做完marker验证再用CellTypist做监督标签迁移。CellTypist是Teichmann团队自己发布的细胞类型注释工具模型库里包含肠道相关的预训练模型。这一步相当于用已有专家注释数据给新数据“投票”能帮你快速锁定那些marker表达不明显或者没见过的稀有群体。import celltypist model celltypist.models.Model.load(modelHealthy_Adult_Colon) predictions celltypist.annotate(adata, modelmodel, majority_votingTrue) adata.obs[celltypist_labels] predictions.predicted_labels提示CellTypist基于逻辑回归实现运行速度飞快百万细胞几分钟就能出结果。但它的输出只能作为参考最终注释仍然以marker验证为准这两者存在冲突时优先检查marker的绝对表达量是否达到可辨识水平。4. 复现过程实录我遇到的三类卡壳问题与排查链路4.1 数据下载遥遥无期断点续传和校验的细节我复现时最先崩溃的点不在分析而在数据。160万细胞的数据体积大约是几十个样本乘以每个样本几个GB总计几十GB。第一次下载我用了普通的浏览器下载下到一半断了重新再来折腾了大半天还没齐。后来换了wget才解决问题。这里有一个容易被忽略的细节wget的-c参数只能对支持断点续传的服务器生效如果服务器不支持它会从头重新下载。判断服务器是否支持先看HTTP响应头里Accept-Ranges字段如果值是bytes那就放心用-c。下载完成之后的校验我做得有点晚导致一个样本的h5ad文件损坏后我把它顺利读进了分析流程直到合并时候才报矩阵维度对不上。从那以后我的习惯是下载完先跑一遍md5sum和官方给的哈希值比对不一致的直接红牌删除重下不抱侥幸心理。4.2 内存不够用分样本处理再合并的工程化思维我最初的想法很简单所有样本读进来concat成一个大h5ad然后跑完整流程。结果就是在我64GB内存的机器上读入时直接把内存吃满系统开始疯狂swap最后被迫kill进程。这里需要算一笔账160万细胞 × 2万基因的表达矩阵如果按普通稠密矩阵float32存储理论占用是160万 × 20000 × 4字节约128GB远超我的机器内存。实际scanpy用稀疏矩阵存储只记录非零条目但这种规模的稀疏矩阵在中间计算环节仍会频繁扩容内存峰值远高于想象。解决办法其实藏在仓库的设计里每个样本先单独做QC和doublet过滤让每个样本的矩阵降到一个可控规模然后再合并。合并后再做高变基因筛选进一步把基因维度从2万降到5000左右。这一步做完稀疏矩阵的内存压力就小了一个数量级。实测路径是这样的# 每个样本先独立过滤保存为干净的中间文件 import scanpy as sc for sample_id in sample_ids: adata sc.read_10x_h5(f{sample_id}.h5) # ... QC和doublet过滤 ... adata.write(f{sample_id}_filtered.h5ad) # 最后循环读入合并 adata_list [] for sample_id in sample_ids: sample_adata sc.read_h5ad(f{sample_id}_filtered.h5ad) adata_list.append(sample_adata) adata adata_list[0].concatenate(adata_list[1:], joinouter, batch_keysample_id)如果你连合并这一步都撑不过去还有一个杀手锏用anndata的backed模式只读模式打开h5ad不把数据整体载入内存需要时按需读取。但backed模式在scanpy的许多操作里支持不完整我的建议是不到万不得已不要依赖它优先把数据瘦身。4.3 环境冲突和运行崩溃完整的错误排查链路复现过程中我遇到过一个让整个聚类环节崩溃的报错现象是运行sc.tl.leiden时突然抛出一个c-level异常提示类似“segmentation fault”日志没有Python堆栈完全没法定位。我当时的第一直觉是数据有问题回头检查了所有QC参数没发现异常。然后怀疑是leidenalg版本问题就重装了两个版本仍然崩溃。最后我从崩溃时机切入发现它总是在处理某个特定的稀疏矩阵时崩溃。再细查问题出在anndata版本和scipy版本的兼容性上。具体根因是老版本anndata生成的稀疏矩阵对象在某些新版本scipy环境中会触发二进制接口不兼容表现就是偶发性的段错误完全没有Python traceback。修复方案很粗暴但有效建立一个干净的conda环境严格按照python 3.9 scanpy 1.8.2 anndata 0.8.0的组合重新安装崩溃就消失了。这次排查给我的教训是单细胞分析栈的版本兼容性远比普通Python项目敏感出现诡异错误时不要急着怀疑数据和代码逻辑先列一遍依赖版本把所有组合对齐到一组已知能跑的存档上。4.4 全量跑之前先拿5万细胞试跑这是我在复现流程中最后悔没有早点做的事。第一次全量运行时我花了几个小时等待最终因为一个参数错误在聚类阶段失败前面所有计算时间全部作废。第二次我学聪明了先用每个样本随机抽样的方式构建一个5万细胞的小数据集把完整流程从QC一路跑到注释跑通一遍确认所有参数没问题后再启动全量任务。试跑数据集的构建逻辑很简单import scanpy as sc # 抽样时保留样本间的相对比例避免抽样后数据偏向某一个体 combined sc.read_h5ad(combined.h5ad) downsampled [] for sample_id in combined.obs[sample_id].unique(): sub combined[combined.obs[sample_id] sample_id].copy() n int(len(sub) * 0.03) # 每样本抽3% sampled_indices np.random.choice(sub.n_obs, n, replaceFalse) downsampled.append(sub[sampled_indices]) small downsampled[0].concatenate(downsampled[1:], batch_keysample_id)试跑的意义不只是验证参数还能精准估算每一步的耗时和内存峰值。全量任务启动之前你心里应该有一张表哪一步预计多久哪一步是内存瓶颈需要预留多少磁盘空间。这些东西不试跑很难提前知道。5. 这套代码真正教会我的东西不止于复现更是分析工程思维5.1 每个步骤独立成脚本别写“上帝笔记本”看完这套仓库后我做的最大的改变是放弃了自己习惯的单notebook跑全流程方式开始把每一步拆成独立脚本。这看起来只是编码习惯问题实际影响很大。独立的步骤意味着每一步都有稳定的输入和输出调试时不用从头跑。基因表达矩阵归一化是一个小时聚类又是一个小时你不想因为最后画图的参数不对把前面几个小时重新来一遍。仓库的做法是一步一产物每个中间h5ad都保存下来路径清晰命名规范。这个习惯我现在也带进了自己的所有项目里。5.2 版本记录与产物管理复现过程中我注意到仓库对版本记录的重视程度远超常规生信项目。每个产物的文件名几乎都带有数据版本或参数标识比如filtered_v2.h5ad、harmony_res0.8_leiden.h5ad。这种命名方式看似随意却让中间结果之间不会互相覆盖复盘时也能明确知道每个cluster是在什么参数下产生的。我后来在自己的流程里也沿用了这套做法并加了一步在保存中间文件的同时把当前session的版本信息写入一个文本文件。scanpy里一行代码就能输出环境清单import session_info session_info.show()别看这行代码不起眼真到了论文投稿被审稿人追问“某一步用的哪个版本”的时候它就是你唯一的救命稻草。5.3 这套流程对你的数据怎么改造复现不是终点把流程迁移到自己的数据上才是输出。如果你有一个中小规模的肠道或其他上皮组织单细胞数据集完全可以照搬这套流程框架只改几个地方QC线粒体阈值要根据自己组织的特征重新校准Harmony的batch key换成你自己的样本分组列marker基因列表换成对应组织的主流markerCellTypist模型也换成和你的组织类型匹配的模型。唯一需要额外注意的是这套流程默认你的数据有足够的生物学多样性来支撑大规模聚类。如果你的数据只来自一个个体、一个部位强行用160万图谱级别的高分辨率参数很可能会把细胞群切得过度碎片化。参数选择永远服务于生物学问题的尺度这一点是比任何代码都更重要的隐性知识。最后再分享一个小技巧。如果你在服务器上同时跑多个任务读h5ad文件时频繁报错“unable to lock file”并且提示errno 11不要怀疑数据损坏这是hdf5文件锁机制在多进程环境下搞的鬼。解决办法是在启动脚本前加export HDF5_USE_FILE_LOCKINGFALSE让h5py跳过文件锁问题立刻消失。这种细节仓库里不会写但实际跑大型图谱时几乎一定会遇到我现在每次跑数据都先把这个环境变量带上省掉了不少无意义的排查时间。