公共数据+复现代码：发育生物学单细胞分析全流程指南

上个月我在整理发育生物学课题时，翻到一篇最新的Cell子刊论文。让我印象最深的并不是它发现了多少新细胞类型，而是打开作者公开的GitHub仓库时，发现从GEO下载原始公共数据，到单细胞质控、聚类注释、拟时序分析，再到最后用来支撑“发育命运决定”故事的那几张关键图表，全部都有可复现代码。这意味着什么？意味着今天做发育生物学研究，拉开差距的不再是“谁手里握着稀缺数据”，而是“谁能把公共数据用更合适的角度重新讲出一个新故事”。

这篇文章想聊的，就是这条完整的思路：如何利用公共数据，以复现代码为主线，搭建一套发育生物学分析流程，并最终把结果组织成一篇逻辑自洽的研究叙事。在正文中，我会结合自己实际跑通这类流程的经验，把数据库选择、分析管线、轨迹推断、调控网络、工程化复现，以及最容易踩到的坑全部拆开来讲。无论你是刚入门的生信小白，还是正在准备一篇基于公共数据论文的在读研究生，这套方法都可以直接套用。

1. 公共数据素材库：发育生物学研究该去哪找数据

1.1 数据库选型：不只GEO一个答案

不少人的第一反应是“公共数据=GEO”，这个思路没错，但发育生物学场景里它往往不是最优起点。发育研究最核心的需求是“动态变化”：同一个器官在不同发育阶段的样子、某个谱系从祖细胞到成熟细胞的过渡状态。这类数据散落在多个平台，而且每个平台的元数据结构差异很大。

我列一张自己常用的素材地图，按发育生物学场景做了归类：

从维护成本来看，我最推荐的是CellxGene + GEO 的组合。CellxGene 的优势在于数据已经经过统一坐标系和基因命名处理，很多发育阶段的细胞图谱都能直接以 h5ad 格式拿到，省去了大量格式清洗工作。但它的缺点也很明显：只能拿处理后的矩阵，拿不到原始测序数据，所以如果你需要重新走比对、定量流程，还是得回到 GEO 下载原始文件。

1.2 判断一个数据集值不值得复现的四个标准

公共数据量很大，但不是每一个都值得花两周时间去复现。我判断标准大概有四条：

• 元数据覆盖度：样本有没有明确的发育阶段、组织部位、生物学重复。很多公共数据集的元数据只写着“胚胎样本”四个字，这种数据跑出来的发育轨迹很难落到具体生物学问题上，故事自然也难讲。
• 时间序列完整性：发育生物学最怕“中间状态缺失”。理想的数据集应该是连续发育阶段，比如 E10.5、E11.5、E12.5 每个阶段至少两个重复。如果只有出生前和出生后两个时间点，轨迹推断大概率会失真。
• 技术平台一致性：同一批数据最好来自同一测序平台、同一建库方案。10x 的 v2 和 v3 化学试剂在捕获效率和基因检出数上差异明显，混在一起分析会产生严重的批次效应。如果作者已经提供了整合好的对象文件，优先用作者给的版本。
• 原始数据可得性：要复现代码，就必须能看到原始数据或至少处理后的 count 矩阵。很多文章只提供 R 语言 Seurat 对象，不提供 raw count，这会让复现者在基因表达定量环节失去控制权。

这四个标准筛下来，真正值得投入的数据集可能只剩两三个。但在这两三个里深入做下去，产出比远高于去“大而全”地把所有数据先跑一遍。我自己的习惯是：筛选阶段花三天，确定主数据集后，后续分析基本全部围绕它可以延伸出两到三个子问题。

2. 从原始公共数据到细胞状态矩阵：管线搭建与选型逻辑

2.1 为什么每个复现项目都要区分上游处理和下游分析

拿到公共数据到正式分析，中间隔着一条很容易被忽略的环节：上游处理。很多拿到公共数据的同学，直接下载作者上传的矩阵或 Seurat 对象就开始聚类。这在很多时候可行，但如果你要讲一个“新故事”，我建议还是把上游掌握在自己手里，因为这决定了后续分析的自主权。

