教你用clusterProfiler实现其它来源富集结果的可视化

阿凡亮

发布于 2020-04-14 15:04:40

3.4K2

发布于 2020-04-14 15:04:40

引言

在对功能富集分析的结果进行可视化的时候，大家肯定都听过Y叔的R包clusterProfiler，这个包可以说是富集分析结果可视化的神器，不仅画出来的图好看，而且种类繁多，可以满足各种需求；但要想用clusterProfiler进行可视化就必须用它进行富集分析；对于其他软件的富集结果，是不能够进行可视化的。那么能不能通过clusterProfiler对其他软件的富集结果进行可视化分析呢？答案是肯定的，本期就给大家分享一下怎么通过clusterProfiler对其他软件的富集结果进行可视化。

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，也就是富集到的基因数目。

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)