复现一篇单细胞图谱的分析流程

跟完这篇，您能挑一篇公开了数据的单细胞（scRNA-seq）图谱论文，把它的数据取回来，按它的口径重跑一遍质控、整合、聚类与注释，最后复刻出和原文趋势一致的 UMAP、marker 点图与细胞比例图——并且在对不上的时候，知道该去哪里找原因。

复现一篇单细胞图谱，是把「Methods 里一句话带过的流程」拆成能照做的动作的最好练习。它比复现一张普通差异图更考验人，因为单细胞多了整合、聚类分辨率、注释这几个「有选择、就有分歧」的环节。这篇就带您把一篇图谱论文完整走一遍，用它的几张关键图当靶子。

先读懂论文的方法与数据口径

动手之前，先把论文的方法部分（Methods）和图注读透。单细胞复现能不能对上，八成取决于您是否抄准了下面这几项。建议边读边填一张表：

要素	在论文哪里找	为什么关键
数据编号	Data availability / 附录，形如 GSExxxxx、E-MTAB-xxxxx，或 Zenodo / figshare 链接	决定您去哪儿取数据
样本与批次	Methods、附表	有几个样本、几个批次，直接关系到要不要整合
建库平台	Methods（如 10x 3′ v3、5′、Smart-seq2）	影响基因检出量与质控阈值的合理范围
质控阈值	Methods（nFeature_RNA、percent.mt、双细胞剔除）	决定留下哪些细胞
是否整合、用什么方法	Methods（Harmony、Seurat CCA/anchors、scVI 等）	决定聚类结构，是复现对不上的高发点
主成分数与分辨率	Methods（用了前几个 PC、resolution 取值）	决定分群的粗细
注释依据	Methods、图注、附表的 marker 列表	每个 cluster 该贴什么标签
软件与版本	Methods（如 Seurat v4 / v5、Scanpy 版本）	默认行为会随版本变化

几个读论文时容易漏、却很要命的细节：

1. 质控是「每个样本各卡」还是「合并后统一卡」。 不同样本的测序深度、线粒体本底可能差很多，卡法不同，留下的细胞集就不同。
2. 到底整合了没有、用的哪种方法。 是简单 merge 拼在一起，还是用 Harmony / CCA / scVI 做了批次校正，这一步几乎决定了 UMAP 长什么样，也是复现最容易分道扬镳的地方。
3. 细胞总数与最终的类别数。 论文正文或图注常有一句「共 N 个细胞、M 种细胞类型」。把这两个数记下来，是您复现到位与否最直接的对照。

复现的目标是趋势一致、结论可重现，不是像素级一模一样。单细胞里 UMAP 的具体坐标、cluster 的编号顺序，会随随机种子和软件版本变化，本来就对不齐。只要细胞类型的种类、各类的相对比例、marker 的表达格局对得上，就算复现成功。

取到对应的单细胞数据

拿到数据编号后，去公共数据库把数据取回来。单细胞数据的落脚点通常是 GEO、ArrayExpress、Zenodo / figshare，有些图谱还会上 cellxgene 这类专门门户。先判断走哪条路，能省下大量工作：

• 论文放了处理好的对象（Seurat 的 .rds、AnnData 的 .h5ad，或已经注释好的矩阵）：如果您只想复刻下游的图和结论，直接读这个对象最省事。
• 论文放了每个样本的 filtered 矩阵（10x 的 barcodes/features/matrix 三件套）：想从质控开始完整复现流程，取这个。
• 论文只给了原始测序 reads（fastq）：需要自己跑 Cell Ranger 之类的定量，工作量最大、也最容易引入版本差异。非必要不走这条。

判断和下载的思路：

1. 打开数据编号对应的页面，先看有没有 Supplementary file 里的处理好对象或每样本矩阵。优先取现成的，别从 fastq 重头跑。
2. 数据取回后，务必核对样本数、每个样本属于哪一组是否和论文一致。GEO 的样本描述（characteristics）里通常写清了分组，照它建好样本信息表。
3. 若拿到的是每样本的 10x 目录，逐个读入、各自建对象，之后再合并；这样才能保留「哪个细胞来自哪个样本」，整合时要用到。

