Seurat 抽样聚类后，怎样回到全量细胞计算各群占比？

在同一份外周血单细胞数据里，一个细胞群在抽样图上占 11.75%，回到全量数据后占 7.77%。两次计算用的是同一批受试材料，也沿用同一套聚类标签。差别从哪儿来？前一个比例的分母是入选建图的 800 个细胞，后一个是质控后可以参与分析的 2,638 个细胞。做图、给细胞贴标签、统计细胞组成，三步各自用了哪些细胞，决定了最后一列百分数能够回答什么。

下面用 10x Genomics 的公开 PBMC 3k 数据从头做一遍。它足够小，可以同时查看抽样和全量结果。先在代表性细胞上发现群体，再把标签转回原始细胞，最后用标志基因核对。学会这套做法后，同样的检查也适用于更大的图谱。Seurat 的 PBMC 数据教程

用一个小数据集把两份细胞名单分清

下载过滤后的基因与条码矩阵，解压后应得到包含 matrix.mtx、genes.tsv 和 barcodes.tsv 的 hg19 子目录。这份矩阵有 2,700 列，每一列对应一个条码。以下三段代码在同一个 R 会话顺序运行。下文图表记录首次复算，使用 R 4.5.2 和 Seurat 5.4.0。下载区另列 Linux 从头运行此脚本的结果。下载本文所用的 10x 原始计数矩阵

library(Seurat)
set.seed(41)
counts <- Read10X("filtered_gene_bc_matrices/hg19")
obj <- CreateSeuratObject(counts = counts, min.features = 200)
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")
obj <- subset(obj, subset = nFeature_RNA < 2500 & percent.mt < 5)
obj <- NormalizeData(obj, verbose = FALSE)
obj <- FindVariableFeatures(obj, nfeatures = 1500, verbose = FALSE)
c(input = ncol(counts), after_qc = ncol(obj))

本例保留了 2,638 个细胞。这里的基因数和线粒体比例阈值是为了复现这份人外周血示例。换用自己的组织时，先按文库查看检测基因数和线粒体比例的分布，再决定阈值。例如在一个细胞本来就富含线粒体转录本的组织里，直接照搬外周血阈值，可能在抽样之前就删掉了感兴趣的群。使用基因编号或小鼠基因符号时，还要检查线粒体基因的匹配模式是否真的命中了基因。

在抽样细胞上学标签，让模型能投回全量

SketchData() 的 LeverageScore 方法倾向于选入能代表不同表达结构的细胞。一些数量少、表达特征鲜明的细胞因此有机会在图上被看见。这个特点有助于发现小群体，也意味着不能假定每种细胞进入抽样集合的概率相同。数据里没有提供给流程的细胞类型标签，抽样也不会替你判断某一群究竟是什么细胞。作者对 sketch 抽样与投影的说明

obj <- SketchData(obj, ncells = 800, method = "LeverageScore",
                  sketched.assay = "sketch", seed = 41)
DefaultAssay(obj) <- "sketch"
obj <- FindVariableFeatures(obj, nfeatures = 1500, verbose = FALSE)
obj <- ScaleData(obj, verbose = FALSE)
obj <- RunPCA(obj, npcs = 20, verbose = FALSE)
obj <- FindNeighbors(obj, dims = 1:15, verbose = FALSE)
obj <- FindClusters(obj, resolution = 0.6, random.seed = 41, verbose = FALSE)
obj <- RunUMAP(obj, dims = 1:15, return.model = TRUE,
               seed.use = 41, verbose = FALSE)

obj <- ProjectData(obj, assay = "RNA", sketched.assay = "sketch",
                   sketched.reduction = "pca", full.reduction = "pca.full",
                   umap.model = "umap", dims = 1:15,
                   refdata = list(cluster_full = "seurat_clusters"))
DefaultAssay(obj) <- "RNA"

