单细胞转录组入门：PBMC数据从读取到聚类全流程解析

先说句实在话：单细胞转录组分析入门，最怕的不是代码报错，而是不知道自己每一步在干什么、为什么这么干。网上教程一大堆，但大多是“跟着抄一遍就完事”，换个数据、换个样本量就懵了。我当年从Bulk RNA-seq转到单细胞，第一次跑PBMC数据时也被一堆概念砸得晕头转向：标准化和归一化到底是不是一回事？PCA选多少个主成分才算合理？UMAP图上细胞聚在一起，为什么不能直接说它们是同一类？

这篇教程我就用经典的PBMC数据集，把从原始数据读取到聚类分群的完整流程走一遍。每一步除了给代码，还会尽量讲清楚背后的原理和判断逻辑，以及我实际跑数据时踩过、也看别人踩过的坑。内容和Seurat官方PBMC教程的主线一致，但会更偏向“实操中真正需要留神的地方”。如果你是刚入门单细胞分析、手上还没有自己的数据，或者跑了几个教程但知其然不知其所以然，这篇应该能帮上忙。

1. 为什么入门首选PBMC：这个数据集好在哪里

PBMC的全称是外周血单个核细胞，简单说就是血液里除去红细胞、血小板之后，剩下那些带细胞核的细胞群体，主要包括T细胞、B细胞、NK细胞、单核细胞和树突状细胞。做单细胞测序，它几乎是所有新手的第一个练手数据集，原因很朴素：细胞类型清楚、标记基因成熟、数据量适中，跑起来不费劲。

这套数据里的细胞类型是“已知中的已知”，这非常重要。你跑完聚类之后，拿已知标记基因去验证，如果每一群都能和预期对上，说明整个分析流程是通的。如果连PBMC都对不上，那多半是代码、参数或者数据读取出了问题，得回头检查。

从实际运行角度说，PBMC 3k大约有2700个细胞、2万个基因左右，一台普通笔记本8G内存就能跑完完整流程，几分钟到十几分钟不等。相比动辄几万、十几万细胞的数据集，试错成本极低。我建议所有刚接触单细胞的人，不管最终课题是不是做免疫，都先把这个数据集完整跑通一遍，跑通了再碰自己的数据，心里才有底。

需要明确的是，这篇教程主线是在Seurat框架下走一遍标准流程。Seurat是目前单细胞分析里生态最成熟、资料最多的R包，社区讨论量大，遇到报错基本都能搜到解决方案。后面所有代码默认你用的是4.x或5.x版本，个别函数在新旧版本里有差异，我会在对应位置提醒。

2. 环境准备：R版本、Seurat依赖和那些莫名其妙的加载冲突

先说版本。Seurat对R版本有要求，建议装R 4.2以上。装的时候别用太老的R，否则很多依赖包装不上。RStudio不是必需的，但强烈建议装一个，查看数据框、画图、看帮助文档都方便很多，尤其对新手来说，纯命令行操作容易劝退。

安装Seurat本身不复杂：

install.packages("Seurat")

但国内网络环境下经常下载到一半就失败，或者卡在某个依赖包上。直接设置国内镜像可以省掉很多烦恼：

options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) install.packages("Seurat")

如果是旧项目或服务器上已经装了老版本，升级之前最好先看一眼当前环境里的版本，有些代码在版本之间是不兼容的：

packageVersion("Seurat")

自己学习用的话，我建议直接装最新稳定版，不要用开发版。开发版虽然新功能多，但文档不齐全，报错也没处查。

还有一个新手几乎都会遇到的坑：加载Seurat之后，再加载dplyr，会提示大量函数被“masked”，比如filter、select都会被dplyr覆盖。如果你习惯了用dplyr的过滤逻辑，后面很容易写出错误的代码。我自己的习惯是，分析脚本里只加载Seurat和ggplot2、patchwork这些可视化包，需要用数据处理的函数时用dplyr::filter()这种带包名的方式显式调用，避免命名空间冲突。

