非模式生物的GO和KEGG富集

最近在做非模式生物的GO和KEGG富集分析，参考了网上的一些帖子和知乎专栏，发现代码总有一些小问题，于是自己摸索修改终于跑通了= =，这里做个记录。

本篇笔记代码主要参考自：

比较基因组学分析3：特异节点基因家族富集分析（非模式物种GO/KEEG富集分析） - 简书 (jianshu.com)

小众或非模式生物的自建库GO/KEGG富集分析 - 知乎 (zhihu.com)

1. eggNOG注释

因为我们要做的物种是非模式生物，所有背景基因需要我们亲自注释，这里只推荐eggNOG，速度快，可以本地运行或者在线运行，在线运行有10万条序列的限制。在线网页的介绍可以看我的这篇笔记：基因组注释（6）——在线版eggNOG-mapper注释功能基因 - 我的小破站 (shelven.com)

假设我们已经拿蛋白序列做了注释，在一大堆结果文件中我们只需要out.emapper.annotation。用notepad++打开，去掉带“#”符号的前三行，把表头的query前面的“#”注释也给去掉，这些处理做好后文件如下：

后续都是从这个文件中找到基因名与GO和KEGG号之间的对映关系。

假设我们从基因家族收缩/扩张分析、转录组等等数据已经拿到了我们想要的序列，创建一个新文件，一行一个基因名放入我们要分析的序列，文件命名为gene.txt如下所示：

2. GO富集

首先进行GO注释库的构建：

# 下载最新的GO库注释
wget http://purl.obolibrary.org/obo/go/go-basic.obo
# 处理注释文件，主要是帮我们找对应关系
grep "^id:" go-basic.obo |awk '{print $2}' > GO.id
grep "^name:" go-basic.obo | awk -F ': ' '{print $2}' > GO.name
grep "^namespace:" go-basic.obo |awk '{print $2}' > GO.class
paste GO.id GO.name GO.class -d "\t" > GO.library

构建的GO注释库如下：

第一列是GO号，第二列是描述（注意要完整的描述，而非一个单词，否则后续因为描述重复而无法富集！），第三列是GO分类。

library(clusterProfiler)
library(dplyr)
library(stringr)

setwd("D:\\zhuomian\\基因家族进化\\GO")	# 设置自己的工作目录

## GO注释生成
options(stringsAsFactors = F)
egg <- read.delim("out.emapper.annotations",header = T,sep="\t")
egg[egg==""]<-NA
gterms <- egg %>%
  dplyr::select(query, GOs) %>% na.omit()
gene_ids <- egg$query
eggnog_lines_with_go <- egg$GOs!= ""
eggnog_annoations_go <- str_split(egg[eggnog_lines_with_go,]$GOs, ",")
gene2go <- data.frame(gene = rep(gene_ids[eggnog_lines_with_go],	# 一个基因可能有多个GOterm，需要拆分成1：1对应关系
                      times = sapply(eggnog_annoations_go, length)),
                      term = unlist(eggnog_annoations_go))
go2name <- read.delim('GO.library', header = FALSE, stringsAsFactors = FALSE)
names(go2name) <- c('ID', 'Description', 'Ontology')

# GO富集
gene_select <- read.delim(file = 'gene.txt', stringsAsFactors = FALSE,header = F)$V1
go_rich <- enricher(gene = gene_select,
                    TERM2GENE = gene2go[c('term','gene')], 	# 这两项都要注意顺序
                    TERM2NAME = go2name[c('ID', 'Description')], 
                    pvalueCutoff = 1, 
                    pAdjustMethod = 'BH', 
                    qvalueCutoff = 1
)
tmp <- merge(go_rich, go2name[c('ID', 'Ontology')], by = 'ID')
tmp <- tmp[c(10, 1:9)]
tmp <- tmp[order(tmp$pvalue), ]
write.table(tmp, 'GO_enrichment.xls', sep = '\t', row.names = FALSE, quote = FALSE)

看一下生成的GO_enrichment.xls文件如下：

我们一般拿到这个文件后，根据qvalue值从小到大筛选可信度最高的富集结果，再手动去掉不合理的（比如你做植物，却注释到动物和微生物的）。

做到这一步生成GO富集文件就可以了~~（再做就不礼貌了）~~，GO富集结果可以在R做可视化，但是我觉得没必要，真的，调参数真的太难受了，而且每个人作图的要求都不一样。我个人比较建议在一些在线网站上手动画，比如这个网站：ChiPlot