上游处理指的是：原始测序文件（FASTQ）或者原始计数表，经过质量评估、比对、基因定量，得到基因-细胞表达矩阵。这一步使用的工具取决于数据来源。如果是 10x Genomics 数据，可以用 Cell Ranger 或 STARsolo；如果是 Smart-seq2 这类全长转录组，通常用 salmon 或 RSEM 做定量。工具选型的底层逻辑很简单：每类建库方案对应的转录本生物信息特征不同，定量算法必须与之匹配。

下游分析则是从表达矩阵开始：质控过滤、归一化、高变基因识别、降维、聚类、注释、拟时序等。这一层我基本都在 ScanPy 或 Seurat 之间二选一。个人偏好是 ScanPy，其中很重要的一个原因是：ScanPy 的数据结构 AnnData 在生产实践中更容易做模块化处理，中间结果保存、跨步骤共享都更干净。对于复现别人工作场景，ScanPy 还具备一个额外优势——它天然适合从 h5ad 格式直接开始分析，能快速接入 CellxGene 和很多新发布的公共单细胞图谱数据。

2.2 质控指标里的门道：过滤阈值不能照搬

质控是复现项目里最容易被“无脑照搬”的一步。很多人直接照抄教程里的min_genes=200, percent_mt<20，然后发现自己的数据聚类结果乱成一团。原因很简单：不同数据集的技术噪声水平、样本保存条件、组织解离难度差异非常大。发育组织样本尤其容易遇到解离应激，线粒体基因比例普遍偏高。

我会用一套“分布驱动式”的质控策略，而不是死阈值。核心思路是：先看基因数、UMI 数、线粒体占比的分布，再用数据的分布特征来判断阈值。

import scanpy as sc adata = sc.read_10x_h5("sample_filtered_feature_bc_matrix.h5") adata.var_names_make_unique() # 计算质控指标 adata.obs["mt_pct"] = ( adata[:, adata.var["gene_symbols"].str.startswith("MT-")].X.sum(axis=1) ).A1 if hasattr(adata.X, "A") else adata[:, adata.var["gene_symbols"].str.startswith("MT-")].X.sum(axis=1) sc.pp.filter_cells(adata, min_counts=1000) sc.pp.filter_genes(adata, min_cells=10) # 观察分布，再决定阈值，而不是套用固定值 sc.pl.violin(adata, ["n_genes_by_counts", "total_counts", "mt_pct"], multi_panel=True)

这个过程的经验点是：先做宽松过滤（只去掉明显空滴和极低表达基因），然后可视化分布，再根据每个样本的峰值位置收紧阈值。另外提醒一点，发育生物学数据里往往存在生理性高线粒体比例的细胞类型，比如早期胚胎细胞线粒体活跃度就很高。如果机械地按 20% 卡阈值，很可能把关键的前体细胞全部过滤掉——这种情况我踩过一次，重新读完原始论文的补充方法才发现，对方是因为有大量线粒体高占比的“合法细胞”才特意放宽容纳度。

双细胞检测也值得纳入常规流程。Scrublet 是性价比很高的选择，它通过模拟双细胞来估计每个细胞的双细胞分数，阈值可以结合分数分布手动调整。公共数据的双细胞污染情况远比想象中严重，特别是在细胞类型异质性高的胚胎组织中。

3. 把“时间”装进数据：拟时序、RNA速率与轨迹推断的复现细节

3.1 拟时序工具到底在算什么

发育生物学故事的叙事核心永远是“从A到B”。在单细胞数据里，这种“从A到B”是通过拟时序来体现的。拟时序的基本假设是：细胞在分化过程中存在一个连续状态空间，我们观测到的每一个单个细胞，都只是这个连续过程中的一张瞬时快照。

Monocle3、Slingshot 和 scVelo 是三个最常用轨迹推断方案，但它们解决的是不完全相同的问题。