内存方面，跑PBMC 3k这种规模，8G内存完全没问题。但如果你后面要跑几万细胞的数据，建议提前了解一个参数：options(future.globals.maxSize = 4 * 1024^3)，这个可以把后续并行计算允许传输的对象上限调大，否则某些步骤会报“size limit exceeded”的错误。这个我现在几乎每个大项目都会先设一下。

3. 数据读取：目录结构就是一个大坑

Seurat官方教程用的PBMC 3k数据来自10X Genomics官网，是标准的Cell Ranger输出格式。下载解压后你会看到三个核心文件：

• barcodes.tsv.gz：细胞条形码列表，每一行代表一个细胞。
• features.tsv.gz：基因列表，新版是features，老版本叫genes。
• matrix.mtx.gz：表达矩阵，记录每个细胞中每个基因的UMI计数。

这三个文件必须放在同一个目录下，而且文件名不能随意改动。用Read10X函数读取时，你传入的是这个目录的路径，而不是单个文件名。这是很多新手第一次读取时报错的原因：搞不清楚到底是传文件夹还是传文件。

library(Seurat) # 假设数据解压到了工作目录下的 pbmc3k 文件夹中 pbmc.data <- Read10X(data.dir = "pbmc3k/")

读进来的pbmc.data是一个稀疏矩阵，行是基因，列是细胞。可以用dim()看一下维度：

dim(pbmc.data)

正常情况下你会看到类似32738 x 2700的结构，2700就是细胞数。看到这个数字和官方教程对得上，说明数据没读错。

还有一个关键步骤，是把稀疏矩阵转成Seurat对象：

pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)

这里面的min.cells = 3和min.features = 200是第一个需要理解的过滤参数。

min.cells = 3表示一个基因至少在3个细胞里有表达才会被保留。原因是单细胞数据噪音很大，如果一个基因只在极少数细胞里检出，大概率不是真实生物学信号，而是测序或建库过程中的背景污染。保留它们只会增加后面标准化和PCA的计算噪音。

min.features = 200表示一个细胞至少检测到200个基因才会被保留。这个阈值主要过滤掉空液滴或破碎细胞。什么是空液滴？10X平台是让单个细胞和带有条形码的微珠一起进入油滴，但总有一些液滴里没有细胞，或者细胞已经破裂，这些液滴里能测到的基因数很少。设置一个最低基因数，能把这些垃圾“液滴”去掉。

拿到Seurat对象后，可以看一个非常基础但很实用的信息：pbmc这个对象本身会显示有多少样本、多少细胞、多少基因。关于这部分，我想多说一句：很多人拿到数据后急着往下跑，其实应该先看一眼这些基础数字，它们往往是后续所有判断的起点。

4. 质控不只是一步代码：nFeature、nCount和线粒体比例到底在过滤什么

下一步是质控（Quality Control, QC）。这一步是单细胞分析里最容易被糊弄过去的，因为代码就那么几行，但背后涉及的问题很多。

先搞清楚QC在看什么指标：

• nFeature_RNA：每个细胞中检测到的基因数。基因数过低，说明细胞可能破裂或空液滴；过高，则可能是双细胞（两个细胞被当成一个测了）。
• nCount_RNA：每个细胞中检测到的总UMI数（可以近似理解为总转录本数）。这个数和基因数正相关，测序深度越深，总UMI数越高。
• percent.mt：线粒体基因的UMI占总UMI的比例。线粒体基因比例高，通常说明细胞状态不好——细胞破裂后，胞浆RNA大量流失，而线粒体相对完整，所以线粒体基因占比会异常升高。

计算线粒体比例用Seurat的现成函数：

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

注意这个正则表达式是区分物种的。人类是^MT-，小鼠是^mt-。如果你用的是其他物种，要先去确认线粒体基因的命名规则，否则算出来全是什么也不像的结果，等于没过滤。