这里有两处值得专门核对。第一，聚类时默认使用 sketch，所以得到的是抽样细胞上的聚类结果。第二，拟合 UMAP 时要保存模型，后面才能把其余细胞放到同一个坐标系里。refdata 指明将抽样群标签转移到全量元数据的新列 cluster_full。一张画得很像原图的 UMAP 只能说明坐标可视化成功，统计前还要检查这列是否覆盖了所有待分析细胞。ProjectData 函数说明

用条码对齐后再做两张频数表

同一个 Seurat 对象可以包含不同细胞集合的 assay。写 ncol(obj) 时很容易忘记当前看的是哪一层，直接从两张表达矩阵取列名更清楚。下面把两套条码分别用于索引同一张元数据表，也要求每个条码都有对应标签。如果检查报错，先找出缺失的条码和标签来源，避免让 table() 默认忽略缺失值后给出看似完整的百分数。

meta <- obj[[]]
full_id <- colnames(obj[["RNA"]])
sketch_id <- colnames(obj[["sketch"]])
full_lab <- as.character(meta[full_id, "cluster_full"])
fit_lab <- as.character(meta[sketch_id, "seurat_clusters"])
stopifnot(!anyNA(full_lab), !anyNA(fit_lab))

lev <- sort(unique(c(full_lab, fit_lab)))
n_fit <- table(factor(fit_lab, levels = lev))
n_full <- table(factor(full_lab, levels = lev))
data.frame(cluster = lev, sketch_n = as.integer(n_fit),
           full_n = as.integer(n_full),
           sketch_pct = round(100 * as.numeric(n_fit) / length(fit_lab), 2),
           full_pct = round(100 * as.numeric(n_full) / length(full_lab), 2))

看群 3，如果把抽样比例直接写进结果，会多出约四个百分点。看群 0，变化方向却相反。表的左侧取建图标签，右侧取投影标签，因此两列差异还可能包含边界细胞被重新赋值的影响。想拆开这两个因素，可以只对入选的 800 个条码比较 table(meta[sketch_id, "seurat_clusters"], meta[sketch_id, "cluster_full"])。集中在对角线外的标签对值得优先回看。同一次分析里的数字编号只是索引，换随机种子或聚类参数后，要按标志基因重新识别对应群。

什么时候可以把“群 2”写成一种细胞？

本例在全量 RNA 的归一化表达中按投影标签取平均，MS4A1 和 CD79A 均在群 2 最高。两个标志物给出了把该群作为 B 细胞候选的依据。可以再运行 DotPlot(obj, features = c("MS4A1", "CD79A", "CD3D", "NKG7"), group.by = "cluster_full")，同时查看表达细胞比例与群内平均表达。如果只有极少数细胞高表达一个标志物，或者两个通常互斥的谱系标志在同一批细胞里同时很强，就要回看单细胞表达、双细胞和群边界，再给群体命名。

用于论文比较时，还要把分母推进一步。本例的 7.77% 表示该文库通过质控、纳入投影的细胞中，7.77% 被分到了群 3。组织解离、捕获和质控都会影响这个比例。对于有多位供体的研究，应先在每位供体内部统计各群数量和占比，同时保留每位供体的总细胞数。这样能看见变化是否在多位供体中同向，也能发现某个大文库对合并图的影响。条件间的推断需要以独立生物学样本为单位。

一个文库捕获的细胞多，会怎样改变合并比例？

把这个例子换成两个供体，问题会更容易看清。假设甲提供了一千个通过质控的细胞，其中一百个属于目标群。乙只提供一百个，其中四十个属于目标群。把细胞先合在一起，会算出一千一百个细胞中有一百四十个属于目标群，约占百分之十二点七。先在每个供体内计算再取平均，得到的却是百分之二十五。这两个数字用了相同的标签，区别完全来自给供体分配的权重。前者描述这批捕获细胞的组成，后者描述这两位供体各占同等权重时的平均组成。给图标注“平均占比”之前，应能说出取的是哪一种平均。