网站的开发者在B站传过可视化富集结果的教学视频：【ChiPlot】绘制KEGG/GO富集图_哔哩哔哩_bilibili

省流不看版

根据富集结果，筛选molecular function、biological process和cellular component三个类别中qvalue值前5的GO term，计算第四列的GeneRatio值。

创建一个新的excel表，需要5列内容，如下：

type那一列不需要，只是为了方便展示，id这一列的列名可以随意，Description、Qvalue、Count、GeneRatio四列列名必须要一样，顺序无所谓。

再做一个展示图层的excel，作用是区分GO term属于哪个类型：

传到上面网站的KEGG/GO富集作图的模块，手动调整参数（自己多试一下）再加点AI作图，反正比我R做的好一些：

3. KEGG富集

思路和GO注释是类似的，相比GO，KEGG可选择的物种会多一些，我们这里仍然当作是非模式生物，以第一步eggNOG注释的基因集为背景基因集。

构建KEGG注释库：

# 下载KEGG注释库 https://www.genome.jp/kegg-bin/get_htext?ko00001，点击 Download json，下载得到ko00001.json文件

library(clusterProfiler)
library(dplyr)
library(stringr)
library(jsonlite)

setwd("D:\\zhuomian\\基因家族进化\\KEGG")

pathway2name <- tibble(Pathway = character(), Name = character())
ko2pathway <- tibble(Ko = character(), Pathway = character())
kegg <- fromJSON("ko00001.json")

for (a in seq_along(kegg[["children"]][["children"]])) {
  A <- kegg[["children"]][["name"]][[a]]
  for (b in seq_along(kegg[["children"]][["children"]][[a]][["children"]])) {
    B <- kegg[["children"]][["children"]][[a]][["name"]][[b]] 
    for (c in seq_along(kegg[["children"]][["children"]][[a]][["children"]][[b]][["children"]])) {
      pathway_info <- kegg[["children"]][["children"]][[a]][["children"]][[b]][["name"]][[c]]
      pathway_id <- str_match(pathway_info, "ko[0-9]{5}")[1]
      pathway_name <- str_replace(pathway_info, " \\[PATH:ko[0-9]{5}\\]", "") %>% str_replace("[0-9]{5} ", "")
      pathway2name <- rbind(pathway2name, tibble(Pathway = pathway_id, Name = pathway_name))
      kos_info <- kegg[["children"]][["children"]][[a]][["children"]][[b]][["children"]][[c]][["name"]]
      kos <- str_match(kos_info, "K[0-9]*")[,1]
      ko2pathway <- rbind(ko2pathway, tibble(Ko = kos, Pathway = rep(pathway_id, length(kos))))
    }
  }
}
colnames(ko2pathway) <- c("KO","Pathway")

write.table(pathway2name,"KEGG.library",sep="\t",row.names = F)

构建的KEGG注释库如下：

第一列是ko号，第二列是ko号对应的代谢通路描述，同样的，需要注意一个基因可能对应多个ko号，在对映关系上需要做一些处理：

## KEGG注释生成
options(stringsAsFactors = F)
egg <- read.delim("out.emapper.annotations",header = T,sep="\t")
egg[egg==""]<-NA
gene2ko <- egg %>%
  dplyr::select(GID = query, KO = KEGG_ko) %>%
  na.omit()

pathway2name <- read.delim("KEGG.library")
colnames(pathway2name)<-c("Pathway","Name")

ko2gene <- tibble(Ko=character(),GID=character())
for (Query in gene2ko$GID){
  ko_list <- strsplit(gene2ko$KO[which(gene2ko[,1]==Query)],split = ',')
  for (ko in ko_list){
    if (length(which(ko2gene[,1]==ko))==0){
      tmp <- data.frame(Ko=ko,GID=Query)
      ko2gene <- rbind(ko2gene,tmp)
    }
    else{
      old_Query <- ko2gene$GID[which(ko2gene[,1]==ko)]
      ko2gene$GID[which(ko2gene[,1]==ko)] <- paste(old_Query,Query,sep = ',')
    }
  }
}

pathway2gene <- tibble(Pathway = character(), GID = character())