然后可以画一个小提琴图看看分布：

VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

这一步不要直接跳过画图就开始过滤。看分布的目的是确定合适的阈值。很多教程上来就写nFeature > 200 & nFeature < 2500 & percent.mt < 5，但这个阈值真的适合所有数据吗？

答案是未必。

线粒体比例5%这个阈值，更多是细胞系或健康外周血样本的经验值。不同组织、不同处理条件下的样本差异很大，比如肝脏或心肌这类线粒体丰富的组织，5%基本会把大多数细胞过滤掉。正确做法是先看小提琴图的分布形态，找到明显偏离主群的那个“尾巴”，在尾巴位置切一刀。用Seurat官方教程的话说，这不只是选阈值，而是在识别异常细胞群。

常规过滤代码是：

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

但我要强调一点：这个过滤只是粗筛，真正的“双细胞”问题靠这个阈值是处理不彻底的。nFeature特别高的细胞不一定是双细胞，可能是分裂期细胞或者代谢特别旺盛的细胞。如果项目对细胞纯度要求很高，后面需要额外的双细胞检测工具（比如DoubletFinder），尤其是10X上机时细胞浓度偏高的情况下。这个话题后面有机会单独写，这里先留个印象。

过滤完之后，我习惯再画一次小提琴图，看看过滤后的分布是不是干净了很多。这一步能帮你确认阈值设得合不合理，而不是无脑套用教程参数。

5. 标准化和高变基因：为什么默认选2000个而不是越多越好

过滤完细胞和基因，接下来是单细胞分析中最重要的数据变换任务。这一步包括NormalizeData()和FindVariableFeatures()，很多人会把它们笼统地叫“标准化”，但实际上它们做的是不同层级的事情。

NormalizeData默认方法是LogNormalize，公式本质上是：

[ log1p(\frac{counts}{totalCounts} \times 10000) ]

翻译成人话就是：先求每个基因的count数占该细胞总count数的比例，然后乘以10000这个缩放因子，最后取log1p（也就是ln(1+x)）。为什么要做这个变换？因为不同细胞的测序深度差异很大，一个细胞测了5万UMI，另一个只测了1万UMI，直接比较原始count数没有意义。归一化之后，数值代表的是“相对表达量”，可以跨细胞比较了。

pbmc <- NormalizeData(pbmc)

紧接着是找高变基因：

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

默认选2000个基因。为什么是2000，不是1000也不是5000？

高变基因就是那些在不同细胞间表达量波动最大的基因。这些基因承载了区分细胞类型的主要信息。选太多（比如全部基因），会把大量表达量平稳的“管家基因”也纳入后续计算，这些基因对区分细胞类型没有贡献，只会增加噪音和计算量。选太少，则可能丢掉一些低表达但关键的标记基因。2000个是实践中验证过的平衡点：信息量足够，计算量可控。

选完之后可以用一个可视化检查基因选择效果：

top10 <- head(VariableFeatures(pbmc), 10) plot1 <- VariableFeaturePlot(pbmc) plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE) plot1 + plot2

如果点的形状呈现标准的“火山”形态，左上角那些高方差基因被标出来，说明选择过程正常。

接下来还有一个ScaleData步骤，这个很容易被误解。它的作用是中心化和标准化：让每个基因在所有细胞中的表达均值为0、方差为1，也就是做Z-score变换。这一步是必须的吗？如果后续要做PCA，是必须的，否则表达量高的基因会主导主成分计算，而ScaleData让每个基因处于同一量纲上。

pbmc <- ScaleData(pbmc)

注意，默认只对前面选出的2000个高变基因做缩放，而不是所有基因，原因就是省内存。对一个几万基因、几万细胞的数据集，如果对所有基因做ScaleData，会生成一个极其巨大的稠密矩阵，内存直接爆掉。所以不要没事找事去对全基因做ScaleData。