ex <- data.frame(donor = c("甲", "乙"), target = c(100, 40), total = c(1000, 100))
transform(ex, pct = round(100 * target / total, 1))
round(100 * c(pooled = sum(ex$target) / sum(ex$total),
              donor_mean = mean(ex$target / ex$total)), 1)

在自己的元数据里，应当保留“供体、组织来源、采样时间、文库”四者中实际存在的层级。假如同一供体的一份样本分在两个上机孔，先在样本内合并技术重复的细胞计数。把两个孔当成两个受试者会凭空增加样本量。相反，同一供体治疗前后的两次取样应保留时间点，可以画一条线连接同一人的两个点。这样能看出是大多数人都向同一方向移动，还是几条线很陡、其他人几乎没有变化。对于只采了一份文库的外周血示例，只报告组成描述即可。

还要区分零与缺失。如果某位供体通过质控的细胞中没有一个被标为目标群，它应出现在占比表里，目标数填零。如果该供体根本没有可用的细胞，总细胞数为零，占比没有定义。把后一种情况也写成零，会使捕获失败的文库看起来像是完全没有这类细胞。整理频数表时先保留所有已知样本，再标出哪一行能计算比例。决定是否剔除有效细胞极少的样本之前，先看每一组将损失多少位独立供体。

给少数群体做一次针对性的投影复核

可以把抽样条码与其余条码分成两组，在同一个候选细胞群内比较标志基因检出率。分子是至少检出一次该标志的细胞数，分母是该部分被分进候选群的细胞数。若抽样中两个相互补充的标志都很清楚，投影后新增的一批细胞却几乎都只有管家基因，先核对这些细胞的文库、检测基因数和测序深度。也可以选取相邻群的反向标志，一起画在原始表达层上。只在用于建图的低维坐标中看点与点的远近，无法替代这一步表达复核。

另一个有用的对照是保持质量门槛、基因集合、维数及其他参数一致，仅改变抽样种子重跑。给两次结果按标志组合命名，再按条码计算两次都被判到候选群的细胞数，以及只有其中一次入群的细胞数。如果两次都能识别同一核心群，一圈低检测深度的细胞反复跨界，报告时可以说明群核心的证据与边缘标签的敏感性。若整个小群只出现在一个种子里，下一步应回到原始表达和样本组成，在给它起生物学名字之前先厘清它出现的条件。

下载代码和示例数据，在 Linux 上复现

下载本篇 R 代码（ZIP，4.8 KiB） · 下载配套示例数据（ZIP，7.3 MiB）

将两个 ZIP 解压到同一位置，进入 001-seurat 目录，运行 Rscript --vanilla run.R。打开 results/metrics.tsv 查看结果，代码包中的 expected/metrics.tsv 可供逐行核对。

本次从下载包复现的步骤是“质控、抽样、聚类及全量标签投影”。Linux x86_64、R 4.5.2、Seurat 5.4.0 的输出为原始条码 2700，抽样细胞 800，质控后全量细胞 2638。

下载版另存 cluster-counts.csv。先分别把 sketch、full 两列求和，确认两张名单各自闭合，再用每群人数除以各自的列和。原文首次运行中，编号 3 为 94/800 和 205/2638。此次 Linux 运行中，同编号为 97/800 和 208/2638，分别为 12.125% 和 7.88%。跨运行核对细胞群要按条码与标志表达匹配，群号只在一次聚类内部有意义。

实跑后的经验总结

这次在实际 PBMC 文库中，群三在八百个建图细胞中占比高于它在全量投影细胞中的占比。复核时最有效的是先对齐条码，把谁参与了建图、谁收到了投影标签、谁进入了统计分母分别写清。迁移到多供体时，再把每人的总数和目标细胞数列出来。只要这三层名单还在，统计口径改变后就能重算，也能找出是哪些细胞改变了结果。

如果这套流程已经在小样本上跑通，下一步要把整批文库投影到同一张图，云生信也提供可自行运行 R 脚本的环境。了解云生信生物信息学分析服务器租赁