for (ko in ko2pathway$KO){
  pathway_list <- ko2pathway$Pathway[which(ko2pathway[,1]==ko)]
  for (pathway in pathway_list){
    if (paste('ko:',ko,sep='') %in% ko2gene$Ko){
      ko <- paste('ko:',ko,sep='')
      if (length(which(pathway2gene[,1]==pathway))==0 ){
        ko2gene$GID[which(ko2gene[,1]==ko)]
        tmp <- data.frame(pathway=pathway,GID=ko2gene$GID[which(ko2gene[,1]==ko)])
        pathway2gene <- rbind(pathway2gene,tmp)
      }
      else{
        old_Query <- pathway2gene$GID[which(pathway2gene[,1]==pathway)]
        Query <- ko2gene$GID[which(ko2gene[,1]==ko)]
        pathway2gene$GID[which(pathway2gene[,1]==pathway)] <- paste(old_Query,Query,sep=',')
      }
    }
  }
}
	
new_pathway2gene <- data.frame()	# 遍历pathway2gene每一行，拆分gene

for (i in 1:nrow(pathway2gene)) {
  pathway <- pathway2gene$pathway[i]
  genes <- strsplit(as.character(pathway2gene$GID[i]), ",")[[1]]
  for (gene in genes) {
    new_row <- data.frame(pathway = pathway, gene = gene)
    new_pathway2gene <- rbind(new_pathway2gene, new_row)
  }
}

#KEGG富集
gene_select <- read.delim ('gene.txt', stringsAsFactors = FALSE,header = F)$V1
kegg_rich <- enricher (gene = gene_select,
                       TERM2GENE = new_pathway2gene[c ('pathway','gene')], 
                       TERM2NAME = pathway2name[c ('Pathway','Name')], 
                       pvalueCutoff = 1, 
                       pAdjustMethod = 'BH', 
                       qvalueCutoff = 1, 
                       maxGSSize = 500)

write.table (kegg_rich, 'KEGG_enrichment.xls', sep = '\t', row.names = FALSE, quote = FALSE)

最终得到的KEGG富集结果KEGG_enrichment.xls如下：

同样是根据qvalue筛选，剔除不合理的富集通路（比如你做的植物，富集到动物、微生物等代谢通路），作图是和GO富集作图一模一样的，这里不赘述了。

剔除不合理的代谢通路

顺便说一下，KEGG结果中会有不少不合理的地方，比如我注释的植物基因代谢通路，会有挺多注释到非植物相关的KEGG通路……

当然，我们可以一个个在KEGG官网查询富集的通路信息，也可以在R中把KEGG的所有植物代谢通路整理出来。这里需要用到一个R包KEGGREST快速获取KEGG数据库信息，需要注意，通过BiocManager下载的KEGGREST包已经不能正常使用了（2022年6月1日调用API地址发生了变化），需要在github上下载最新版：

# 国内或许会出现网络问题，需要一点魔法（摊手）
## 不要找别人备份的库，可能也会有问题，认准官方Bioconductor
library(devtools)
devtools::install_github("https://github.com/Bioconductor/KEGGREST")
library(KEGGREST)

参考代码：玩转KEGG (二)——植物富集到了动物通路？难不成咱研究的植物人儿😂 - 知乎 (zhihu.com)

具体思路是用KEGGREST包抓取植物物种信息，获取每个KEGG库中植物的代谢通路，取并集：

org <- data.frame(keggList("organism"))
plants <- org[grep("Plants", org$phylogeny), ]# 选取植物的物种信息
pathways_tot <- vector()

for (i in 1:length(plants$organism)) {
  try({
    pathways <- keggLink("pathway", plants[i,2])
    pathways <- sub(paste(".*",plants[i,2], sep = ""), "", pathways)
    pathways <- unique(pathways)
    pathways_tot <- append(pathways_tot,pathways)
    pathways_tot <- unique(pathways_tot) })
}

pathways_tot <- paste0("ko", pathways_tot)
writeLines(pathways_tot, "plants.kegg.txt")

每一个物种都要抓取ko号，所以时间会有点慢，大概10分钟左右。

可以看到现在KEGG中一共有152个植物，取代谢通路的并集后，每行输出结果。

最后用我熟悉的python简单处理过滤一下：

with open('plants.kegg.txt', 'r') as file:
    ko_plant = file.read().splitlines()

with open('KEGG_enrichment.xls','r') as input, open('KEGG_enrichment.filter.xls','w') as output:
    for line in input:
        if line.startswith('ko'):
            id = line.split('\t')[0]
            if id in ko_plant:
                output.write(line)
        else:
            output.write(line)