6. PCA降维不是玄学：主成分数量应该这样选

ScaleData之后，进入降维环节。PBMC数据集有2万多个基因，如果直接在高维空间里算细胞之间的距离，不仅计算量巨大，而且高维空间里绝大多数维度都是噪音。PCA的作用就是把这些基因表达信息压缩成少数几个主成分，每个主成分都是原始基因表达量的线性组合，按解释方差的大小排序。

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

跑完之后，可以看每个主成分的基因载荷，也就是哪些基因对PC1、PC2贡献最大：

VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

还可以把细胞投影到前两个主成分上看一眼分布：

DimPlot(pbmc, reduction = "pca")

PCA跑完之后，接下来最重要的决定是：到底选多少个PC用于后续聚类？

这个问题没有标准答案，但有几个判断套路。

第一个方法是看碎石图：

ElbowPlot(pbmc)

图上每个点代表一个主成分能解释的标准差，随着主成分数量增加，解释力度快速下降，然后趋于平缓。那个“拐点”之后的主成分基本就没什么信号了。PBMC 3k数据跑出来，通常拐点都在10~20之间。

第二个方法是官方教程里用过的JackStraw，本质是通过随机置换检验每个PC的显著性。但这个方法计算量大，而且对大数据集来说非常慢，我自己在实际项目中基本不用，ElbowPlot足够。

第三个方法是最朴素的：直接对比不同PC数量下的聚类结果。比如分别用dims = 1:10和dims = 1:30跑聚类，然后用已知标记基因看哪组结果更干净。这个方法虽然有点笨，但最靠谱。

我个人的习惯是：先选20个PC跑第一遍，画出UMAP看整体分群情况，再用标记基因检查。如果分群结果和已知细胞类型对不上，再回去调整PC数量。比起一开始就在理论上纠结选多少个PC，这种快速试错的方式效率更高。说到底，PC数量不是一锤子买卖，后面的聚类步骤完全可以反复调整。

7. 聚类的坑：UMAP图看着分开了，不代表聚类就成功了

聚类这一步，很多人拿到UMAP图觉得“分得挺开”，就以为分析成功了。但实际上，UMAP只是把高维空间的关系投射到二维平面上的可视化结果，它不是聚类的依据，更不是定量分析的依据。

我见过一个很典型的误区：有人直接对标准化后的表达矩阵跑UMAP，然后用UMAP上的坐标做聚类。这个做法在概念上就是错的。Seurat的正规做法分两步：

第一步，构建细胞间的关系图：

pbmc <- FindNeighbors(pbmc, dims = 1:20)

这里用前面选好的PC维度，在高维PCA空间里计算细胞之间的距离，并为每个细胞找到它的近邻，构建一个K近邻图。

第二步，用图聚类算法把细胞分成群：

pbmc <- FindClusters(pbmc, resolution = 0.8)

默认用的是Louvain算法，新版Seurat 5里还可以选Leiden算法。这两个算法区别在于社区检测方式不同，Leiden在理论上能避免一些Louvain的“坏连接”问题，实际使用中结果差异不会太大。如果要指定算法，在FindClusters里传algorithm = 4（Leiden）或者直接让它用默认。

这里的核心参数是resolution，可以理解为“分辨率”：数值越大，分出的群体数越多。0.5会分出比较粗的群体，1.0以上会分出更多细小的亚群。官方PBMC教程用的0.5，能分出9个左右群体，对应T细胞、B细胞、NK、单核、树突细胞等主要类群。

当你调大分辨率时，某些大群会被继续细分。比如CD8+ T细胞可能被分成效应记忆亚群、中央记忆亚群，CD4+ T可能分出初始、记忆、调节等。新手看到一群CD4+T被分出好几个小群，第一反应往往是“算法有问题”。但其实这恰恰是单细胞分析的正常状态——细胞类型之下还有细胞状态，分辨率越高，看到的细节越多。