• Slingshot：先在降维空间中找到聚类后的细胞群，再通过最小生成树构造主路径。优点是稳健、抽象程度低，适合你“大概知道有哪些细胞群、想验证它们之间的先后关系”的场景。
• Monocle3：把细胞嵌入到一个由分区图构成的流形中，再学习根节点到终末细胞状态的最短路径。它的最大优势是学习到的轨迹更贴合数据本身的结构，但参数敏感度较高，对根细胞的选择非常敏感。
• scVelo：基于 RNA 剪接动力学（未剪接/剪接 mRNA 比例），直接估计细胞的“速度向量”，不需要事先指定根细胞。它最接近真实的发育时间，但前提是数据里能测到足够的未剪接 reads。

选型上我的经验是：先用 Slingshot 或 scVelo 做“粗确认”，再用 Monocle3 做“精修”。如果你从一开始就用 Monocle3，且对根细胞指定错误，后续所有伪时间轴的生物学解读都会被带偏。

3.2 轨迹推断了，怎样验证它可信

轨迹推断出来不代表故事成立。你还需要三个维度的证据来交叉验证：

• 基因表达连续性：真正沿着轨迹分化的基因，表达变化通常平滑，而不是跳变。选定一组谱系标记基因，把它们按伪时间轴绘制出来，观察是否存在符合预期的梯度表达。
• 与已知生物学一致：如果轨迹推断说 T 祖细胞来源于某个内皮样群体，那么该群体里应该能检测到内皮相关的转录特征。如果完全矛盾，优先怀疑轨迹推断参数。
• RNA 速率方向一致性：scVelo 估算出的速度向量在嵌入空间中应该是“平滑指向终末状态”的。如果速度箭头混乱，一种可能是切换率矩阵学习不充分，另一种可能是数据本身捕获的转录动力学信号太弱。

代码上，用 scVelo 做 RNA 速率的大致流程是这样：

import scvelo as scv adata = scv.read("processed_adata.h5ad") scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000) scv.pp.moments(adata, n_pcs=30, n_neighbors=30) scv.tl.recover_dynamics(adata, n_jobs=8) scv.tl.velocity(adata, mode="dynamical") scv.tl.velocity_graph(adata) scv.pl.velocity_embedding_stream(adata, basis="umap", color="cell_type")

跑通 RNA 速率有个非常容易被忽略的要点：recover_dynamics是一个非常耗时的步骤，在 10 万细胞级别数据上如果不开多线程，可能要跑一晚上。而且它不是每次都收敛良好，如果发现速度流场一片混乱，可以先尝试提高n_top_genes或换用更严格的高变基因选择，而不是急着调速度模型参数。

4. 从聚类到生物学故事：注释、富集与调控网络串联证据链

4.1 细胞类型注释不能只靠一个 marker 列表

发育生物学数据里的细胞类型注释，比成年组织的注释要棘手得多。原因在于发育过程中细胞“身份”是高度可塑的，同一个基因在 E11.5 可能是某个前体细胞的 marker，到了 E14.5 可能变成了另一个分化细胞的 marker。如果只凭一篇文献里的 marker 列表直接打标签，很容易把处于中间状态的细胞强行归到最终命运。

我目前常用的注释策略是“三层证据法”：

1. 已知 marker 集合打分：用 AddModuleScore 或 AUCell 对一组代表不同谱系的 marker 基因打分，得到每个细胞的谱系倾向。
2. 公共参考图谱映射：如果研究组织有已经发表的胚胎图谱，用 celltypist 或 scArches 把查询数据映射到参考数据上，用参考注释作为第二层证据。
3. 差异基因人工复核：自动注释结果最终一定要回到差异表达基因列表里人工过一遍。特别是对发育细胞，我会重点看转录因子表达，比如 Foxc2 出现意味着淋巴内皮前体特征，GATA 家族高表达常与造血谱系相关。

这三层注释经常会出现结果不一致的情况。以我的经验，取交集并保留“状态型注释”（例如“向XX谱系分化的前体”而不是绝对化的“XX细胞”）是更稳妥的做法。