# 若论文提供每样本的 10x 矩阵,通常是这样一组文件(具体名字以页面为准)
# sample1/barcodes.tsv.gz  features.tsv.gz  matrix.mtx.gz
# sample2/barcodes.tsv.gz  features.tsv.gz  matrix.mtx.gz
ls sample1/            # 先看一眼:三件套齐不齐?

矩阵行列的整理规范（谁是行、谁是列、基因 ID 用哪套），见 表达矩阵与分组表整理规范；从 GEO 取数据的完整做法见 从 GEO 下载表达矩阵。

读入矩阵、按样本建对象

把每个样本的矩阵读进来，各自建成 Seurat 对象，并在元数据里记下它来自哪个样本、属于哪一组。这一步就是后面整合与分组比较的地基。

library(Seurat)

# 逐个样本读入 10x 目录,建对象时把样本名写进 project
samples <- c("sample1", "sample2", "sample3")
obj_list <- lapply(samples, function(s) {
  counts <- Read10X(data.dir = s)
  seu <- CreateSeuratObject(
    counts       = counts,
    project      = s,
    min.cells    = 3,     # 至少在 3 个细胞里检出的基因才留(官方常用默认)
    min.features = 200    # 至少检出 200 个基因的细胞才留(官方常用默认)
  )
  seu$sample <- s         # 记下样本来源,整合与分组比较都要用
  seu
})

# 合并成一个对象(此时只是拼在一起,还没做批次校正)
seu <- merge(obj_list[[1]], y = obj_list[-1], add.cell.ids = samples)

• min.cells = 3、min.features = 200 是官方教程的常用默认，用来先滤掉几乎不表达的基因和明显是空液滴的条形码。这是惯例、非铁律，稀有细胞多的数据可以放宽。
• 如果论文直接给了合并好的对象，这一步可以跳过，但仍要确认对象里有「样本 / 分组」这类元数据列，否则后面无法按论文口径整合与比较。

质控：按论文的口径卡阈值

单细胞数据里混着空液滴、双细胞和死细胞，必须先清掉。复现阶段的关键是——照抄论文的阈值，别自作主张换一套。常用三个指标：

指标	含义	太高常意味着	太低常意味着
nFeature_RNA	检出基因数	双细胞	空液滴 / 低质量细胞
nCount_RNA	总 UMI 数	双细胞	文库偏浅
percent.mt	线粒体基因占比	细胞破损 / 濒死	——

先算线粒体比例、看分布，再按论文写的阈值卡：

# 人类线粒体基因以 MT- 开头;小鼠用 "^mt-"
seu[["percent.mt"]] <- PercentageFeatureSet(seu, pattern = "^MT-")