问题是：多高的分辨率才是“对的”？

我的标准是三分法结合：A. 分出的群体能不能用已知标记基因解释；B. 每个群体的细胞数是不是少到统计上没意义（比如小于总细胞数的1%）；C. 调高或调低一个档位，关键细胞群是否依然稳定存在。如果三个条件都满足，那这个分辨率就是可用的，而不是纠结于“0.8比0.6更正确”。

聚类完成后，画UMAP：

pbmc <- RunUMAP(pbmc, dims = 1:20) DimPlot(pbmc, reduction = "umap", label = TRUE)

顺便提醒一点：RunUMAP的dims参数要和FindNeighbors保持一致。有人前面用1:15聚类，后面画图却用了1:30，结果图和cluster对应不上，又找不到原因。这个低级错误我见过不止一次。

每次跑完聚类和UMAP，我强烈建议保存一下中间结果：

saveRDS(pbmc, "pbmc_seurat.rds")

后面不管是继续分析还是重新调整参数，直接从RDS读回对象就行，不用从头再跑一遍。这个习惯能省下大量时间。

8. 分群之后怎么验证：标记基因检查和手动注释入门

聚类完成只是“数学上把细胞分开了”，不等于“知道每一群是什么细胞”。把cluster编号对应到具体细胞类型，是单细胞分析里最需要生物学功底、也最考验耐心的一步。这也正是“手动注释”最近讨论度很高的原因——自动化注释工具虽然快，但经常混淆相似细胞群，最终判断还是要靠人。

手动注释的基本思路说起来很简单：每个cluster都有相对特异的标记基因，把这些标记基因的表达情况拉出来看看，就能推断它是什么细胞。但你不可能一个一个基因去查图表，效率太低。标准做法是先用FindAllMarkers把每个cluster的高表达差异基因找出来：

markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

这个函数会比较每个cluster和所有其他细胞，找出在这个cluster里显著高表达的基因。得到结果后，按cluster分组、按avg_log2FC排序，看每个cluster的top基因。然后把这些top基因和已知细胞类型的标记基因做匹配。

PBMC数据里最常见的几类免疫细胞标记如下：

注意几个容易踩的坑。第一个是不要用单个基因定义细胞类型。NK细胞和CD8+ T细胞都表达NKG7，单看这一个基因分不清；但结合CD8A和CD3D来看，NK细胞不表达CD3D，而CD8+ T细胞表达。第二个是要注意某些基因在不同细胞群都有低水平表达，比如CD14在单核细胞里高表达，但在树突状细胞里也可能有信号，所以要看相对表达量，而不是“有或没有”。

验证标记基因最常用的三种可视化方式：

# 1. 小提琴图：看基因在每个cluster的分布 VlnPlot(pbmc, features = c("CD3D", "MS4A1", "NKG7")) # 2. FeaturePlot：把表达量映射到UMAP上 FeaturePlot(pbmc, features = c("CD3D", "CD14", "MS4A1", "NKG7")) # 3. 热图：每个cluster的top marker基因整体展示 top10_markers <- markers %>% group_by(cluster) %>% slice_max(n = 10, order_by = avg_log2FC) DoHeatmap(pbmc, features = top10_markers$gene)

我个人做注释的时候，通常会先把FindAllMarkers的结果按cluster整理出来，每个cluster看15~20个top基因，自己先形成一个初步判断，再用上面的图去横向验证。这个过程比想象中花时间，但值得做。最近很多人在讨论自动注释工具，比如SingleR、Celldex、Garnett，也有基于数据库的打分注释。我的观点是：自动注释可以当参考，但不能当结论，它输出的标签最终还是要拿已知标记基因去人工复核。原因很简单，自动注释依赖参考数据集，参考数据集本身可能不包含某些稀有亚群，或者不同平台之间的批次效应会干扰判断。

