利用clusterprofiler包对非模式物种进行GO和KEGG富集分析

阅读数: 次 2020-08-02

写在文章之前，首先感谢Y叔写出来的clusterprofiler包。但是从事生物信息分析工作的我们都知道，很多非模式物种并没有构建好的物种注释包
这个时候就需要我们自己构建，作为初学者的我，也仿照着Y叔的个人博客来进行的分析

获取非模式物种的注释信息

这里我们需要使用的软件是eggnog-mapper

1.首先安装eggnog-mapper软件

使用git下载最新版的emmapper，最新版的软件已经集成了他所需要各个依赖

1	git clone https://github.com/eggnogdb/eggnog-mapper.git

2.使用软件得到注释文件集合

进入软件的目录，使用软件自带的脚本进行下载，请确保你的系统环境拥有python2环境

1 2	$ cd data $ ../download_eggnog_data.py

物种注释包有点大，需要花不少时间去下载，下载完成后需要进行解压缩
当然你也可以通过文件的网址直接进行进行wget下载。也不是不行不是吗？当然要放在data文件夹下面哦

然后拥有了这些准备文件之后，我们需要一个拿来构建的物种的全部蛋白序列，这里就不提供了。你可以去ncbi或者其他地方进行下载
本次使用的蛋白质数据大于1000条序列，所以我们使用diamond进行比对；

运行之前，请查看一下你是否已经安装了指定版本的biopython，如若没有，还请安装一下

$ emapper.py -m diamond \
           -i input.pep \
           -o diamondoutput \
           --cpu 20
$  sed -i '/^# /d' diamond.emapper.annotations 
$  sed -i 's/#//' diamond.emapper.annotations

耐心等待，记得给这一步留出足够多的内存哦。待其结束，我们就获得了一个自己物种的注释集合，可以拿来用于构建orgDb了。

PS

事后才发现这个软件居然有一个网页版的http://eggnogdb.embl.de/#/app/emapper
啊，早有网页版的我装什么本地的。白嫖他不香吗？我哭了

进行GO分析

本人至今对于在无图形化的服务器上运行R程序一直感到深恶痛绝，所以基本上我的R程序都是本地运行的。
在运行clusterprofiler之前，我们需要提取注释信息

如果没有安装AnnotationHub和clusterProfiler可以使用BioManager安装，十分推荐

#加载所需要的包
$ library(AnnotationHub)
$ library(stringr)
$ library(dplyr)
$ library(clusterProfiler)
#读取数据
$ egg_f<-"diamond.emapper .annotations"
$ egg <- read.csv(egg_f, sep = "\t")
$ egg[egg==""]<-NA
#提取数据
$ gterms <- egg %>%
  dplyr::select(query_name, GOs) %>% na.omit()
$ gene2ko <- egg %>%
  dplyr::select(GID = query_name, Ko = KEGG_ko) %>%na.omit()
$ gene_info <- egg %>%
  dplyr::select(GID = query_name, GENENAME = `eggNOG.free.text.desc.`) %>% na.omit()
$ gene2go <- data.frame(term = character(),
                      gene = character())
$ for (row in 1:nrow(gterms)) {
  gene_terms <- str_split(gterms[row,"GOs"], ",", simplify = FALSE)[[1]]
  gene_id <- gterms[row, "query_name"][[1]]
  tmp <- data_frame(gene = rep(gene_id, length(gene_terms)),
                    term = gene_terms)
  gene2go <- rbind(gene2go, tmp)
}
#输出数据，如果你需要的话
$ write.table(file="gene2go.table",gene2go，sep="\t")
#构建clusterprofiler的GO富集数据并进行GO富集
$ gene_list<-gene2go$gene[1:40]
$ term2gene<-gene2go[,c(2,1)]
$ df<-enricher(gene=gene_list,
             pvalueCutoff = 0.05,
             pAdjustMethod = "BH",
             TERM2GENE = term2gene)
#转换坐标
$ df<-as.data.frame(df)
$ df1<-go2term(df$ID)
$ write.table(file="go2ont.table",df2,sep="\t")
$ df2<-go2ont(df$ID)
$ df2<-go2ont(df$ID)
$ df$Ont<-df2$Ontology
$ df3<-df%>%
  select(c("term","Ont","pvalue")) #根据需要选择自己想要绘制的列进行绘制
$ library(ggplot2)
$ ggplot(df3,aes(x=term,y=-log10(pvalue)))+
  geom_col(aes(fill=Ont))+
  coord_flip()+labs(x="")+
  theme_bw()

~~KEGG的富集分析类似于GO富集，只是你需要kegg_info.RData这个文件来提供pathway2name这个信息~~
~~剩下得内容基本和GO富集一致，我在这里就不继续写下去了~~
参考资料：
https://guangchuangyu.github.io/cn/
https://www.jianshu.com/p/9c9e97167377
https://www.jianshu.com/p/9c9e97167377