做单细胞测序分析的同学应该都有这种体会:经过聚类、UMAP降维、差异表达分析之后,我们手里攒了一大堆marker基因,每个cluster动辄几十上百个基因,看表头看得眼花缭乱,但真要写进论文、解释生物学意义的时候,却发现这些基因列表根本不"说话"——它们不会告诉你这个cluster是T细胞还是巨噬细胞,也不会告诉你这些细胞正在执行什么样的生物学功能。于是就有了单细胞测序流程里的关键一步:marker基因转化和GO富集分析。简单说,就是把每个cluster的marker基因从"基因符号列表"变成"功能注释结果",通过统计检验看看这些基因在哪些生物学通路、分子功能或细胞组分上显著富集,这一步是串联"差异基因"和"细胞类型/生理机制"的核心桥梁,也是很多审稿人必看的分析环节。这篇文章我会把ID转换、富集原理、实操代码和排查经验一次讲透,适合已经跑完单细胞上游分析、正在做下游功能解读的科研人员和生信初学者。
1. 单细胞测序marker基因:从细胞聚类到功能洞察的关键一步
1.1 为什么marker基因需要"二次加工"
先想一个问题:我们从Seurat或Scanpy里拿到的那张marker基因表,本质上是什么?它是每个cluster相较于其他cluster差异表达显著的基因集合。比如cluster 0高表达CD3D、CD3E、IL7R,那我们可以初步推断这是T细胞;cluster 2高表达MS4A1、CD79A,那可能是B细胞。这种判断靠的是"人肉背诵"经典的细胞标志基因,经验丰富的老师一眼就能看出来。但面对大量非典型marker基因,或者遇到罕见细胞群、疾病状态下的异常细胞时,这种"人肉识别"就完全行不通了。
这时候就轮到GO富集分析上场了。它的逻辑很朴素:一个cluster里的marker基因如果在"免疫应答调节"这个GO条目上显著富集,那即便我们没有认出一个已知的T细胞标志物,也能从功能层面推断这个细胞群高度参与免疫调控。这种从"基因列表"到"功能模块"的提升,就是marker基因二次加工的核心价值。顺带说一句,GO富集分析不仅能辅助细胞类型注释,在后续的拟时序分析、细胞通讯解读中也能提供重要的背景支持。
另外,从工作流程的角度看,marker基因转化这一步还承担着"数据规范化"的作用。不同平台、不同版本注释文件下得到的基因symbol格式可能不一致,有的带版本号(比如ENSG00000167286.17),有的不带;而富集分析工具(尤其是R生态的clusterProfiler)往往对输入格式有严格要求。如果直接把乱七八糟的基因格式塞进去,报错不说,还可能因为基因注释不完整导致假阴性结果。所以,这里的"转化"包含了两个层次:一是标识符格式的统一,二是基因ID与数据库注释之间的映射。
1.2 这一环节在整个单细胞分析流程中的定位
单细胞测序流程走到marker基因转化和GO富集,已经是整个链条比较靠后的位置了。完整走一遍大致是这样的:原始测序数据质控(fastp、Cell Ranger等工具)→ 表达定量与细胞过滤 → 降维聚类 → 差异表达分析鉴定marker基因 → 细胞类型注释 → 功能富集分析(GO/KEGG等)→ 下游高级分析(拟时序、细胞通讯、SCENIC等)。也就是说,我们在这篇文章中做的事,是用"聚类注释后的静态结果"去回答"这些细胞群到底在干什么"的问题。
这一步做得好不好,直接影响后面几件事。第一,如果GO富集结果准确、有生物学意义,你写论文的结果部分会非常有底气,Discussion部分也有据可依。第二,富集结果反过来可以验证细胞注释是否合理——如果某个cluster注释为CD8+ T细胞,但GO富集结果显示"脂肪酸代谢"相关通路异常高富集而免疫相关的GO类别很少,那这个注释十有八九有问题。第三,富集结果通常也是后续做拟时序、转录因子分析时挑选感兴趣基因模块的重要依据。可以说,这一步是单细胞转录组分析"从数据到生物学结论"的必经之路。
2. 基因ID转化:别让小细节卡住大流程
2.1 常见的基因ID格式与转化场景
很多人第一次跑GO富集时,最崩溃的不是算法跑不动,而是第一步就报错——"Error in check_gene_id()"。原因很简单:clusterProfiler的enrichGO函数默认要求输入的基因是ENTREZ ID,而不是我们在Seurat的marker表里最常见的基因symbol(比如TP53、CD3D这种)。
那常见的基因ID有哪些呢?我梳理一下,方便新手建立整体概念:
| ID类型 | 例子 | 说明 |
|---|---|---|
| SYMBOL | TP53, CD3D | 人类易读的基因简称,日常最常用 |
| ENTREZ ID | 7157, 916 | NCBI维护的数字ID,很多R包的首选输入格式 |
| ENSEMBL ID | ENSG00000141510 | Ensembl数据库的基因ID,可能带版本后缀 |
| GENCODE ID | 一般也以ENSG开头 | 通常与ENSEMBL ID相对应,但版本管理方式不同 |
| RefSeq ID | NM_000546 | 基于转录本序列的ID,常用于核酸序列层面 |
| UniProt ID | P04637 | 蛋白质层面的ID,蛋白质组学常用 |
为什么偏偏要ENTREZ ID?因为clusterProfiler的内部逻辑是:将输入的基因ID映射到OrgDb(如org.Hs.eg.db)中预先构建的GO注释关系,而OrgDb的中央枢纽就是ENTREZ ID。也就是说,无论你传SYMBOL还是ENSEMBL ID,它都会先转换成ENTREZ ID再去检索GO注释。与其让它转,不如我们自己在前面控制好这一步,至少我们能清楚看到有多少基因转换成功、有多少丢失了。
还有一个常见场景:有的同学拿到的是Ensembl ID转录本版本号(比如ENST开头),还有一些单细胞注释工具输出的是基因全名。这些格式在GO富集分析中都不是直接可用的,预处理阶段必须统一。
2.2 两种最常用的转化方案:bitr()函数和biomaRt
方案一:clusterProfiler自带的bitr()函数。这是我最推荐新手使用的方式,因为代码少、依赖少、报错信息友善。基本用法如下:
library(clusterProfiler) library(org.Hs.eg.db) # gene_list是你从Seurat的marker表里提取的基因symbol向量 gene_list <- c("TP53", "CD3D", "MS4A1", "CD79A", "IL7R") gene_entrez <- bitr(gene_list, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)运行完之后,gene_entrez是一个两列的数据框,左边是原始的SYMBOL,右边是对应的ENTREZID。如果某些基因在OrgDb中找不到注释,bitr会自动丢弃它们并在返回结果里体现为"没出现"。要注意,bitr不会报错告诉你丢了哪些基因,你需要自己对比一下输入和输出的基因数量,做好记录。这一步看似简单,但直接关系到后续富集的基因数量和质量。
方案二:biomaRt包。这个是从Ensembl数据库在线取数的方式,灵活性更高,能自定义选择物种、数据库版本,还能批量获取多种注释信息(比如基因在染色体上的位置、序列信息等)。基础用法:
library(biomaRt) # 选择人类基因数据集,注意这里的host可以根据Ensembl版本切换 ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl") gene_map <- getBM(attributes = c("hgnc_symbol", "entrezgene_id"), filters = "hgnc_symbol", values = gene_list, mart = ensembl)说句实话,如果只是单细胞marker基因的GO分析,用biomaRt反而有点"杀鸡用牛刀"了。它更适合那种需要大量、多物种、跨数据库ID映射的复杂场景,比如做进化比较或需要自定义基因注释信息的时候。对于大多数场景,bitr + OrgDb已经足够,而且不依赖网络,跑起来快得多。
2.3 转化时容易踩的三个坑
第一个坑是物种匹配错误。人源数据一定要用org.Hs.eg.db,小鼠用org.Mm.eg.db,大鼠用org.Rn.eg.db,斑马鱼用org.Dr.eg.db。如果你人的marker基因拿小鼠的OrgDb去转换,结果不仅是转换率极低,更可怕的是某些基因会映射到同源基因上,导致后续富集结果完全错误。这种错误挺隐蔽的,因为它不一定报错,只是结果很怪。
第二个坑是gene symbol的大小写。不同数据库对大小写的敏感度不同,尤其是人类基因symbol规范是大写,但部分小鼠基因是首字母大写、其余小写(比如Cd3d)。如果你的输入格式和参考数据库不一致,转换率会断崖式下降。处理办法是先用toupper()或Hmisc::capitalize()规范一下格式。但要注意,不是所有基因都是纯大写或首字母大写,极少数基因像"Sept7"这类拼写比较特殊,处理时要留个心眼。
第三个坑是基因名称的版本后缀。如果你从Ensembl下载的注释文件里带版本号(如ENSG00000141510.17),直接拿去转ENTREZ ID往往匹配失败。需要用gsub("\..*", "", x)把后缀去掉。这个我在实际项目中踩过一次,当时差点把几万个基因的版本后缀当成了正常格式,前几次转换率只有30%,排查了半天才发现是小尾巴闹的。
3. GO富集分析的原理:为什么是GO而不是KEGG
3.1 GO的三大类(BP/MF/CC)及适用场景
基因本体论(Gene Ontology,GO)是把基因功能用一套结构化的、有向无环图(DAG)形式的词汇表来描述的体系。它分成三大类,对应生物学功能的不同层面:
- BP(Biological Process,生物学过程):描述基因参与的生命过程,比如"免疫应答"、"炎症反应调整"、"细胞周期"、"DNA修复"。简单理解就是"细胞正在做什么事"。
- MF(Molecular Function,分子功能):描述基因在分子层面的活性,比如"激酶活性"、"转录因子结合"、"钙离子结合"。这是"基因作为分子有什么功能"。
- CC(Cellular Component,细胞组分):描述基因产物在细胞中的位置,比如"细胞膜"、"线粒体"、"核小体"。这是"在哪儿干活"。
在单细胞marker基因注释场景中,我个人用得最多的是BP。因为marker基因反映的是细胞类型特异性的功能状态,而BP最能直接对应生物学过程和表型。MF用得相对少一些,除非你做的是特定分子功能的精细对比。CC常用于补充验证,比如你发现某个cluster的marker基因高度富集在"核糖体"和"线粒体内膜",那可能提示这个cluster是代谢活跃的细胞,或者在质控时线粒体基因污染偏高。
3.2 超几何检验:富集分析的统计学底座
现在聊一聊GO富集分析的统计原理。很多同学把enrichGO函数一跑,出来一堆p值就完事了,但理解背后的统计逻辑对解读结果非常重要。
GO富集分析的经典方法是基于超几何分布(hypergeometric distribution)的富集检验,也可以用Fisher精确检验来等价实现。它的问法很直接:给定一个基因总数(背景,通常是基因组所有编码基因),其中有M个基因被注释到了某个GO条目上;我们的marker基因列表有K个基因,其中X个正好落在该GO条目里。那么,如果marker基因是随机抽出来的,抽出X个甚至更多基因落在该GO条目的概率有多大?这个概率就是p值。
# 以超几何检验为例,直观感受一下 # 想象总基因N=20000,某GO条目相关基因M=100,marker基因K=200,X=15 # p = sum_{i=X}^{min(K,M)} (choose(M,i) * choose(N-M,K-i)) / choose(N,K)p值小,说明富集不是随机发生的,越小的p值意味着越显著的富集。但要注意,单次检验的p值只对照一个GO条目,而GO有成千上万个条目,每次都算一遍p值,必然会遇到大量假阳性。所以实际分析中必须做多重假设检验校正。clusterProfiler默认用BH(Benjamini-Hochberg)方法校正,得到adjusted p值,也就是q值,然后用q值阈值(默认0.2)筛选显著富集的GO条目。这里有两个阈值参数:pvalueCutoff和qvalueCutoff。pvalueCutoff通常设0.05,qvalueCutoff可以稍微放宽到0.1甚至0.2,具体看你的富集结果数量和研究的探索性程度。
3.3 为什么建议先做GO再做KEGG
这个环节也有不少新手会直接问:能不能只做KEGG不做GO?我的回答是,如果你想要更全面的生物学解释,GO是不能跳过的,原因有两点。第一,GO的覆盖范围远大于KEGG。KEGG分析的是代谢通路和信号通路,很多基因没有对应的KEGG条目,尤其是非编码调控和未充分研究的基因。GO在定义上覆盖所有有注释的基因功能,覆盖率更高。第二,GO提供了BP、MF、CC三个层面的解释,而KEGG偏重于"通路"这一种线性视角。单细胞转录组数据分析中我们常常需要从多个维度解释一个细胞群的功能,GO的这一结构优势很难替代。
当然,KEGG也没有必要完全放弃。KEGG通路图直观、有严格的通路层级关系,在论文里展示的时候可读性非常好。常见的做法是:先用GO富集从功能层面"圈出"生物学过程,再用KEGG验证这些过程有没有具体的通路支撑。这样两种富集结果互相印证,结论更稳健。
4. 完整实操:从marker基因列表到GO富集结果
4.1 准备工作与数据输入
在我们开始跑代码之前,先把环境和数据准备好。你需要有:R(版本4.0以上比较稳妥)、必要的R包(Seurat、clusterProfiler、org.Hs.eg.db、ggplot2、dplyr等),以及一份已经跑完FindAllMarkers的Seurat对象。如果你手里暂时没有自己的数据,可以用Seurat官方教程里的pbmc_small或者pbmc3k来练习,效果是一样的。
# 安装所需包(如果还没装的话) # BiocManager::install(c("clusterProfiler", "org.Hs.eg.db")) library(Seurat) library(clusterProfiler) library(org.Hs.eg.db) library(dplyr) library(ggplot2)关于数据的输入格式,需要提醒的是:marker基因列表最好来自差异表达分析的结果表格,而不是你自己随便选的几个基因。因为GO富集分析是统计检验,它需要足够的基因数量才能产生可靠的结果。如果只拿十来个marker基因去做富集,统计功效很低,结果基本不可用。一般建议至少每个cluster选50~200个显著差异基因做GO富集。
4.2 差异marker提取与基因ID转化实操
假设你的Seurat对象叫pbmc,已经跑完了NormalizeData、FindVariableFeatures、ScaleData、RunPCA、RunUMAP、FindNeighbors、FindClusters这些步骤,接下来就是FindAllMarkers:
# 找所有cluster的marker基因 all_markers <- FindAllMarkers(pbmc, only.pos = TRUE, # 只保留上调基因 min.pct = 0.25, # 在至少25%的细胞中表达 logfc.threshold = 0.25) # log2FC阈值返回的all_markers是一个数据框,每一行是一个基因在一个cluster中的差异表达检验结果。列包括p_val、avg_log2FC、pct.1、pct.2、p_val_adj、cluster、gene这些字段。我们会按cluster分组,每个cluster取top50的基因:
top50_per_cluster <- all_markers %>% group_by(cluster) %>% top_n(n = 50, wt = avg_log2FC) # 取出所有去重后的基因symbol all_genes <- unique(top50_per_cluster$gene) # ID转化:SYMBOL -> ENTREZID gene_entrez <- bitr(all_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) cat("原始基因数:", length(all_genes), "\n") cat("成功转换数:", nrow(gene_entrez), "\n") cat("转换率:", round(nrow(gene_entrez) / length(all_genes) * 100, 2), "%\n")跑完之后会看到一个转换率。正常情况下人源数据转换率在75%~95%之间,如果转换率低于60%就要怀疑前面那几个"坑"了。这里有一个细节:bitr的结果里没有转换成功的基因不会出现在输出中,所以你想知道丢了哪些基因,可以用setdiff(all_genes, gene_entrez$SYMBOL)导出看看。
4.3 enrichGO函数参数详解与运行
ID转换完成后,下一步就是富集分析的核心函数enrichGO()。这里参数比较多,我逐个讲一下常用参数的含义和一个合理的推荐配置:
# 以单个cluster为例演示,先取cluster 0的数据 cluster0_symbols <- top50_per_cluster %>% filter(cluster == 0) %>% pull(gene) cluster0_entrez <- bitr(cluster0_symbols, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) ego_cluster0 <- enrichGO(gene = cluster0_entrez$ENTREZID, # 输入基因的ENTREZID向量 OrgDb = org.Hs.eg.db, # 物种注释包 keyType = "ENTREZID", # 输入基因ID类型 ont = "BP", # 本体类别:BP/MF/CC/ALL pAdjustMethod = "BH", # 多重检验校正方法 pvalueCutoff = 0.05, # p值阈值 qvalueCutoff = 0.2, # q值阈值 readable = TRUE) # 结果中显示原始SYMBOL参数含义拆解:
- gene:必填,基因ID向量。注意这里的ID类型要和keyType参数保持一致,千万别传SYMBOL却设置keyType = "ENTREZID"。
- OrgDb:物种相关的OrgDb对象,决定了基因和GO注释之间的映射来源。
- ont:可选"BP"、"MF"、"CC"或"ALL"。实际分析中建议分别跑,最后再汇总展示。
- pAdjustMethod:多重假设检验校正方法,BH是默认也是比较经典的,但如果你样本量特别大或者GO条目特别多,也可以考虑"bonferroni",结果会更严格。
- pvalueCutoff和qvalueCutoff:这两个阈值决定了结果中保留哪些GO条目。通常pvalueCutoff设0.05,qvalueCutoff设0.2。如果你的结果太少,可以适当放宽;结果太多,可以收紧。
- readable:设为TRUE后,结果表格里会在ENTREZID旁边附加基因SYMBOL,方便阅读。这个默认是FALSE,我强烈建议打开。
跑完enrichGO之后,得到的是一个enrichResult对象。你可以用head()直接查看结果表格,也可以用as.data.frame()转成常规数据框。每一行代表一个显著富集的GO条目,关键列有ID、Description、GeneRatio、BgRatio、pvalue、p.adjust、qvalue、geneID和Count。其中GeneRatio是"命中的基因数 / 输入的基因数",BgRatio是"该GO条目的背景基因数 / 总背景基因数"。简单说,GeneRatio越高,说明这个GO条目在你关注的基因集中覆盖度越高。
4.4 结果表格解读:从p值到富集分数
看到结果表格后,怎么判断哪些条目值得放进论文?我的经验是先看q值排序(从小到大),然后关注GeneRatio和Count。q值最小说明最显著,GeneRatio高说明这个功能在你的基因集里的覆盖度高。一般我会选Top 5~10的GO条目来解读。
举个例子,假设cluster 0的结果里有"GO:0002376 immune system process"和"GO:0006955 immune response"排在最前面,q值都远小于0.05,GeneRatio分别是0.2和0.15,这就比较清晰:这个cluster高度参与免疫应答。配合经典T细胞marker(CD3D、CD3E)的高表达,基本可以断定这个cluster是T细胞。但如果你发现一个cluster的富集结果全是"核糖体生物合成"和"蛋白质靶向膜"这类典型的"管家功能",那就要回头看看是不是这个cluster质量有问题,或者是不是降维聚类参数不合适,导致没有真正分离出有特殊功能的细胞群。
4.5 可视化:barplot与dotplot的能力边界
富集分析结果只用表格展示既枯燥又不直观,论文里通常要配图。clusterProfiler自带几种可视化方式,最常用的是barplot和dotplot。
# 柱状图 barplot(ego_cluster0, showCategory = 15) # 气泡图 dotplot(ego_cluster0, showCategory = 15)两个图各有特点。barplot展示的是富集得分(一般是-log10(p.adjust))从高到低排列的前若干条目,适合快速呈现富集显著性;dotplot在垂直轴上放GO条目、水平轴放GeneRatio,气泡大小对应Count、颜色对应p.adjust,信息含量更丰富,我一般优先选择dotplot。
不过要提醒的是,这两种图都有"只显示Top N"的限制,默认showCategory = 10或15。如果你的结果特别多,最好根据q值排序后手动挑选更有生物学意义的条目来呈现,不要盲目展示数量过多。还有一个细节:默认的dotplot颜色渐变用的是p.adjust,在黑白打印时不容易区分,可以手动调整color参数,或者后期用AI/PPT加标注。此外,如果你对不同cluster的富集结果做横向比较,建议用compareCluster()函数先合并再可视化,这样的对比图在单细胞论文中非常常见,也相当出彩。
5. 实操中常见的五类问题与排查实录
5.1 ID转换失败率过高怎么办
转换率低是大家问得最多的问题。我总结下来原因通常是这几类:物种用错、gene symbol大小写不规范、基因本身是老的synonym(别名)而非官方symbol、带Ensembl版本后缀。排查步骤可以按顺序来:先核对物种,再统一大小写,然后用bitr的多区块转换尝试。bitr还支持从一个类型同时映射到多种ID,例如fromType = "SYMBOL", toType = c("ENTREZID", "ENSEMBL"),这可以帮助一次性检查多个ID体系的映射情况。如果某基因在组织里找不到,还可以用alias2SymbolTable()(来自limma包)把别名转换成官方symbol再转一轮。
5.2 富集结果为空或显著项太少
这种情况通常有两种可能:一是输入的基因数量太少,统计功效不足;二是qvalueCutoff设得太严格。先看输入基因数量,如果只有30个基因,那没有显著富集是正常的。一种缓解的思路是把marker基因数量增加,比如从top50扩到top200,但这样也会引入噪声,需要平衡。另一种思路是适当放宽qvalueCutoff到0.1或0.2,同时结合pvalueCutoff=0.05来综合判断。需要注意的是,如果放宽后仍然无显著富集,那就要检查背景基因(universe)设置是否合理了——有些人会把整个基因组作为背景,有些人使用所有检测到的基因,这会对统计结果造成比较大的影响。在clusterProfiler中,默认背景是OrgDb里的全库注释基因,一般问题不大,但如果你想自定义背景,enrichGO支持universe参数,建议以自己数据中检测到的所有基因作为背景。
5.3 结果多到看不过来:如何筛选与精简
另一个"幸福的烦恼"是富集结果太多,动辄几十上百个显著GO条目。这时不要贪多,建议按以下思路筛选:先锁定最能代表生物学问题的三个方向,比如免疫应答、细胞增殖、代谢重编程,然后再从这些方向里找一个具体的显著GO条目深挖。展示时常用的门控方式是只保留q值最小的前10~20个条目,同时优先选择那些覆盖度较高的BP条目。一旦发现某一类结果特别多,比如"蛋白转运"相关条目占到一半,那大概率是marker基因列表质量偏低,混入了大量管家基因,建议重新调整差异表达阈值或者换一种marker筛选策略。
5.4 不同版本数据库之间的差异
这是一个不可忽视的问题。同一个cluster的marker基因,用org.Hs.eg.db 3.14版本和3.19版本去做GO富集,得到的显著GO条目很可能不一样——因为数据库会不断地添加新的基因和新的注释。不少同学在复现别人论文结果时发现对不上,很多时候就是这个原因。建议在发表论文时,在Methods里明确写出使用的包版本号,最好还加上sessionInfo()的信息。另外一个实用建议是尽量使用统一的版本环境来分析同一个项目的所有样本,不要把不同时间的分析结果直接放一起比较。
5.5 背景基因(universe)设置的正确姿势
最后说说背景基因。在富集分析里,背景基因就是"所有可能出现的基因集合"。如果背景里混入了差异性很大的物种基因,或者把一些根本不在你实验体系里表达的基因也塞进去,会稀释统计功效。最好的做法是:使用你单细胞数据中实际检测到的所有基因作为背景。不要用全基因组注释。我遇到过一些同学不设置universe参数,结果表面看富集得很漂亮,但实际上很多富集到的GO类别是"假阳性",根源就是背景选得不对。clusterProfiler中设置背景的做法是:
# 假设你检测到的所有基因是 all_detected_genes universe <- bitr(all_detected_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID ego_cluster0 <- enrichGO(gene = cluster0_entrez$ENTREZID, universe = universe, OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2, readable = TRUE)6. 一些经验之谈:怎样让富集结果真正服务于你的课题
这个环节走到最后,我特别想分享的一点是:GO富集分析不是终点,而是帮你把"基因列表"变成"生物学故事"的枢纽。在实际工作中,我不会一上来就机械地对每个cluster跑一遍enrichGO,而是先快速浏览marker表,挑出几个明显有“身份特征”的基因做初判,再针对性做富集。比如先看CD3D、CD79A、LYZ这些泛系别标志物,猜出大体细胞类型,然后用GO富集去验证和补充——这样拿到的手结果可信度更高,解读时也更有的放矢。
另外,多cluster对比我曾经踩过不少坑。随着cluster数量增加,逐一跑enrichGO再手动拼图非常麻烦,代码也容易出错。建议从一开始就用compareCluster()把所有cluster的marker基因一次性放进去做统一分析,输出一个可以画成分组dotplot的整体结果,效率高、自己看图也直观。脚本是这样:
# 构建一个以cluster分组的长格式数据框 marker_list <- split(top50_per_cluster$gene, top50_per_cluster$cluster) # 批量GO富集并比较 comp_ego <- compareCluster(geneCluster = marker_list, fun = "enrichGO", OrgDb = org.Hs.eg.db, keyType = "SYMBOL", ont = "BP", pvalueCutoff = 0.05, qvalueCutoff = 0.2, readable = TRUE) dotplot(comp_ego, showCategory = 10)compareCluster会用SYMBOL作为输入,内部自动处理ID映射,输出结果按cluster分别富集,方便横着对比。我在多个单细胞项目里一直沿用这套做法,配合每个cluster单独跑一遍enrichGO做精细检查,两相对照,既能提高效率又能保证关键结论的稳健。最终你在论文里放的那些GO富集图,往往不用做得花里胡哨,只要筛选出的功能条目清楚、逻辑自洽,审稿人基本就不会在这个环节上挑毛病。