单细胞RNA-seq分析这条流水线，越往后越容易让人迷失。聚类那几步跑完，UMAP图看起来漂漂亮亮，但一追问“这几群细胞到底差在哪里”，就卡住了。NBIS的单细胞系列教程我一路跟到第五篇，终于走到差异基因（DEG）这一环，才发现前面做的所有降维、聚类、注释，最后都要落到这张基因列表上。这篇内容适合正在跑单细胞数据分析、尤其是已经拿到聚类结果但不知道下一步怎么处理的人，也适合那些用Seurat做了初步分析、却在“找差异基因”和“把结果解释清楚”之间反复横跳的朋友。我会把整套差异基因分析拆开讲清楚，从方法选型到手动注释，再到环形热图这种比较进阶的可视化，全部串起来，给你一条可以直接照着走的路径。
1. 从聚类到差异基因：第四篇结束之后，这一篇到底在算什么
1.1 前面四篇做了什么铺垫
NBIS这套单细胞教程的前面几篇，基本把标准流程走了一遍：从原始测序数据质控、细胞过滤、双细胞去除，到归一化、找高变基因、PCA降维、UMAP/tSNE可视化，再到基于图的聚类。到第四篇结束时，每个细胞都已经被打上了cluster标签，UMAP上能看到一个个分群，但这些群的身份还是未知的。
大多数人在这个节点会干一件事：打开Seurat的FindAllMarkers()，跑一遍，把padj小于0.05的基因全捞出来，然后对照文献找marker。这个流程没错，但很容易忽略一个问题——差异基因分析不是单纯跑一行代码，它背后牵扯到“你拿什么当分组”“你用什么检验”“你怎么看结果”这三个决策。
1.2 差异基因的三个比较逻辑
单细胞差异分析至少有三层不同的比较逻辑，混着用很容易得出自相矛盾的结果。
第一层是cluster之间的比较，这是最常用的。你想知道cluster 0和cluster 1在转录组层面有什么区别，于是把这两个cluster里的所有细胞作为两组，做差异检验。这种比较回答的问题是“这两群细胞的表达状态哪里不同”。
第二层是样本/个体之间的比较，比如对照组和处理组。这里要小心的是一般单细胞数据里有多个样本，每个样本里又有很多细胞，如果直接把细胞当重复，会把同一只小鼠的不同细胞当成独立重复，导致统计假阳性膨胀。
第三层是细胞类型内部的扰动比较，比如CD8 T细胞在肿瘤组和正常组之间的差异。这种分析通常先做细胞类型注释，然后指定某一种细胞类型，在不同条件下做差异分析。
NBIS教程的第五篇，核心其实是把第一层做扎实，并且带出手动注释的思路——找到差异基因不是终点，把这些基因和细胞身份对应起来才是。很多人在这一步翻车，是因为跳过注释直接拿cluster编号讲生物学故事，结果发现cluster 3既表达T细胞marker又表达髓系marker，根本没法解释。
2. 差异基因计算之前的三个准备动作，直接影响结果质量
2.1 分组设计：先问清楚“谁和谁比”
FindMarkers()在Seurat里最简单的用法就是指定ident.1和ident.2，默认分组是cluster。但实际操作里，你需要先确认几个前提。
前提一：cluster是稳定的吗？如果聚类分辨率换了一下，cluster就被拆开或合并了，那么基于这个cluster的差异基因结论是不稳的。我一般建议先跑几个分辨率（0.4、0.8、1.2），确认目标cluster在哪个分辨率下稳定存在，再做后续分析。
前提二：样本组成会不会干扰差异结果？假如cluster 0主要来自样本A，cluster 1主要来自样本B，那么cluster之间的差异基因可能反映的是样本效应，而不是细胞类型差异。这时要检查每个cluster里的样本构成比例，如果发现明显的样本偏倚，需要用FindMarkers()里的latent.vars参数（配合MAST方法）或者先做样本整合（Harmony、CCA等）来纠正。
2.2 归一化方式与高变基因：换了方法结果会漂
很多人在数据预处理时用的是默认的LogNormalize()，也就是Seurat最传统的标准化。但如果你前面用的是SCTransform，那么差异分析的输入数据格式就要对应调整——SCTransform之后，FindMarkers()默认会跑在SCT的assay上，表达值和对数倍数变化的计算会不一样。实测下来，同一条数据用LogNormalize和SCTransform跑同一组差异基因，拿到的top基因排名会有不小出入，尤其是高表达基因的排名变动很大。这不是谁对谁错，而是两种归一化对高表达基因的压缩方式不同。
高变基因的选择也会影响后续结果。如果你当初用FindVariableFeatures()选了2000个高变基因，那么FindMarkers()默认只在这2000个基因里跑差异检验——因为Seurat的默认features参数是NULL时，会跑在VariableFeatures()上。如果你只关心全部基因里的某个特定通路，一定要显式传入features，否则结果里根本不会出现非高变基因。
2.3 过滤环节的默认参数不能照抄
Seurat的FindMarkers()默认有个关键参数：min.pct = 0.1，意思是某个基因至少在两组中某一组的10%细胞里表达，才会纳入检验；另一个是logfc.threshold = 0.25，要求平均log2FC的绝对值大于0.25。这两个阈值的作用是把那些“只有极少数细胞表达、表达差异也不明显”的基因过滤掉，降低多重检验的负担。
但实际场景里，这两个默认值经常不合适。比如你研究的是稀有细胞亚群，一个cluster只有200个细胞，那么10%就是20个细胞，这样的min.pct设置勉强合理，但如果某类marker本身就是低频表达基因，比如转录因子，表达比例可能不足5%，默认参数会直接把这些基因过滤掉。所以我的习惯是先跑一遍全参数（min.pct = 0，logfc.threshold = 0），拿到全量结果后，再用自己的阈值做筛选，而不是依赖函数内部过滤。
3. FindMarkers之外：常用差异检验方法怎么选，参数怎么调
3.1 四种方法的适用场景
Seurat的FindMarkers()支持多种检验方法，常用的有这几种，我放在一起对比：
| 方法 | 适用场景 | 特点 | 注意事项 |
|---|---|---|---|
| Wilcoxon秩和检验 | cluster间快速筛选 | 速度快，不要求正态分布，对零膨胀有一定容忍度 | Seurat默认方法，但容易把“表达比例有差异但表达量差异小”的基因也算进去 |
| MAST | 考虑检测率（dropout）的影响 | 使用hurdle模型，把基因表达分为“是否检测到”和“检测到后的表达量”两部分 | 对细胞数量和计算资源要求更高，支持latent.vars矫正混杂因素 |
| DESeq2 | 样本层面的差异分析（伪bulk） | 需要原始counts，基于负二项分布，适合有生物学重复的实验设计 | 不建议直接用在单细胞层面，因为细胞之间不是独立样本 |
| edgeR/limma | 伪bulk分析的另一选择 | 速度和稳定性都不错，可以处理小样本量 | 也需要先聚合到样本/个体层面，再跑差异 |
如果你只是快速看看cluster之间有哪些候选基因，Wilcoxon足够。如果后续要发文章、需要更严格地控制假阳性，并且你的数据存在明显的批次效应或样本构成不均，建议至少用MAST跑一遍对比结果。如果研究对象是不同处理组之间的同类型细胞，那最优路径是PseudobulkExpression()聚合样本，再交给DESeq2。
3.2 FindMarkers实操中的几个关键参数
拿常见的二群比较举例，我的代码通常是这样的：
library(Seurat) # 假设seu已经完成聚类，celltype列是手动注释或cluster编号 Idents(seu) <- "seurat_clusters" # 比较cluster 0 和 cluster 1 degs <- FindMarkers( seu, ident.1 = "0", ident.2 = "1", test.use = "wilcox", min.pct = 0.1, logfc.threshold = 0.25, only.pos = FALSE )几个容易被忽略的参数：
only.pos = TRUE：如果你想找的上调marker，只返回阳性的差异基因，结果列表会更干净。但这个参数自身不加p值过滤，筛选逻辑完全靠min.pct和logfc.threshold。min.diff.pct：这个参数在Seurat 4之后的版本里可以设置，要求两组间的表达细胞比例差异达到某个阈值，能有效过滤掉那些两组表达细胞比例差不多、只是个别细胞表达量极高的“假差异”。做免疫细胞亚群分析时，这个参数特别好用。assay：确认你跑差异分析的那个assay是什么。如果做过SCTransform，却还在RNA这个assay上跑，结果会和SCT有差异。我遇到过有人导入旧代码，没指定assay，跑出来的基因列表和自己之前对不上，查了一下午才发现。
3.3 阈值怎么定，怎么防止“显著但没意义”
拿到差异基因列表之后，真正的筛选才刚刚开始。Seurat返回的结果里有三列最有用：p_val、avg_log2FC、p_val_adj。p_val_adj是校正后的p值，默认用Bonferroni方法，因为单细胞基因数量多，这个校正极其严格。你可以观察到有些基因原始p值可能是1e-30，但校正后变成0.05，就是因为基因总数有两万多个。
我的筛选标准通常是这样一套组合拳：
p_val_adj < 0.05，这是底线；|avg_log2FC| > 0.5，这个阈值比默认0.25更严，选出来的基因在后续实验里更容易被qPCR或免疫组化验证到；- 手动检查
pct.1和pct.2两列：如果一个基因在cluster 0的80%细胞里表达，在cluster 1只有5%，即使log2FC不算大，这个基因也值得关注，因为它反映的是“表达细胞比例”的差异，而不是单纯的平均表达量差异。
做到这一步，你会得到一个几十到几百个基因的列表。但列表本身不能回答“这群细胞是什么”，这时候就要进入手动注释环节。
4. 差异基因和细胞身份对不上？把手动注释这一步补齐
4.1 为什么做到差异基因这一步要回头做手动注释
自动注释工具（SingleR、CellTypist、Garnett等）已经很好用了，但它们在单细胞数据分析里的位置更像“预注释”，而不是终审。原因很简单：参考数据库不覆盖所有组织状态，尤其是疾病样本或少见细胞类型，自动注释经常会给出一个模糊的、甚至是错误的标签。
手动注释的逻辑是：拿差异基因列表和已知的细胞类型marker交叉验证。比如你是做肿瘤免疫的，cluster 5的top差异基因里有CD3D、CD3E、IL7R，那基本可以认定是T细胞；如果同时出现CD8A、GZMB、NKG7，可以进一步判断是细胞毒性T细胞或NK样T细胞。这种判断无法完全交给算法，因为不同文献里对同一群细胞的命名都不完全一致。
4.2 手动注释的完整流程
我做手动注释的习惯分四步：
第一步，拿到候选marker。用FindAllMarkers()输出每个cluster的top差异基因，通常每个cluster选top 20就够了。multi-marker验证比单个marker靠谱得多。
第二步，特征图叠加看表达模式。不要只盯着表格里的数字，用FeaturePlot()把关键marker打到UMAP上，看它的表达是否特异地集中在某一群细胞里。如果某个marker在UMAP上到处都亮，那它就不是好的细胞类型标志物，即使统计上显著也没用。
FeaturePlot(seu, features = c("CD3D", "CD8A", "GZMB", "NKG7"), cols = c("lightgrey", "firebrick"))第三步，做点图/气泡图。DotPlot()比FeaturePlot()更适合同时对比多个cluster的多个marker，因为它能直观显示每个cluster的阳性细胞比例和平均表达量。
DotPlot(seu, features = c("CD3D", "CD3E", "CD8A", "GZMB", "NKG7", "MS4A1", "CD79A", "LYZ", "FCGR3A"))第四步，结合文献给cluster重命名。这一步是整个单细胞分析里最依赖经验的地方。看到marker组合能对上某一类已知细胞，就把seu$celltype赋值上去，然后再跑一次差异基因，验证新分组下的差异基因是否符合预期。
4.3 手动注释和差异基因分析之间有个顺序陷阱
这里有个常见的路线问题：到底是先注释再跑差异，还是先跑差异再注释？
我的建议是：先跑cluster间的差异基因，拿到候选marker后再手动注释，然后基于注释结果再跑一轮细胞类型间的差异分析。第一次差异分析是为了注释服务的，第二次差异分析才是真正为了找“生物学意义”的差异基因。如果你一开始就把某个cluster命名为“T细胞”，又从T细胞marker里挑差异基因，这就是循环论证——你已经在marker里选了基因，然后又拿这些基因去证明这群细胞是T细胞。这个顺序问题，是我见过最多的一个方法论错误。
5. 环形热图：当常规热图放不下几十个marker时的另一种画法
5.1 环形热图和普通热图的差别
差异基因可视化最常见的是火山图和热图。火山图适合展示全部基因的分布，热图适合展示一组特定基因在不同cluster/细胞类型中的表达模式。但用普通热图会遇到一个问题：当marker基因数量多（比如50个）、细胞类型也多（比如10种以上）时，DoHeatmap()画出来的图会变得非常宽，基因名挤成一团，根本没法看。
环形热图（circular heatmap）把矩形热图的“行”映射到圆周上，基因沿着圆周排列，不同细胞类型用同心圆环表示。它的优势是信息密度高、结构紧凑，特别适合展示几十个marker在多种细胞类型中的表达模式，在文章里作为summary figure很受欢迎。
5.2 circlize绘制环形热图
我用的是circlize包里的circos.heatmap()函数。它的核心思路是把表达矩阵按基因（行）和细胞类型/样本（列）排好，然后映射到圆周上。
给一个可直接跑通的示例。首先准备一个表达矩阵，行是marker基因，列是细胞类型，每个值是平均表达量（z-score归一化后更直观）：
library(Seurat) library(circlize) library(ComplexHeatmap) # 目标marker列表，来自差异基因筛选和文献 markers <- c("CD3D", "CD3E", "CD8A", "GZMB", "NKG7", "MS4A1", "CD79A", "LYZ", "FCGR3A", "CST3") # 按手动注释后的celltype分组求平均表达量 avg <- AverageExpression(seu, features = markers, group.by = "celltype")$RNA avg <- as.matrix(avg) # 转置：行变成基因，列变成细胞类型 mat <- t(avg) # 按行（基因）做z-score归一化 mat <- t(scale(t(mat))) # 绘图 circos.clear() circos.par(gap.after = c(rep(5, ncol(mat) - 1), 15)) circos.heatmap( mat, col = colorRamp2(c(-1.5, 0, 1.5), c("#2166AC", "white", "#B2182B")), split = factor(colnames(mat), levels = colnames(mat)), rowname.side = "outside", cluster = TRUE ) # 添加图例 lgd <- Legend(title = "Z-score", col = colorRamp2(c(-1.5, 0, 1.5), c("#2166AC", "white", "#B2182B"))) draw(lgd, x = unit(1, "cm"), y = unit(1, "cm"), just = c("left", "bottom"))这段代码有几点要注意：
gap.after控制不同细胞类型在圆周上的间隔，最后一段留15°空间，是为了给legend和细胞类型标签留位置。cluster = TRUE会对基因做聚类，如果希望保持marker的分组顺序，可以设成FALSE或者自定义顺序。- 颜色映射用z-score的对称区间更合理，这样高表达和低表达在视觉上是对称的。
- 如果细胞类型之间表达谱差异不大，环形热图的“环”会显得非常均匀，看起来没有区分度，这时说明你选的marker或分组不合适，需要回到差异基因列表重新挑选。
5.3 环形热图之外的备选方案
环形热图不是万能的。如果你的目的是跟审稿人展示“cluster 3高表达XX、低表达YY”这种简单结论，普通的分面热图或者气泡图反而更直接。DoHeatmap()的普通热图配上group.by参数，可以把细胞类型和样本信息同时标在顶部，阅读顺序更符合直觉。
此外，如果你的差异基因列表里有大量基因属于某个通路（比如TNF signaling、interferon response），这时候更适合做通路富集分析（GSEA、clusterProfiler），而不是画一张密密麻麻的热图。环形热图是一种“锦上添花”的展示方式，不要指望它替代统计检验和功能富集。
6. 单细胞差异基因分析里，我反复见到和踩过的坑
6.1 伪复制问题：细胞不是独立样本
这是单细胞差异分析里最大的统计陷阱。假设你有3个对照样本、3个处理样本，一共1万个细胞。如果直接在细胞水平上跑差异检验，等于把同一个样本里高度相似的细胞当成独立重复，p值会非常低，几乎任何基因都会显著。这个问题的专业说法叫“伪复制”（pseudoreplication）。
正确的做法是把表达量聚合到样本水平再做差异分析。Seurat里可以用AggregateExpression()或者PseudoBulkExpression()，然后交给DESeq2或edgeR。这样样本量就变成了真实的n=3 vs n=3，统计检验的保守程度大幅提高，但结果的可信度也大幅提升。做细胞类型内部的组间比较时尤其要这样，直接拿细胞跑出来的显著基因列表，很难重复出来。
6.2 注释结果不能只看p值
我在手动注释阶段吃过一次亏：某个cluster的top差异基因里有一个非常漂亮的marker，叫XX（假设是一个成纤维细胞标志物），p值小到离谱，FeaturePlot表达也特异地集中在那一团。我当时直接把这个cluster注释成了成纤维细胞。后来把同一组的其他marker逐个看了一遍，发现这个cluster其实同时高表达内皮细胞基因和成纤维细胞基因，再回头查原始数据，发现是双细胞没去干净。
这个教训是：单个marker的p值再显著，也不能替代多个marker的交叉验证。手动注释时至少要看3-5个marker，并且要关注pct.1和pct.2的表达比例，而不是只看log2FC。如果一个cluster里几乎所有细胞都同时表达两套不同谱系的marker，先不要急着注释，回到双细胞过滤环节排查更实际。
6.3 复现性和可视化导出的经验
差异基因分析的结果要可复现，关键是设置好随机种子，并且在代码里固定Seurat版本。Seurat 4和Seurat 5在FindMarkers()的结果列名上就有差异（比如avg_logFCvsavg_log2FC），如果你换了版本，旧代码跑出来的结果很可能对不上。
我个人的习惯是：跑差异分析前先把sessionInfo()存下来，把Seurat、R的版本号以及关键参数写在一个config文件里。这样即使三个月后回来复现，也能知道当时用的什么参数。
最后，导出可视化的图片时，建议直接保存PDF或SVG矢量格式，不要只存PNG。无论环形热图还是普通热图，审稿阶段对清晰度的要求非常高，矢量图在任何缩放比例下都不会糊。用pdf()包裹绘图代码，或者用ggsave(..., device = "pdf")，这个习惯早养成早省事。