4.2 差异表达、GO/KEGG 与 GSEA：证据链怎么串

注释完成后，故事要往前推一步：这一群细胞和那一群细胞之间，到底发生了什么生物学变化？差异表达是这个环节的基础。但发育故事里，我们看到的大量差异基因其实是连锁反应的结果，它们背后的调控程序才是关键。

所以我的固定流程是：先做差异表达，再从中筛选潜在关键转录因子，然后用 GSEA 看整条通路级别的变化趋势。GO/KEGG 富集解决的是“哪些功能被激活”，GSEA 解决的是“哪些连续基因集整体发生偏移”，后者对发育分化事件更敏感。

import scanpy as sc sc.tl.rank_genes_groups(adata, groupby="cell_type", method="wilcoxon") sc.tl.gsea(adata, groupby="cell_type", gene_sets="KEGG")

实际项目里，别只对“细胞类型”做差异，建议同时对“发育阶段”这个维度的差异做一次。两次差异结果放在一起对比，经常能发现很有意思的规律：某个基因没有在细胞类型间出现差异，但在发育阶段间出现显著变化——这类基因很可能参与的是阶段转换调控，而不是谱系维持，这本身就是很好的故事切入点。

4.3 转录调控推断：SCENIC 是如何让发育故事拥有“因果感”的

发育生物学故事的终极问题是：“是什么决定了这个细胞做出这个命运选择？”单靠差异表达不能回答这个问题，因为大量基因表达变化只是结果，并不是原因。转录调控推断的价值就在这里。

SCENIC 这类工具的做法分两步：第一步，基于共表达和 DNA 基序信息推断出转录因子-靶基因调控模块（regulon）；第二步，用 AUCell 计算每个细胞中每个 regulon 的活性。活性高的 regulon，往往能对应到驱动细胞状态转换的核心转录因子。

# 假设已有 pySCENIC 环境 pyscenic grn loom/adata.loom \ hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather \ --output grn_adj.csv \ --num_workers 8 \ --seed 42 pyscenic ctx grn_adj.csv \ hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather \ --annotations_fname motifs-v9-nr.hgnc-m0.001-o0.0.tbl \ --expression_mtx_fname loom/adata.loom \ --output ctx_output.csv pyscenic aucell loom/adata.loom \ ctx_output.csv \ --output scenic_output.loom \ --num_workers 8

SCENIC 跑完后的关键产出是 regulon 特异性的“活性热点”：比如某个 regulon 在祖细胞群活性很高，在分化成熟细胞中活性消失，这就是一个非常漂亮的“窗口期信号”。把这个信号和拟时序分析结合起来，你就能在轨迹上看到“那扇命运之门在什么时候打开、什么时候关上”。

这里要提醒一下：SCENIC 的计算成本很高，公共数据的细胞量动辄十几万，建议先随机抽 5 万细胞做初步分析，确认 regulon 候选后，再在全量数据上做 AUCell 打分。这样能节省大量时间和内存，且初步分析的结论通常不会因为抽样而改变。

5. 工程化复现：代码结构、环境和“一键跑通”的经验

5.1 一个可复现仓库里，真正值钱的文件是什么

复现代码这件事，很多人理解成“把分析代码传上 GitHub 就行”。但经历过真正复现的人都知道，一个仓库能否被他人顺利跑通，取决于几个隐性因素。

• 依赖锁定：一个environment.yml或requirements.txt只要写了包名而没锁版本，时间过去半年后大概率已经跑不起来。优秀的复现仓库会把核心分析包版本精确锁定。
• 数据路径说明：公共数据下载来源、下载日期、下载时所用版本。GEO 里的文件偶尔会被重新处理，如果不记录，复现结果就会对不上。
• 中间产物管理：分析链路的每一步都保存中间产物，是一个被低估的好习惯。它既方便你自己断点续跑，也方便别人在任意环节接入自己的数据。

我见过的复现仓库最糟糕的状态是：只有一个run_all.R，从头跑到尾，中途任何一个包升级都会导致全盘崩溃。稍微好一点的做法是分模块、分脚本、有明确命名。更好一点，用 Snakemake 或 Nextflow 这类流程化管理工具。

