企业官网建设流程全解析

拿到cellranger输出的filtered矩阵那一刻，大多数新手的第一反应是：接下来呢？读进来、画两张图、然后呢？我见过太多人卡在这一步——矩阵读进来了，看着那个几万行乘几千列的表格，完全不知道哪些细胞该留、哪些该扔，更不知道后面那堆聚类分群代码跑出来到底在说什么。这篇就接着上一篇的进度，把scRNA-seq数据分析里最关键的“数据清洗”和“细胞群初步分型”这部分讲透。说白了，这一阶段做的所有事情，都是为了让下游结论站得住脚：如果输入的是脏数据，后面UMAP上任何一个“分群”都可能是假象。这篇文章适合刚跑完上游比对、准备自己完成“从表达矩阵到细胞类型注释”这条完整链路的人。我会按实际操作的顺序，把每一步为什么要这么做、参数怎么定、踩过哪些坑，一次说清。

1. 从计数矩阵到分析就绪：清洗前先看懂数据长什么样

1.1 上游产物到底给了你什么

不管用的是10x Genomics的Cell Ranger还是其他平台，最终拿到手的通常是一个基因×细胞的计数矩阵。以Cell Ranger的filtered_feature_bc_matrix为例，它已经帮你过滤掉了一部分明显空液滴，但“过滤掉空液滴”和“数据干净”之间还差着十万八千里。这个矩阵里的每个数字代表一个UMI（Unique Molecular Identifier）计数——就是这个细胞里，某个基因被捕捉并测序到了几次。

读取时用Seurat的话，代码不复杂：

library(Seurat) library(dplyr) library(Matrix) data_dir <- "path/to/filtered_feature_bc_matrix" counts <- Read10X(data.dir = data_dir) # barcode是列名（细胞），基因是行名 obj <- CreateSeuratObject(counts = counts, project = "my_sample", min.cells = 3, min.features = 200)

这里有两个容易忽略的参数：min.cells = 3表示一个基因至少在3个细胞里检测到才保留，目的是去掉那种只在单个细胞里零星出现的背景噪声（也可能是测序错误）；min.features = 200则是最起码有200个基因表达的细胞才保留，这一步会直接扔掉一些明显异常的空泡。这两个值不是死规矩，但作为起步标准非常稳。我个人的习惯是往后看分布再决定要不要收紧，而不是一开始就卡死。

1.2 三个质控指标：UMI数、基因数、线粒体比例

读进来之后，下一步就是看每个细胞的质量。Seurat会自动帮你算nCount_RNA（这个细胞所有基因UMI的总和）和nFeature_RNA（这个细胞检测到了多少个基因），但线粒体基因比例得自己算：

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

注意这个正则表达式是有物种差异的：人用^MT-，小鼠用^mt-。大小写搞错了，比例算出来全是0，我身边至少两个人在这上面浪费过一下午。

这三个指标背后都是有生物学意义的，不是随便挑出来的统计量：

• nCount_RNA反映测序深度和RNA捕获效率。如果某个细胞的UMI总数低得离谱，它很可能是空液滴或细胞已经破裂；如果高得反常，就要怀疑是不是两个细胞被当成一个测了。
• nFeature_RNA反映细胞转录组的复杂程度。测序深度低会导致基因检出数低，但反过来，一个细胞如果检出的基因数异常高，也可能是双细胞。
• percent.mt是最经典的质量标签。线粒体基因在活细胞中通常占一小部分（一般是5%~15%），当细胞膜受损、胞浆RNA泄漏后，细胞质里的mRNA丢失，剩下的主要是线粒体RNA，所以这个比例会飙升。看到percent.mt超过20%的细胞，基本可以直接视为垂死状态。

1.3 阈值怎么定：画图比翻教程有用

新手最喜欢问“阈值到底填多少”，答案永远是“看你的数据分布”。直接用VlnPlot和FeatureScatter看全景：

