1. 为什么你需要compareCluster?一个真实的故事
做生物信息分析的朋友,尤其是经常和转录组数据打交道的,肯定对“富集分析”这四个字不陌生。我们辛辛苦苦做完了差异表达分析,拿到了一堆上调和下调的基因,下一步自然就是想知道:这些基因到底在干什么?它们共同参与了哪些生物学过程、信号通路或者与哪些疾病相关?
最开始,我们都是老老实实地用 enrichGO 或者 enrichKEGG,一个基因列表一个基因列表地跑。比如,我有三组样本:对照组、药物处理A组、药物处理B组。常规操作就是,分别提取三组各自的差异基因(通常是上调和下调分开),然后得到6个基因列表,再分别对这6个列表做6次富集分析。这还没完,做完之后,你面前会摆着6份长长的结果表格,里面密密麻麻的通路和P值。这时候最头疼的问题来了:我怎么比较?药物A和药物B影响的通路有什么异同?上调基因和下调基因的功能倾向性又分别是什么?用眼睛一行行比对吗?那简直是灾难,不仅效率低下,而且很容易遗漏重要信息。
我当年就踩过这个坑。一个项目里比较两种抑制剂的效果,手动比对了半天,在PPT里画了三个分开的气泡图,老板看完直接问:“所以,核心的、只有A药能激活而B药不能的通路是哪条?” 我当场语塞,又得回去重新翻结果。从那次之后,我就下定决心,一定要找到一个能“并排”比较分析结果的方法。直到我发现了 clusterProfiler 包里的 compareCluster 函数,简直像打开了新世界的大门。它解决的核心痛点就是:一次性、系统性地对多个基因集进行富集分析,并将结果整合在一起,直接进行可视化比较。你不用再自己手动合并结果、自己画图去对齐,它全部帮你搞定,输出一个可以直接用于发表的可视化图形。
简单来说,compareCluster 就是为“比较”而生的。无论你是想比较不同实验条件下的基因集(比如处理vs对照),还是比较同一条件下不同方向的基因(上调vs下调),或者是你有多个完全独立的基因集合(比如来自不同疾病的标志物基因),它都能帮你一站式完成分析、整合和可视化,让你对数据的生物学解读维度,从“单个列表有什么功能”提升到“多个列表的功能差异在哪里”。这对于机制探索、药物作用模式比较、疾病亚型分型等研究场景来说,是至关重要的分析步骤。
2. 从零开始:准备你的多组基因数据
工欲善其事,必先利其器。用好 compareCluster,第一步就是把数据准备成它喜欢的格式。这一点看似简单,但很多新手容易在这里出错,导致后续分析报错。别担心,我们一步步来。
compareCluster 函数主要接受两种形式的数据输入,这也是它的两个核心应用场景。
第一种,也是最直接的一种:一个包含多个基因集的列表(List of gene vectors)。 这非常适合你已经有了若干个独立的基因集合的情况。比如,你从三篇不同的文献里收集了三个疾病相关的基因签名(Gene Signature),想看看它们富集到的通路有何异同。在R里,你需要构建一个列表,列表的每个元素都是一个字符型向量,代表一个基因集,而每个基因最好使用 Entrez Gene ID(这是最稳妥、兼容性最好的格式)。
让我们用 clusterProfiler 自带的示例数据 gcSample 来感受一下。这个数据本身就已经是完美的格式了。
# 加载必要的包和示例数据
library(clusterProfiler)
data(gcSample)
# 看看它长什么样
str(gcSample)
# 输出:List of 8 ... 说明这是一个包含8个元素的列表
# 看看每个基因集的前几个基因
lapply(gcSample, head)
运行上面的代码,你会看到 gcSample 是一个包含了8个字符向量的列表,每个向量里都是一串Entrez ID。这就是 compareCluster 需要的标准“多基因集”格式。你自己的数据也可以很容易地做成这样,比如:
# 假设你有三个差异基因列表,分别是上调、下调和所有差异基因
up_genes <- c("1234", "5678", "9012") # 这里替换成你的Entrez ID
down_genes <- c("3456", "7890")
all_de_genes <- c("1234", "5678", "9012", "3456", "7890")
# 将它们组合成一个列表
my_gene_list <- list(
Upregulated = up_genes,
Downregulated = down_genes,
All_DE = all_de_genes
)
第二种输入方式,也是功能更强大的一种:基于数据框的公式(Formula)输入。 这种方式特别适合你手头有一个包含所有基因及其分组信息的大表格。比如,你做完差异分析,得到了一个数据框,里面每一行是一个基因,有基因ID列,还有一列或多列表示这个基因属于哪个分组(例如,根据表达变化方向分组“up/down”,或者根据表达模式聚类分组“Cluster1/2/3”)。
我们再用包里的另一个示例数据 geneList 来演示如何构建这种格式。geneList 是一个带有名字的数值向量,名字是Entrez ID,值是log2FoldChange。
# 加载数据
data(geneList, package="DOSE")
head(geneList)
# 将它转换成一个数据框
mydf <- data.frame(Entrez = names
432



被折叠的 条评论
为什么被折叠?