# 先看分布,确认和论文描述的量级对得上
VlnPlot(seu, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

# 阈值照抄论文;下面的数字只是示意,请换成论文里写的那套
seu <- subset(
  seu,
  subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 15
)

关键提醒：

• nFeature_RNA 上限、percent.mt < 15 这类都是按数据和组织而定的惯例、非铁律。线粒体阈值常见落在 5%–20% 之间，代谢旺盛的组织（如心、肝）本底就高、应放宽。复现时以论文写的数字为准；论文没写具体值，就先看小提琴图按分布卡，并在复现说明里注明您的选择。
• 论文若提到用 DoubletFinder、scDblFinder 等专门工具剔除双细胞，您也要做同一步，否则细胞数会对不上。质控三指标怎么读、阈值怎么定，展开见 单细胞质控三指标怎么设阈值。

标准化、整合、聚类：还原它的分群

质控之后的这一段，是单细胞复现最容易分道扬镳的地方。核心是忠实还原论文的两处选择：整合方法、以及主成分数与分辨率。

先做标准化和找高变基因，这几步各家流程差别不大：

seu <- NormalizeData(seu)                                   # LogNormalize,官方默认
seu <- FindVariableFeatures(seu, nfeatures = 2000)          # vst,默认取 2000 个高变基因
seu <- ScaleData(seu)
seu <- RunPCA(seu)

接着是是否整合这个岔路口——照论文来：

• 论文没做批次校正（样本同质、直接合并）：跳过整合，PCA 之后直接建近邻图。
• 论文用了整合：按它写的方法做。Harmony、Seurat 自带的 anchor（CCA）、scVI 是最常见的三种，它们都属于被广泛接受的标准做法，选哪种取决于论文，而不是哪种「更好」。

以论文常用的 Harmony 为例，按样本做批次校正，之后的近邻图与 UMAP 都建在校正后的空间上：

library(harmony)

# 按 sample 做批次校正,输出一个名为 "harmony" 的降维结果
seu <- RunHarmony(seu, group.by.vars = "sample")

# 用肘图判断取前多少个维度(拐点之后信息增量很小)
ElbowPlot(seu)

# dims 与 resolution 照抄论文;下面的数字只是示意
seu <- FindNeighbors(seu, reduction = "harmony", dims = 1:30)
seu <- FindClusters(seu, resolution = 0.5)
seu <- RunUMAP(seu, reduction = "harmony", dims = 1:30)

DimPlot(seu, reduction = "umap", group.by = "seurat_clusters", label = TRUE)

对着论文校准三件事：

1. 整合的分组变量：论文按样本、还是按供体 / 批次校正，group.by.vars 就写对应那列。校正的粒度不同，结果会不同。
2. 维度数 dims：用前多少个主成分（或 Harmony 维度），照论文；论文没写就看 ElbowPlot 的拐点，常取 20–30 个，按肘图定、非固定值。
3. 分辨率 resolution：它直接决定分群粗细，调小合并、调大细分。这是复现里最需要对齐的一个参数——先用论文写的值，再对照它报告的 cluster 数微调。

提醒：就算 dims、resolution 都抄对了，您的 UMAP 长相和 cluster 编号顺序也未必和原图一致，这由随机种子和软件版本决定，属正常现象。要看齐的是分出几群、每群对应什么细胞，不是图形摆位。降维聚类与整合的思路总览，见 单细胞分析全流程：Seurat 标准 pipeline。

注释：对上论文的细胞类型

聚类只给出 0、1、2…… 这样的编号，复现的最后一步是把每个 cluster 对应到论文里的细胞类型。做法是找出每群的 marker，再对照论文附表给的 marker 面板来判断——尽量用论文自己列的那套 marker，而不是另起炉灶。

# 找每个 cluster 相对其余细胞高表达的基因(marker)
markers <- FindAllMarkers(
  seu,
  only.pos        = TRUE,   # 只看上调
  min.pct         = 0.25,   # 至少在 25% 的细胞里表达
  logfc.threshold = 0.25
)

# 把论文附表里列的经典 marker 画成点图,肉眼对上每个 cluster
DotPlot(seu, features = c("CD3D", "MS4A1", "LYZ", "NKG7", "PECAM1")) +
  RotatedAxis()

对上之后，把 cluster 编号重命名成细胞类型，力求和论文用的名字一致：

# 顺序对应当前的 cluster 编号 0,1,2,...;标签沿用论文的叫法
new.ids <- c("T cell", "B cell", "Monocyte", "NK cell", "Endothelial")
names(new.ids) <- levels(seu)
seu <- RenameIdents(seu, new.ids)

seu$cell_type <- Idents(seu)
DimPlot(seu, reduction = "umap", label = TRUE)

注意：

• 上面几个基因（CD3D T 细胞、MS4A1 B 细胞、LYZ 单核、NKG7 NK / 杀伤、PECAM1 内皮）是外周血 / 常见组织里被广泛使用的经典 marker；换组织就换论文对应的 marker 面板。
• 注释本质上靠领域知识和 marker 判断，SingleR 这类自动注释可以做参考，但结果需人工核对，不能全盘照收。手动与自动注释各自的做法，见 细胞注释：手动 marker 与自动 SingleR。
• 复现时的对照标准是：您分出的细胞类型，能不能和论文列的那几类对上，且各类的相对比例量级接近。 名字、数量、比例三者大体吻合，注释就算复现到位。

复刻它的关键图

细胞类型注释好，就能照着原文把几张标志性的图搭起来。单细胞图谱最常复刻的是这三张：

图	用什么画	复刻时要对齐
带细胞类型的 UMAP	DimPlot(group.by = "cell_type")	细胞类型的种类、配色
marker 点图 / 热图	DotPlot / DoHeatmap	论文列的 marker、每类的表达格局
细胞比例堆叠图	按 sample×cell_type 计数后用 ggplot2	各组各类细胞的相对占比

细胞比例图是图谱论文的常见主结论之一（「某类细胞在疾病组显著增多」），做法是先按样本和细胞类型交叉计数，再画成堆叠柱：

library(ggplot2)
library(dplyr)

# 每个样本内,各细胞类型的占比
prop <- seu@meta.data %>%
  count(sample, cell_type) %>%
  group_by(sample) %>%
  mutate(freq = n / sum(n))

ggplot(prop, aes(x = sample, y = freq, fill = cell_type)) +
  geom_col() +                       # 堆叠柱,各类占比之和为 1
  labs(y = "细胞比例") +
  theme_bw()

对着原图逐项校准：配色尽量沿用原文（多数图谱给每类细胞固定一个色）、marker 点图里基因的排列顺序照论文、比例图里样本 / 分组的先后也和原图一致。这样两张图并排，趋势对不对得上一目了然。

对不上怎么排查

第一次复现，几乎不会一次对齐——这在单细胞里尤其正常。按下面几类原因从大到小逐个排除，它们对结果的影响通常也是从上往下递减：

可能原因	典型表现	怎么查、怎么办
整合方法 / 分组不同	UMAP 结构、cluster 数明显不同	核对论文是否整合、用哪种方法、按什么变量校正，照它重做
质控口径不同	细胞总数和论文差很多	对齐 nFeature、percent.mt 阈值，以及双细胞剔除是否做了
分辨率不同	细胞类型大体对得上，但分群偏粗或偏细	把 resolution 往论文写的值靠，再微调到 cluster 数接近
样本没对齐	样本数、分组标签对不上	查论文有没有剔除样本、是否只用了部分批次
注释判据不同	类别名对不上、比例有出入	尽量用论文附表的 marker 面板，而不是自选基因
软件版本差异	结构接近但 UMAP 摆位、编号顺序不同	记下 Seurat 等版本；随机种子和版本会让图形摆位变化，属正常

排查时的两条实用心法：

1. 从影响最大的往下查。 整合、质控、分辨率这类「结构性」选择，会让整张图面目全非；随机种子、软件小版本只带来图形摆位的差异。先排大的，再抠小的。
2. 记录您的环境。 跑完 sessionInfo() 把 R 和各个包的版本存下来，既方便自己回溯，也是将来别人复现您结果时的依据。

sessionInfo()   # 记下 R 与 Seurat、harmony 等的版本号

诚实地说：单细胞图谱想和原文的 UMAP 分毫不差，往往做不到——当年的软件版本、整合方法的默认参数可能已难以完整还原。这不是您做错了。复现追求的是「细胞类型、比例趋势、关键 marker 能重现」，而不是复刻每一个点的坐标。把差异查清、能解释清楚，本身就是有价值的产出。若下游还要做组间差异或富集，思路见 统计方法怎么选 与 富集分析：GO / KEGG / GSEA。

在平台上跑一遍

上面这套「取数据 → 质控 → 整合聚类 → 注释 → 复刻关键图」，在 百沐一下 不写代码也能走完：上传论文提供的每样本矩阵（或处理好的对象）和样本信息表，用一句话说清「按论文口径做质控、用 Harmony 按样本整合、resolution 用它写的值聚类，再对照这份 marker 注释、画 UMAP 和细胞比例图」，组学智能体会零代码替您跑完全流程，并连同图形、结果表与可复现的 R 代码、包版本一起交付，方便您和原文逐项比对。完整说明见 组学智能体使用指南。

相关

• 单细胞分析全流程：Seurat 标准 pipeline —— 从矩阵到注释的主流程总览
• 单细胞质控三指标怎么设阈值 —— 三个指标怎么读、阈值怎么定
• 细胞注释：手动 marker 与自动 SingleR —— marker 面板与人工核对
• 从 GEO 下载表达矩阵 —— 把论文的公共数据取回本地