单细胞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主要来自样本Acluster 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 0logfc.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。这样样本量就变成了真实的n3 vs n3统计检验的保守程度大幅提高但结果的可信度也大幅提升。做细胞类型内部的组间比较时尤其要这样直接拿细胞跑出来的显著基因列表很难重复出来。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)这个习惯早养成早省事。
单细胞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主要来自样本Acluster 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 0logfc.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。这样样本量就变成了真实的n3 vs n3统计检验的保守程度大幅提高但结果的可信度也大幅提升。做细胞类型内部的组间比较时尤其要这样直接拿细胞跑出来的显著基因列表很难重复出来。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)这个习惯早养成早省事。
