第一步是使用eggnog-mapper做功能注釋

conda activate emapper
python emapper.py -i orgdb_example/GCF_000002945.1_ASM294v2_protein.faa --output orgdb_example/out -m diamond --cpu 8

將注釋結(jié)果下載到本地，手動(dòng)刪除前三行帶井號(hào)的行，第四行開(kāi)頭的井號(hào)去掉，文件末尾帶井號(hào)的行去掉。

利用R語(yǔ)言將注釋結(jié)果整理成enricher函數(shù)需要的輸入格式

library(stringr)
library(dplyr)
egg<-read.table("out.emapper.annotations",sep="\t",header=T)
egg[egg==""]<-NA
gterms <- egg %>%
 select(query_name, GOs) %>%
 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)
}
head(gene2go)

> head(gene2go)
# A tibble: 6 x 2
 gene           term
 <chr>          <chr>
1 NP_001018179.1 GO:0003674
2 NP_001018179.1 GO:0003824
3 NP_001018179.1 GO:0004418
4 NP_001018179.1 GO:0005575
5 NP_001018179.1 GO:0005622
6 NP_001018179.1 GO:0005623

獲得一個(gè)兩列的數(shù)據(jù)框，有了這個(gè)數(shù)據(jù)框就可以做GO富集分析了

在 https://www.jianshu.com/p/9c9e97167377 這篇文章里的評(píng)論區(qū)有人提到上面用到的for循環(huán)代碼效率比較低，他提供的代碼是

gene_ids <- egg$query_name
eggnog_lines_with_go <- egg$GOs!= ""
eggnog_lines_with_go
eggnog_annoations_go <- str_split(egg[eggnog_lines_with_go,]$GOs, ",")
gene_to_go <- data.frame(gene = rep(gene_ids[eggnog_lines_with_go],
                                  times = sapply(eggnog_annoations_go, length)),
                        term = unlist(eggnog_annoations_go))
head(gene_to_go)

> head(gene_to_go)
           gene       term
1 NP_001018179.1 GO:0003674
2 NP_001018179.1 GO:0003824
3 NP_001018179.1 GO:0004418
4 NP_001018179.1 GO:0005575
5 NP_001018179.1 GO:0005622
6 NP_001018179.1 GO:0005623

用這個(gè)代碼替換for循環(huán)，確實(shí)快好多。

首先準(zhǔn)備一個(gè)基因列表，我這里選取gene2go中的前40個(gè)基因作為測(cè)試 還需要為TERM2GENE=參數(shù)準(zhǔn)備一個(gè)數(shù)據(jù)框，第一列是term，第二列是基因ID，只需要把gene2go的列調(diào)換順序就可以了。

library(clusterProfiler)
gene_list<-gene2go$gene[1:40]
term2gene<-gene2go[,c(2,1)]
df<-enricher(gene=gene_list,
            pvalueCutoff = 0.05,
            pAdjustMethod = "BH",
            TERM2GENE = term2gene)
head(df)
barplot(df)
dotplot(df)

y軸的標(biāo)簽通常是GO term (就是文字的那個(gè))而不是GO id。clusterProfiler包同樣提供了函數(shù)對(duì)ID和term互相轉(zhuǎn)換。go2term()go2ont()

df<-as.data.frame(df)
head(df)
dim(df)
df1<-go2term(df$ID)
dim(df1)
head(df1)
df$term<-df1$Term
df2<-go2ont(df$ID)
dim(df2)
head(df2)
df$Ont<-df2$Ontology
head(df)
df3<-df%>%
 select(c("term","Ont","pvalue"))
head(df3)
library(ggplot2)
ggplot(df3,aes(x=term,y=-log10(pvalue)))+
 geom_col(aes(fill=Ont))+
 coord_flip()+labs(x="")+
 theme_bw()

 

這里遇到一個(gè)問(wèn)題：數(shù)據(jù)框如何分組排序？目前想到一個(gè)比較麻煩的辦法是將每組數(shù)據(jù)弄成一個(gè)單獨(dú)的數(shù)據(jù)框，排好序后再合并。

library(stringr)
library(dplyr)
library(clusterProfiler)
egg<-read.table("out.emapper.annotations",sep="\t",header=T)
egg[egg==""]<-NA
gene2ko <- egg %>%
 dplyr::select(GID = query_name, Ko = KEGG_ko) %>%
 na.omit()
head(gene2ko)
head(gene2go)
gene2ko[,2]<-gsub("ko:","",gene2ko[,2])
head(gene2ko)
#kegg_info.RData這個(gè)文件里有pathway2name這個(gè)對(duì)象
load(file = "kegg_info.RData")
pathway2gene <- gene2ko %>% left_join(ko2pathway, by = "Ko") %>%
 dplyr::select(pathway=Pathway,gene=GID) %>%
 na.omit()
head(pathway2gene)
pathway2name
df<-enricher(gene=gene_list,
            pvalueCutoff = 0.05,
            pAdjustMethod = "BH",
            TERM2GENE = pathway2gene,
            TERM2NAME = pathway2name)
dotplot(df)
barplot(df)

以上最開(kāi)始的輸入文件是eggnog-mapper軟件本地版注釋結(jié)果，如果用在線(xiàn)版獲得的注釋結(jié)果，下載的結(jié)果好像沒(méi)有表頭，需要自己對(duì)應(yīng)好要選擇的列。