5.2 本地跑通全流程的工程习惯

在实际动手复现的时候，我的目录结构一般是这样：

project/ ├── data/ │ ├── raw/ # 原始公共数据 │ ├── processed/ # 质控后矩阵 │ └── metadata/ # 样本信息 ├── scripts/ │ ├── 01_preprocess.py │ ├── 02_clustering.py │ ├── 03_trajectory.py │ ├── 04_differential_expression.py │ └── 05_regulatory_network.py ├── results/ │ ├── figures/ │ ├── tables/ │ └── models/ ├── env/ │ └── environment.yml └── README.md

每个分析脚本只做一件事，输入输出路径全部通过配置文件或命令行参数传入。一开始做这样的工程拆分，会感觉多花了一点时间，但当数据集从 3 万细胞扩展到 20 万细胞时，你会发现能模块化重跑的价值巨大。

还有一个小习惯：跑完关键步骤后，立刻把输出图的 PDF 和高分辨率 PNG 同时保存。很多期刊要求矢量图，你不可能每次为了插图再重新跑一遍全流程。

5.3 记录复现过程中的“差异”

复现一个公共数据集，最常遇到的问题不是“跑不通”，而是“结果和原文不完全一致”。有几个地方最容易出现数字对不上：

• 基因版本：同一个 gene symbol 在不同参考基因组注释版本下可能对应不同 ID。
• 过滤阈值：文章补充方法里经常只写“按照标准流程”，标准流程本身就有不同版本。
• 批次校正算法：Harmony 的参数不同，整合后 UMAP 明显不同，这不算错误，但会影响图表视觉呈现。

我的建议是把这些差异当作正常现象。复现工作的核心目标不是复刻像素级相同的图，而是确认分析结论的可重复性：原文说这群细胞高表达 A、B、C 基因，并且沿着某条轨迹分化，你的结果能不能重现同样的模式。能，说明这个公共数据的故事是稳固的；不能，恰恰说明你有机会挖出“新故事”。

6. 公共数据复现中的高频坑位：批次效应、内存和随机种子

6.1 批次效应不是靠 Harmony 一个按钮就能解决

公共数据最让人头疼的就是批次效应。发育生物学数据尤其明显：不同胚胎、不同解离批次、不同测序深度，带来的技术噪声经常比真实生物学差异还大。

整合工具的选择要视情况而定：

我的核心经验是：**把“整合”视为一种分析选择，而不是默认动作。**如果每个批次本身包含了足够的完整发育阶段，我甚至不会做显式整合，而是用“合并后聚类 + 查看批次分布”的方式评估。只有当不同批次的细胞类型构成高度倾斜时，才考虑对校正矩阵做处理。

做整合后，一定马上画一张“批次分布 UMAP”，再看一眼“已知 marker 的表达是否仍然保留了原有的空间格局”。批次校正若做得好，细胞类型应形成紧凑的 cluster，而不是按批次颜色分堆。

6.2 内存、临时文件和重复计算

公共数据单细胞图谱的规模，远超多数人的本地电脑能承受的范围。我遇到过 60 万细胞的数据集，UMAP 直接吃了一上午 CPU，临时文件把磁盘写满的事也发生过。几个有效的实践：

• 矩阵格式优先用稀疏存储：ScanPy 的.X如果是 numpy 稠密矩阵，内存占用可能直接翻十倍。用.csr_matrix或 AnnData 默认的稀疏格式能省下大量内存。
• 计算时用backed模式：AnnData 支持backed=True，可以把数据留在磁盘上，按需加载到内存。这对超大数据集非常实用。
• 把跑过的中间结果写成.h5ad保存：一次聚类跑了三小时，随手存盘，后面调参数直接用缓存数据，而不是重新加载原始矩阵再重跑全套。

6.3 随机种子与参数锁定