手动注释还有一个进阶做法，就是给每个cluster算一个Marker Score。简单点说，把你认为属于某细胞类型的几个标记基因的表达量取出来，求平均，然后看哪个cluster的平均分最高。这个方法在判断单核细胞和树突状细胞这种“你中有我”的相似群体时特别有用。我在这篇先不展开完整算法，给个思路：Seurat的AddModuleScore函数就是干这个的，你可以用它给每个cluster打一个“T细胞分数”或者“单核细胞分数”，再结合聚类结果做判断。

9. 我替你们踩过的那些坑：配置、函数和批次效应的实操细节

最后这部分，我把自己跑PBMC流程时真实遇到过的坑集中写出来。很多问题单独看很小，但能卡住一晚上。

第一个坑是Read10X报错找不到genes.tsv。旧版Cell Ranger输出的是genes.tsv.gz，新版是features.tsv.gz。如果你的数据是老版本，而Seurat新版本默认找features.tsv，就会报错。解决方法是手动创建features.tsv.gz，把genes.tsv.gz复制一份改名即可。或者用Read10X的替代方式去读原始矩阵。

第二个坑是dplyr::filter和Seurat的subset混淆。我在前面说过，加载dplyr会mask掉很多函数。如果你写了filter(percent.mt < 5)但没指定包名，而当前环境中没有加载dplyr，就会报“找不到filter函数”。反过来，如果你加载了dplyr但用的是Seurat对象的列名过滤，也可能得到奇怪的结果。我的建议是：操作Seurat对象一律用subset()，操作普通数据框可以放心用dplyr。

第三个坑是线粒体基因匹配不上。PercentageFeatureSet的正则表达式必须和物种对应，而且要注意基因symbol在矩阵里的写法。有些公司下线的数据用Ensembl ID或者加了版本号，这时候^MT-根本匹配不到任何基因，percent.mt全为0，质控等于白做。处理前先确认基因名格式，如果发现是Ensembl ID，需要先转换成symbol再做这一步。

第四个坑是DoHeatmap画不出来，或者画出来一大片都是同一个颜色。这通常是因为用了全部基因做ScaleData时的内存问题，或者features里传入的基因名在矩阵里不存在。先用%in%确认基因名列表都在rownames(pbmc)里，再去画图。

第五个坑是高估UMAP的“距离意义”。很多人看到UMAP图上一个群体离另一个远，就说这两个群体“差异大”。UMAP为了在二维平面上展示高维结构，会刻意拉伸和压缩空间，图上的距离并不能直接反映转录组的差异程度。要比较群体间差异，应该用FindMarkers去看差异基因数量，或者算群体间的平均距离。这个观念不纠正，后面分析容易跑偏。

第六个坑是后面做整合分析时的批次效应。这篇教程只分析了单个样本，但实际项目里很少只跑一个样本。如果你有多个样本，不要直接把所有细胞合并起来跑上面的流程，那样批次效应会严重干扰分群。规范做法是用SCTransform加Harmony或者Seurat的IntegrateData做批次校正。这个话题涉及的内容比较多，我计划放在这个系列的后续文章里单独展开。

最后再分享一个我自己的习惯：整个分析过程中，每完成一个关键步骤就saveRDS一次。数据读取后存一份，QC之后存一份，PCA之后存一份，聚类注释完再存一份。因为单细胞分析经常要回头调整参数——可能是QC阈值要改，可能是PC数量要重新选，也可能是注释到一半发现某个群需要更高分辨率细分。如果不存中间结果，每一次调整都要从头跑一遍，几十万细胞的数据集跑下来是很痛的。

这篇先把PBMC的完整流程走通了。下一步其实更关键，也更有意思：怎么从这些cluster里鉴定标记基因、怎么做更严谨的手动注释、多组样本之间怎么消除批次效应做整合分析。那些内容我都会在这个系列里继续展开，讲清楚每一步的判断依据。