企业官网建设流程全解析

做单细胞分析这些年，我一直有个执念：真正能直接跑通的大规模下游代码太少了。大多数论文的代码要么只给定结果，要么绕了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、SoupX；Python那边是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 --id=sample_001 \ --fastqs=/data/fastq/sample_001 \ --sample=sample_001 \ --transcriptome=/ref/refdata-gex-GRCh38-2020-A \ --include-introns \ --localcores=16 \ --localmem=64

有几个地方要特别注意。一是参考基因组，必须下载匹配10x官方的refdata，别自己从UCSC乱抓注释文件，否则基因名对不上；二是--include-introns要不要加，取决于论文代码里用的是什么设置，加了intron会把pre-mRNA也算进去，会导致部分基因表达量系统性偏高；三是肠道样本里经常混入植物、微生物来源的reads，比对后要检查一下人类reads比例，如果明显低，这类样本要标记出来单独讨论。

如果你发现某个样本的比对率特别低，还要警惕是不是barcode泄漏或者建库污染。这个在肠道样本里不算罕见，毕竟取样过程容易混入肠道内容物。处理方式是把样本单独放在一边，先不参与整合，等整体数据跑顺了再回头确认它是不是可以保留。

3. 核心流程与原理拆解

3.1 质量控制：别用固定阈值偷懒

质量控制是单细胞分析里最考验判断力的环节。很多教程就是两个阈值跑到底：min_genes>200、pct_counts_mt<20。放到这个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两条路径，我实际体验如下：

我的经验是先用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_genes=200) sc.pp.filter_genes(adata, min_cells=3) sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], percent_top=None, inplace=True)

pct_mito阈值按样本分组灵活处理。如果你用阈值列表太僵硬，也可以改用"mito比例高于样本中位数三倍MAD"的细胞剔除策略，这在大量样本项目中更好用。接着做标准化、高变基因选择和PCA：

sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, n_top_genes=3000, flavor="seurat_v3") adata.raw = adata.copy() sc.pp.pca(adata, n_comps=50, use_highly_variable=True)

然后跑scVI整合。注意在setup_anndata时指定batch_key为样本ID列：

import scvi scvi.model.SCVI.setup_anndata(adata, batch_key="sample_id") model = scvi.model.SCVI(adata, n_latent=30, n_layers=2) model.train(max_epochs=100, early_stopping=True) adata.obsm["X_scVI"] = model.get_latent_representation()

紧接着用潜在表示建图并聚类：

sc.pp.neighbors(adata, use_rep="X_scVI", n_neighbors=15, n_pcs=30) sc.tl.leiden(adata, resolution=1.0) sc.tl.umap(adata)

如果复现时发现聚类结果偏碎或者偏粗，优先调resolution。leiden聚类的resolution从0.5往上逐步试，直到主要细胞类型都能被分离出来。

4.2 Seurat与Scanpy互转的几个注意事项

有些环节官方代码是在Seurat里完成的，比如SoupX去污染、某些轨迹分析。那我需要把数据在R和Python之间来回搬运。最常用的转换链是SeuratDisk：

library(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对象里可能有多个assay（RNA、SCT），转换时会默认保存第一个assay，你确定自己在用哪个；
• 原来的factor列转成字符串后，整数型metadata会被读成object类型，后续做数值操作前要先转换；
• .uns里如果有特殊对象（比如graph、DBI连接），保存h5ad时可能直接报错，建议先清空或简化obj@misc。

4.3 内存优化与并行调优

160万细胞的矩阵不是拿来就完事儿。我一开始直接跑UMAP，32GB内存瞬间告警，进程被杀，连原始对象都没存住。后来学乖了：

• 只保留需要的列，去掉大段的原始字符型注释；
• 把不需要的矩阵转成稀疏存储，不要轻易转稠密；
• 在跑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 其他小坑：内存爆炸、进程被杀、临时文件把磁盘塞满

最后再补几个环境层面的坑。第一，临时目录别放/tmp，160万细胞的中间文件可能上百GB，最好指定到一个空间充裕的目录。第二，重复运行同一流程时，覆盖文件会引起磁盘碎片，建议每个实验单独开文件夹，统一命名。第三，scvi模型在CPU上训练虽然可行，但别开太多进程，不然内存会被线程池吃空。

写到这里，我顺便把一些高频问题整理成了速查表，方便你在复现卡壳时快速定位问题。

6. 复现之外，还能拿这套代码做什么

6.1 给自己的新数据做注释

这是最直接的收益。论文开源的CellTypist模型训练好了，肠道数据可以直接用现成的模型跑注释，不用自己再从头训练。即使是别的组织（肺、肝、皮肤），只要你能找到一个训练好的参考模型，这套逻辑一样适用。

我自己就试过把新收集的肠道样本直接喂给模型，大部分细胞注释结果和手工判断一致，速度简直可以用秒级形容。需要做的只是确保输入矩阵的基因名和模型特征名对齐。

6.2 与空间转录组数据联合分析

在160万细胞的单细胞参考图谱之上，你可以做空间转录组deconvolution。比如把10x Visium的空间斑点数据映射到参考细胞类型上，估算每个物理位置上细胞类型的组成比例。这样单细胞图谱和空间信息就结合起来了，能回答"哪些细胞靠近哪些细胞"的问题。

6.3 跨数据集、跨疾病的比较

有了统一的参考图谱，你可以把疾病样本（如溃疡性结肠炎、克罗恩病）的细胞向参考图谱映射，计算细胞类型比例偏移、基因表达差异。这套思路在细胞图谱研究中非常流行，实际操作时只需要把参考数据作为anchor，用标准映射算法（比如Seurat的MapQuery或者Harmony）把你自己的数据映射上去。

这篇文章写到这里，其实已经把"复现一套160万细胞单细胞分析"的完整链路拆开了。最后再分享一点个人体会：复现不是目的，理解才是。那些看似简单的QC阈值、批次变量选择、注释验证，才是真正决定分析成败的地方。希望你跑完这套流程后，再面对自己的单细胞数据时，心里能更有底气。