VlnPlot(obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0.1) FeatureScatter(obj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

正常的样本里，大部分细胞的nFeature_RNA会集中在一个主峰附近，旁边拖一条低质量细胞的尾巴。我在多数人源组织样本里的经验是：nFeature_RNA下限设在200~300，上限设在2500~4000；percent.mt上限设在10%~20%。但如果你做的是肝脏、心脏这类代谢活跃的组织，线粒体比例天然偏高，一刀切15%可能会把一群本来健康的细胞全扔掉。这时候宁可把阈值放宽一点，等聚类之后再看哪个cluster的线粒体比例异常高、单独处理。

我自己的经验法则是：先画图，找分布的“膝盖点”（elbow point），而不是拍脑袋填一个数。过滤时也只做一次，别来回反复试阈值——你试一次阈值就会引入一次偏见，后面统计检验会失真。

2. 清洗不止过滤低质量细胞：双细胞、环境RNA和批次效应

2.1 双细胞：藏在“高质量”细胞里的内鬼

很多新手过滤完低质量细胞就觉得万事大吉了，但最危险的一类污染恰恰躲在高UMI的细胞里——双细胞，也就是一个液滴里包了两个细胞。从质控指标看，它可能完全“健康”：基因数高、UMI数高、线粒体比例也不超标。但实际上它是两个人（或者两个不同类型细胞）的转录组掺在一起，在聚类时会在两个真实群体之间制造一个虚假的过渡群。

双细胞的识别需要用专门的工具，我常用的是DoubletFinder。它的核心思路是：人为在真实数据里注入模拟的双细胞，然后用机器学习判断真实细胞更像哪一个“模拟双细胞”。实际操作里有两个坑：

第一个坑是参数pK的选择。新手直接跑默认参数经常报错或者效果很差，正确做法是先做一遍参数扫描：

library(DoubletFinder) # 先做标准归一化和降维（后面会讲） obj_scaled <- NormalizeData(obj) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA() # param sweep sweep.res <- paramSweep(obj_scaled, PCs = 1:20, sct = FALSE) sweep.stats <- summarizeSweep(sweep.res, GT = FALSE) bcmvn <- find.pK(sweep.stats) # 选择BCmetric最高的pK值 optimal_pk <- as.numeric(as.character(bcmvn$pK[which.max(bcmvn$BCmetric)]))

第二个坑是预期双细胞比率（nExp_poi）。10x官方有个大致估算：每1000个细胞大约有0.8%的双细胞率，但这只是文库制备时的物理参数，实际还会受上样量影响。对于3万个细胞左右的样本，我通常会按1%~2%估算。设得低了会漏掉，设得高了会把真实细胞误杀。如果你用的是更复杂的组织（比如肿瘤），双细胞率会更高，可以适当上调到2%~4%。

2.2 环境RNA：背景污染比你想的更普遍

环境RNA（ambient RNA）是液滴里游离的RNA片段，主要来自制备过程中细胞破裂释放的mRNA。它在数据里的表现是：几乎所有细胞都“低表达”很多不该表达的基因。如果你发现某个分化标志基因在所有cluster里都有1~2的log-normalized表达量，却没有任何一个cluster明显高表达它，大概率就是环境RNA污染。

处理环境RNA有专门的工具，比如SoupX和DecontX。但我给新手的建议是：先观察，别急着矫正。环境RNA的影响在差异分析阶段通常可以容忍，因为所有细胞都受到同等污染，组间比较时背景会互相抵消。只有当你要分析某个稀有群体的特异表达、或者文库的背景污染肉眼可见地严重时，再考虑用SoupX做减法。矫正不当反而会引入新偏差，比如把低表达基因的信息抹掉。

2.3 批次效应：合并样本前必须想清楚的问题

如果你的项目里有多于一个样本，逃不掉的话题就是批次效应。所谓批次效应，指的是技术差异（不同天建库、不同测序批次、操作人员不同）导致的系统性表达差异，它和生物学差异混在一起时，会让UMAP上相同类型的细胞因为“出自同一批次”而各自抱团。

判断是否存在批次效应最直观的方法是：跑完聚类后，用样本身份给细胞着色。如果UMAP上的每个cluster都混合着所有样本，那是理想状态；如果出现“这个cluster几乎全是样本A，那个cluster几乎全是样本B”，就要警惕。

处理批次效应有两个主流路线：

• Seurat的CCA整合：适合跨数据集整合和注释转移，它能找到不同批次间共享的细胞状态。Seurat v5的写法是IntegrateLayers，整合后生成一个新assay（比如integrated），下游PCA聚类都基于这个assay。
• Harmony：速度极快，适合大批量样本（几十个以上的10x文库）。我个人在大规模整合时更喜欢Harmony，因为CCA在样本数很多时计算量会爆炸，而Harmony的迭代矫正逻辑非常稳。

不管用哪个，都要记住：整合是在“规范化+高变基因”做完之后、降维之前做的。而且整合完别用原始RNA assay去跑PCA，要用整合后的assay——这是一条新手最容易踩的坑，跑完发现cluster完全跟着样本走，回头一看代码，原来是忘了切换assay。

3. 标准化、高变基因与降维：为聚类铺一条干净的路

3.1 归一化到底在做什么，为什么不直接比较原始计数

归一化这个问题，我见过太多人完全理解反了。原始UMI计数最大的问题是：不同细胞测序深度不同，有的细胞测到了3万UMI，有的只有3000，如果直接拿原始数字比较，测序深度高的细胞里所有基因都“看起来更高表达”。所以标准做法是把每个细胞的总UMI数拉平到同一个水平：

obj <- NormalizeData(obj, normalization.method = "LogNormalize", scale.factor = 10000)

这行代码做的事情是：每个基因的表达量除以该细胞总UMI数，再乘以10000，然后做log1p变换（也就是log(x+1)）。这个变换有两个目的：一是把数据从偏态分布拉成接近正态分布，方便后续假设检验；二是把动态范围压缩几个数量级，让低表达基因也能参与下游比较。这也是为什么我前面说环境RNA的背景表达“1~2”是log-normalized值——它对应的原始比例其实是0.01%的量级。

另一个选择是SCTransform（sct），它在归一化的同时还能建模UMI计数深度对表达的影响，选高变基因也比默认的vst方法更稳健。但它的计算量比LogNormalize大不少，而且在Seurat里一旦用了sct，后面的FindVariableFeatures就可以跳过，因为sct已经选过了变基因并写在模型里了。我的建议是：新手先老老实实用LogNormalize流程跑通整个pipeline，搞清楚每个对象里data、scale.data、counts三个槽分别是什么。等熟练了，再换成sct也不迟。特别提醒：sct不适合用来做细胞周期矫正，也不要在同一个数据集里混用两种流程对比结果。

3.2 高变基因、ScaleData和PCA之间的因果关系

归一化做完，表达矩阵里依然有几万个基因。这里有个核心问题：如果直接拿全基因去算细胞间距离，结果会被那些在不同细胞间变化不大、只在个别细胞里爆发的基因主导（比如线粒体基因），真实的细胞类型信号反而被淹没了。所以需要先筛出高变基因（HVGs）：

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

为什么是2000？这不是生物学定律，纯粹是计算精度和效率的折中。2000个高变基因通常已经能覆盖绝大部分决定细胞身份的信息，再用更多基因，边际收益很低，计算代价却直线上升。你可以扩展到3000或者5000，但很少看到有人用全基因组做下游聚类。

选完高变基因，接下来是ScaleData。这一步做了两件事：对每个基因做z-score标准化（减去均值除以标准差），并把数据存进scale.data槽。为什么要z-score？因为PCA对方差极其敏感，如果一个基因的表达量范围是0到5，另一个是0到1000，不归一化的话后者会主导主成分方向。这里还有一个隐藏问题：z-score会让所有基因的方差基本相同，低表达基因的噪声也随之被放大，所以ScaleData通常只对高变基因做，正好衔接上一步。

后面就是标准的降维三部曲：

obj <- RunPCA(obj, npcs = 30, verbose = FALSE) ElbowPlot(obj, ndims = 30)

RunPCA把高维（2000维）的表达谱压缩成几十个主成分。ElbowPlot画出来之后你会看到一条从高处快速下降然后趋于平坦的曲线，拐点之后的主成分解释的方差占比已经很低。通常我会选拐点附近再加两三个主成分，比如拐点在12~15之间，就用15或者20个PC。不是PC越多越好——加上那些接近噪声的维度，聚类反而会分裂出虚假群体。

4. 聚类与UMAP可视化：细胞群的初始分群

4.1 KNN图聚类与resolution参数的直觉

降维得到的几十个主成分并不是终点，它们是用来算细胞相似度的。FindNeighbors基于PC坐标构建KNN图——每个细胞只跟它最相似的若干细胞连边。FindClusters则在这个图上运行社区发现算法（Seurat默认是Louvain，后来版本也支持Leiden），把连接紧密的节点分成一个个社区，也就是我们说的cluster。

resolution这个参数是新手最容易纠结的，它的作用简单说就是控制分群粒度：数字越大，分的cluster越多越细。我用过的经验值供你参考：

• resolution = 0.1：会把整个数据集分成很少的几个大类，适合第一眼看全局。
• resolution = 0.5~0.8：这是做初步分型最常见的区间，一般能得到8~15个cluster，和细胞类型基本对应。
• resolution = 1.2~2.0：会把已识别的细胞类型进一步切成亚群，适合后续精细分析。

注意：FindClusters跑出来的cluster编号是随机的，不代表任何生物学顺序。cluster 0不等于“最重要的细胞”，它只是社区发现算法遍历时第一个碰到的社区。很多新手会因为“cluster 0”在UMAP上占据最大面积，就想当然认为它是某种主要细胞类型，这完全是误会。

4.2 UMAP和tSNE怎么选

聚类完成后的可视化主流是UMAP和tSNE。两者的区别，用大白话说：UMAP更在意“全局结构”，速度快，在大数据集上能保持不同细胞类型之间的相对距离关系；tSNE则更擅长把局部邻域关系展示得极度清晰，但压缩全局结构、组间距离几乎没有意义。

实际使用上，我通常用UMAP做总览和审阅分群，用tSNE出最终论文图（因为它看起来更“干净”，细胞簇边界更明显）。跑起来也很简单：

obj <- RunUMAP(obj, dims = 1:15) obj <- RunTSNE(obj, dims = 1:15) DimPlot(obj, reduction = "umap", group.by = "seurat_clusters", label = TRUE)

dims的参数要和前面选的PC数目保持一致，否则你等于换了输入数据在画图。这是个细节，但每次有人跟我抱怨UMAP图“怎么长得跟前一次完全不一样”，十有八九是RunPCA之后换过npcs，忘了同步RunUMAP里的dims。

4.3 聚类之后先别急着注释，先“质检”

有太多人跑完聚类就急着找marker、给细胞命名。但在我自己的流程里，聚类完成后的第一件事是检查分群质量。具体看三件事：

• 每个cluster的QC指标是否正常。如果某个cluster的percent.mt中位数明显高于其他cluster，它很可能是一群垂死细胞，应该标记出来而不是硬着头皮注释。
• 每个cluster是否被某个样本主导。画出DimPlot(obj, group.by = "orig.ident")，如果看到颜色按cluster整片划分，说明批次效应没处理好，回头处理批次，而不是继续往下走。
• 有没有“小得可疑”的cluster。UMAP上偶尔会出现那种只有几个细胞的孤立小团，它们通常是低质量细胞或技术噪声，可以人工剔除，也可以暂时保留，但别指望能注释出什么可靠结论。

这个“聚类后质检”的习惯，能帮你省下后面至少一半的返工时间。

5. marker基因识别与细胞群注释：给cluster贴上身份标签

5.1 FindAllMarkers的正确打开方式：三看原则

注释细胞群的第一步是找出每个cluster的特征基因，Seurat的标准接口是FindAllMarkers：

obj_markers <- FindAllMarkers(obj, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25, test.use = "wilcox")

这里的关键参数逐一拆解：

• only.pos = TRUE：只关注在本cluster中上调的基因，因为负marker（某基因在其他cluster高表达）对注释的指示意义弱，而且会让输出量大增。
• min.pct = 0.25：基因至少要在cluster内和cluster外的25%细胞中有检测。这是一个“表达广度”的过滤条件，防止某个基因只靠几个极值细胞撑起一张显著面孔。实际项目里0.25偏宽松，我在粗略初筛时用0.25，筛选候选marker时会提高到0.4~0.5。
• logfc.threshold = 0.25：log2倍数变化的门槛。注意这个值是log2尺度，0.25对应的其实只有约1.19倍的表达差异，也就是个非常宽松的起步值。真正有价值的marker，我一般会在后续人工筛选中要求avg_log2FC > 1甚至更高。
• test.use = "wilcox"：Wilcoxon秩和检验。它是非参数检验，不假设数据正态分布，而单细胞数据普遍是零膨胀的，用t检验很容易被少数细胞带偏。默认用Wilcoxon是合理的选择，也是我在真实项目中的选择。

输出结果里有一列avg_log2FC和一列p_val_adj。很多新手只看p值，觉得p越小基因越“标志性”。但在单细胞数据里，p值受细胞数影响极大——一个cluster有5000个细胞，另一个有50个细胞，两者比较时本来就容易“显著”。我给自己定的判断标准是“三看”：一看p_val_adj是否小于0.05（基本门槛）；二看pct.1与pct.2的差距是否够大（表达覆盖度差异，例如pct.1=0.9而pct.2=0.1，才是真正的cluster特异性，而pct.1=0.9、pct.2=0.85说明这个基因只是普遍表达高，算不上特征）；三看avg_log2FC是否够大（表达量差异，至少大于0.5，严格时用1）。

# 查看cluster 0的top markers top0 <- obj_markers %>% filter(cluster == 0) %>% arrange(desc(avg_log2FC)) %>% head(10)

得到候选marker后，不要直接照着一篇文献就把细胞类型定死。打开一个已知的marker基因列表，去FeaturePlot上视觉确认一遍：高表达的细胞是不是正好落在这个cluster里？表达的细胞范围是否和UMAP上的分群轮廓吻合？这一步虽然“土”，但永远是最可靠的验证。

5.2 常见组织里值得优先查的marker面板

不同组织的marker体系差异很大，但有几个通行的“默认面板”覆盖了最常见的细胞类型。以人和小鼠的血液、肿瘤组织为例，我做了一张速查表：

这里要特别提醒一句：没有哪个基因是100%专一的。比如CD14在巨噬细胞里也很高，NKG7在部分T细胞亚群中也会表达。注释时永远要组合判断，至少两三个marker共同支持，再给结论。单一marker就断言细胞类型，基本是注释事故的前兆。

5.3 人工注释与自动注释的配合

现在也有不少自动注释工具，比如SingleR、Garnett，以及基于参考图谱的映射工具。对新手来说，我建议把这些工具当作“预注释助手”而不是“最终裁判”。以SingleR为例：

library(SingleR) library(celldex) ref <- HumanPrimaryCellAtlasData() singler_pred <- SingleR(test = GetAssayData(obj, assay = "RNA", slot = "data"), ref = ref, labels = ref$label.main) obj$SingleR_classification <- singler_pred$labels

跑完之后，把自动注释结果对照DimPlot看一遍：如果某个cluster里绝大多数细胞被标成同一种类型，那这是个很好的起点；如果同一个cluster被拆得七零八落、标出了五六种完全不同的类型，就要怀疑它根本是个混合群体或双细胞残留。

即使自动注释结果看起来很合理，我也会做两件事做最终确认：一是对每个cluster取出top markers，人工和已知marker面板比对；二是画一个DoHeatmap，把每个cluster的top marker表达量热图展示出来，看看是否“各群界限分明”。格式如下：

top_markers <- obj_markers %>% group_by(cluster) %>% top_n(5, avg_log2FC) DoHeatmap(obj, features = unique(top_markers$gene), group.by = "seurat_clusters")

如果热图上有大片横向条纹——也就是某个marker不仅在目标cluster高表达，其他cluster也一片红——那说明这个cluster的marker特异性不够，注释前要重新审视。

5.4 初步分型不等于最终结论

当所有cluster都有了细胞类型标签之后，通常会把细粒度的cluster合并成生物学上更粗的类别。例如cluster 1、5、9都是T细胞亚群，可以合并成一个T_cell大群。合并之后再回看UMAP，确认每个大群之间边界清晰。这一步其实是在校验“初步分型”的稳定性。我遇到过一种情况：某个cluster在第一次聚类时被分成了两半，合并后重新跑UMAP，发现两半其实分散到不同区域了——这说明之前的切分依据并不是真正的生物学差异，可能只是某个混杂因素（比如细胞周期）在作祟。

另外，初步分型阶段千万不要过度解读。测序深度不均、某些稀有细胞类型本来就只有几十个细胞，它们的注释置信度天然就低。对这类“边缘群体”，我通常在报告里标注为“疑似XX细胞”或“状态待定”，而不是硬给它一个确定的名字。数据清洗和细胞分型本身就是一环套一环渐进的，一次流程跑完能锁定主要细胞群，已经是很大的成功了。

最后分享一个小经验：注释是一个“越做越觉得自己之前是错的”的活。我回头看自己在同一个数据集上三个月前做的注释，总能发现几个当时误定的群体。这不是坏事，反而说明你在进步。分析单细胞数据，最忌讳的就是“得出一个漂亮的结论就停手”，多跑几个分辨率、多试几套参考注释、多画几张FeaturePlot，结论才立得住。建议在这一步把每个版本的cluster标记和marker列表保存成CSV，给打开的每个RDS命名都带上日期。不然等你三个月后再回来看这个项目，数据集还在，但当时的判断逻辑早就忘了。