在对功能富集分析的结果进行可视化的时候,大家肯定都听过Y叔的R包clusterProfiler,这个包可以说是富集分析结果可视化的神器,不仅画出来的图好看,而且种类繁多,可以满足各种需求;但要想用clusterProfiler进行可视化就必须用它进行富集分析;对于其他软件的富集结果,是不能够进行可视化的。那么能不能通过clusterProfiler对其他软件的富集结果进行可视化分析呢?答案是肯定的,本期就给大家分享一下怎么通过clusterProfiler对其他软件的富集结果进行可视化。
01
clusterProfiler的可视化对象格式解析
首先通过clusterProfiler的样例数据进行富集分析,代码如下:
library(clusterProfiler)
#导入示例数据data(geneList)
#得到差异表达的基因
de <- names(geneList)[abs(geneList) > 2]
#进行富集分析
edo <- enrichDGN(de)
#生成柱状图的富集结果,并只显示前10类
barplot(edo, showCategory=10)
柱状图富集结果展示如下:
通过代码可以发现富集分析的可视化是通过变量edo完成的,因此只要分析清楚变量edo的结构,就可以对其他软件的富集结果进行可视化分析。首先通过R里面的str函数对变量edo的格式进行查看,变量的edo的格式如下所示:
从图中可以看到edo的格式是一个名字叫enrichResult的R的class类型,因此我们只要将其他软件的富集结果改写成这个class,就可以用clusterProfiler进行可视化了,通过查找DOSE(clusterProfiler的一个依赖包)包,得到了enrichResult这个class的格式,具体格式如下所示:
setClass("enrichResult",
representation=representation(
result = "data.frame",
pvalueCutoff = "numeric",
pAdjustMethod = "character",
qvalueCutoff = "numeric",
organism = "character",
style="margin: 0px; padding: 0px; color: rgb(221, 17, 68);">"character",
gene = "character",
keytype = "character",
universe = "character",
gene2Symbol = "character",
geneSets = "list",
readable = "logical"
),
prototype=prototype(readable = FALSE)
)
由格式可以看出这个类中最重要的部分是result这个变量,因为它是一个dataframe,所以是用来存储富集结果的。因此继续对result的格式进行查看,方法还是用str函数,得到result的变量格式如下:
从图中可以看到result总共有9列,第一列是ID,也就是富集通路的编号;第二列是Description,也就是富集通路的名称;第三列是GeneRatio,也就是要富集的基因中在对应通路中的比例;第4列是BgRation,也就是对应通过的基因在全基因组注释中的比例;第5,6,7列都是统计检验的结果,可以按对应类型进行填写,如果只有一个可以填写一样的;第8列是geneID,也就是富集到基因的名字,它的格式是以斜线隔开的;第9列是Count,也就是富集到的基因数目。
02
agrigo结果实例演示
接下来以argigo的富集结果为例子进行演示,首先从agrigo的网站上下载得到agrigo的富集结果,具体结果如下所示:
由agrigo的结果可以看出,它包含clusterProfiler可视化结果数据的所有信息,所以可以用clusterProfiler进行可视化。
然后将agrigo的结果改写成enrichClass的格式,具体代码如下所示:
#读取agrigo的富集结果
agrigo <- read.table('375624573.1.txt',header = T,sep='\t',quote="")
#修改agrigo的geneID为clusterprofiler的geneID格式
geneID<-gsub(' // ','/',agrigo$entries)
geneID<-gsub('// ','',geneID)
#得到富集到的所用geneID
geneID_all <- unlist(apply(as.matrix(geneID),1,function(x) unlist(strsplit(x,'/'))))
#生成enrichResult类的result属性格式
Over <- data.frame(ID = as.character(agrigo$GO_acc),
Description = agrigo$Term,
GeneRatio = apply(agrigo[,c("queryitem","querytotal")],1,paste,collapse='/'),
BgRatio = apply(agrigo[,c("bgitem","bgtotal")],1,paste,collapse='/'),
pvalue = agrigo$pvalue,
p.adjust = agrigo$FDR,
qvalue = agrigo$FDR,
geneID = geneID,
Count = agrigo$queryitem,
stringsAsFactors = FALSE)
rownames(Over)<-Over$ID
#生成agrigo结果的enrichResult类
x <- new("enrichResult",
result = Over,
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
qvalueCutoff = 0.05,
gene = geneID_all,
universe = 'Unknown',
geneSets = list(),
organism = "Unknown",
keytype = "Unknown",
style="color: rgb(221, 17, 68); margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px; padding-bottom: 0px; padding-left: 0px; padding-right: 0px; padding-top: 0px;">"agrigo",
readable = FALSE
)
注:代码中的pvalueCutoff、pAdjustMethod、qvalueCutoff的数据在这里只是标识,没有实际含义,各自可以按照自己意愿进行改写。
最后对结果以柱状图的形式进行可视化:
barplot结果如下:
barplot(x, showCategory = 10)
dotplot结果如下:
dotplot(x,showCategory = 10)
upsetplot代码如下:
upsetplot(x)
总结
这里我们基本上完成了clusterProfiler可视化对象的构建,还有一些细节没有考虑,比如富集结果的显著性的的过滤;感兴趣可以自己对代码再进行修改;clusterProfiler对富集结果的可视化方式还有很多,感兴趣的同学也可以试试其他的可视化方式,这里就不在展示了。