单细胞聚类分析里，UMAP 和基于图的社区检测（Leiden 聚类）都依赖随机初始化。同一个数据，两次运行得出的 cluster 编号可能不同，这个不是 bug，而是算法特性。

但复现代码时，这就很要命。所以从一开始就要random_seed=42一以贯之。ScanPy 里常见的种子位置包括：sc.pp.neighbors(..., random_state=42)、sc.tl.leiden(..., random_state=42)、sc.tl.umap(..., random_state=42)。我还会在脚本最顶端统一设置全局随机种子：

import numpy as np import random np.random.seed(42) random.seed(42)

另一个容易忽略的是：**Harmony 和 scVI 这类整合方法也有随机成分，运行前一定要记录种子。**否则今天跑出来的 UMAP 图，明天再跑可能就换了布局，审稿人重新运行你的代码后自然会产生质疑。

7. 怎么把复现结果讲成一篇“发育新故事”

7.1 先设计故事骨架，再画图

很多人分析做完了才开始想“我要讲个什么故事”，这其实是本末倒置。以我复现多篇公共数据论文的经验，一个扎实的发育生物学故事骨架通常长这样：

从一张整体细胞图谱切入，告诉读者“这个器官在某个发育窗口里有哪些细胞状态”，然后锁定某一个特别有意思的谱系细胞群，展开轨迹分析，展示从祖细胞状态到成熟状态的连续分化过程。接下来用差异表达和富集分析指出“分化过程中发生了什么”，最后用转录调控推断回答“为什么是这个时候发生。”这条线索对应到图表上就是：图谱总览 → 谱系追溯 → 轨迹推断 → 关键通路变化 → 调控网络窗口。

这基本是“用公共数据讲新故事”最通用的结构，也和我们前面做的分析模块完全对应。动手分析之前，先把这个结构画在纸上，后面每一步都是往里填内容，而不是漫无目的地探索。

7.2 三条做出新意的路线

基于公共数据的发育研究，最怕的事情是“把别人的结果重新命名之后又发一遍”。要让故事有增量，我建议走这三条路线。

第一条：跨物种比较。同一个器官、同一发育阶段，不同物种的数据放在一起比较，可以回答“哪些发育程序是保守的，哪些是物种特异的”。这类分析只需要公共数据，但结论的普适性很强。

第二条：跨时间整合。同一器官不同发育阶段的数据往往分散在多个数据集里，将它们拼接起来，本身就是新数据集。很多时候，单个数据集的样本量不足以支撑完整轨迹，但把三四个数据集整合后，就能看到一段完整的命运变化。

第三条：重注释与重分析。有些公共数据本身已经发表过，但当时的分析限于工具和方法，没有深挖调控机制。多年后你用 SCENIC、RNA 速率这些新工具重新分析，往往能发现新的调控窗口。这类“旧数据+新方法”的组合，是当前低成本撬动高质量故事的捷径。

我举一个具体例子：假设有两个公共数据，分别覆盖了心脏发育早期和晚期阶段，但原始文章都只是做了细胞类型普查。你把他们整合以后，关注一个此前未被重视的、处于“内皮-间充质过渡”状态的细胞群，用拟时序做出这个过渡的完整轨迹，再结合 SCENIC 发现一个调控这个过渡的关键转录因子家族。这样一来，你既没有做任何实验，却讲述了一个关于特定谱系命运决定窗口的完整故事——这就是公共数据真正的价值。

7.3 复现项目如何给别人留下可用资产

最后强调一点：既然你的工作建立在复现和公共数据上，就更有责任把自己的复现资产也完整开放出来。这包括三件事：代码仓库要写清楚运行顺序和依赖环境；每个关键图表要注明对应分析脚本和参数；公共数据下载信息要精确到数据编号、下载日期和版本。

这样做不仅符合学术共享的通行规则，还能让这个项目在公开意义上“立得住”。别人基于你的复现成果继续往下做，你的工作才有可能被引用，进而成为下一段新故事的起点。从我个人的经验来看，项目初期多花一个周末时间把工程细节收拾整洁，后面节省的时间和沟通成本远远超过这点投入。