寫在前面

• 本文主要參考生信小白2018與多啦A夢(mèng)的時(shí)光機(jī)的簡(jiǎn)書，但由于EggNOG（v5.0）數(shù)據(jù)庫(kù)與eggnog-mapper（v2.1.3）更新后，修改了輸出文件的格式，因此他們二人寫的腳本會(huì)出現(xiàn)一些問(wèn)題，但R學(xué)的好的話完全可以克服，我保留了他們的大部分源碼并進(jìn)行了debug，優(yōu)化了一些步驟，供大家參考。
• 注意：完成下面的步驟至少需要~40GB運(yùn)行內(nèi)存與~50GB的存儲(chǔ)空間，因此沒(méi)有服務(wù)器的同學(xué)不要輕易嘗試??
• 當(dāng)然，如果你無(wú)論如何也構(gòu)建不成功可以和我koukou466478340，注意：并非無(wú)償哈

1. 使用eggnog-mapper生成注釋文件

• 以下代碼在Linux中運(yùn)行

1.1. 下載eggnog-mapper軟件

• 進(jìn)入你想存放的文件夾下下載。由于eggnog-mapper只托管在GitHub上，所以使用git下載，而不能用wget或pip或conda下載,托管到GitHub上的好處就是下載之后就是文件夾，不用解壓、編譯之類的操作。
• 2022年3月27日更新：conda上可以直接下載了 conda install -c bioconda eggnog-mapper
• 2022年10月3日更新：pip上可以直接下載了 pip install eggnog-mapper

git clone https://github.com/jhcepas/eggnog-mapper.git

1.2. 下載EggNOG數(shù)據(jù)庫(kù)

• 下載鏈接：http://eggnog5.embl.de/download/
• 請(qǐng)嚴(yán)格按照?qǐng)D片操作。

選擇最新的文件下載

• 注意：一定要放到eggnog-mapper軟件目錄下的data文件夾中，使用conda裝的放到類似的conda路徑下面/root/miniconda3/envs/python3/lib/python3.7/site-packages/data；其中需要更改的就是你安裝miniconda的路徑以及conda環(huán)境的名稱。因?yàn)槲沂前惭b在root/目錄下面，python3環(huán)境中的，所以路徑如上。
  

  下載這兩個(gè)文件

• 使用gunzip解壓

gunzip eggnog.db.gz eggnog.proteins.dmnd.gz

解壓

1.3. 將eggnog-mapper文件夾添加至環(huán)境變量【conda裝的可忽略】

• 注意：將路徑修改為你自己的

echo 'export PATH=/public/user/software/eggnog-mapper:$PATH' >> ~/.bashrc
source ~/.bashrc

1.4. 配置eggnog-mapper的運(yùn)行環(huán)境（可省略）

• 因?yàn)閑ggnog-mapper是運(yùn)行在python3（>v3.7）環(huán)境下的軟件，可以先使用conda建一個(gè)python3（教程）的環(huán)境，其實(shí)我們做生信的，常備不同的環(huán)境是很必要的

conda create -n python3 python=3.7.8 -y
activate python3

1.5. 運(yùn)行eggnog-mapper

注意：蛋白文件的中基因名必須要和你以后富集分析的基因名格式相同，如果不同，后面的操作都是徒勞，千萬(wàn)注意。

emapper.py -i mm39.pep.fa -o mm39 -d euk --cpu 10 --dbmem

• -i：輸入的蛋白組序列文件，請(qǐng)自行下載你物種的pep序列文件，并將其中的ID換為你常用的ID，不然到時(shí)候根據(jù)ID富集不到也就白忙活了。
• -o：輸出文件的前綴，后綴會(huì)自動(dòng)補(bǔ)充為.emapper.annotations
• -m：默認(rèn)就是diamond，因此不用指定。
• --dbmem:將數(shù)據(jù)庫(kù)全部存入運(yùn)行內(nèi)存，極大提高運(yùn)行速度，但特別吃運(yùn)存，大約要45GB的運(yùn)行空間。
• --cpu：使用的線程數(shù)
• --override：可以覆蓋原有的文件，多次嘗試的時(shí)候可以加上
• -d: 比對(duì)數(shù)據(jù)庫(kù)的類型，作者給了euk（真核）, bact（細(xì)菌）, arch（古菌）三大類。雖然大類更全，但有顧慮就是植物動(dòng)物微生物亂匹配，因此推薦使用--tax_scope
--tax_scope: 有以下參數(shù)可以選擇。
  如小鼠可選Eukaryota下面的Chordata（脊椎動(dòng)物），即emapper.py -i mm39.pep.fa -o mm39 --tax_scope Chordata
  如擬南芥可選Eukaryota下面的Viridiplantae（綠色植物），即emapper.py -i tair.pep.fa -o tair --tax_scope Viridiplantae
  再比如想選擇多個(gè)，那么就把他們寫入到文件中，然后入讀，即emapper.py -i tair.pep.fa -o tair --tax_scope ./example.txt
  更多分類也可以查看該目錄下的文件：/root/softwear/eggnog-mapper-master/eggnogmapper/annotation/tax_scopes

1.6. 漫長(zhǎng)的等待

• 根據(jù)pep序列文件的大小所需要的時(shí)間不同，我用小鼠的6萬(wàn)條序列，估計(jì)用了5個(gè)小時(shí)左右，當(dāng)然加入--dbmem的時(shí)間在20分鐘左右。這點(diǎn)時(shí)間，你可以大致閱讀下面的內(nèi)容。

1.7. 運(yùn)行結(jié)果

• 運(yùn)行結(jié)束后，eggnog-mapper會(huì)生成三個(gè)文件，以小鼠mm39為例:

結(jié)果文件

• mm39.emapper.hmm_hits: 記錄每個(gè)用于搜索序列對(duì)應(yīng)的所有的顯著性的eggNOG Orthologous Groups(OG). 所有標(biāo)記為"-"則表明該序列未找到可能的OG
• mm39.emapper.seed_orthologs: 記錄每個(gè)用于搜索序列對(duì)的的最佳的OG，也就是mm39.emapper.hmm_hits里選擇得分最高的結(jié)果。之后會(huì)從eggNOG中提取更精細(xì)的直系同源關(guān)系(orthology relationships)
• mm39.emapper.annotations: 該文件提供了最終的注釋結(jié)果，也是我們后續(xù)構(gòu)建OrgDb的文件
  它一共有21列，列名如下:
• query：檢索的基因名或者其他ID，就是提供蛋白組序列文件中的ID，所以這個(gè)ID和以后做GO/KEGG分析的ID要相同，不然就沒(méi)辦法做了，已經(jīng)說(shuō)了好多次了。
• seed_ortholog
• evalue
• score
• eggNOG_OGs
• max_annot_lvl
• COG_category
• Description
• Preferred_name：預(yù)測(cè)的基因名，特別類似AP2有含義的基因名，可以直接在NCBI 的 Gene數(shù)據(jù)庫(kù)中查看，一般是其他物種確定了功能的基因名
• GOs：推測(cè)的GO的詞條
• EC
• KEGG_ko：推測(cè)的KEGG KO詞條
• KEGG_Pathway：該條目顯示的是K中包含的Ko以及map，沒(méi)用，后面我們會(huì)自己生成
• KEGG_Module
• KEGG_Reaction
• KEGG_rclass
• BRITE
• KEGG_TC
• CAZy
• BiGG_Reaction：BiGG
• PFAMs

1.8. 處理數(shù)據(jù)

• 首先刪掉【mm39.emapper.annotations】開頭與末尾沒(méi)用的行和表頭的#即可，然后使用awk提取需要的列
  重輸出為【mm39.annotations】

id=mm39  #剛才的前綴
sed '/^##/d' ${id}.emapper.annotations| sed 's/#//g'| awk -vFS="\t" -vOFS="\t" '{print $1,$9,$10,$12}' > ${id}.annotations

• 然后就可以使用【mm39.annotations】開始在R中構(gòu)建OrgDb了

2. 構(gòu)建OrgDb

• 以下代碼在R中運(yùn)行

2.1. 安裝并導(dǎo)入所需R包

#install.packages("")
library(dplyr)
library(stringr)
library(jsonlite)
library(AnnotationForge)
#順手設(shè)置一下options
options(stringsAsFactors = F)

2.2. 讀入生成的annotations文件

emapper <- read.table("mm39.annotations", header=TRUE, sep = "\t",quote = "")
#將空值替換為NA，方便后續(xù)使用na.omit()函數(shù)提出沒(méi)有注釋到的行
emapper[emapper==""]<-NA

2.3. 提取GO信息

• 把emapper中的【query】列、【Preferred_name】列【GOs】列提取出來(lái)(其中【%>%】相當(dāng)于linux中的管道符【|】)

# library(dplyr)
gene_info <- emapper %>% dplyr::select(GID = query, GENENAME = Preferred_name) %>% na.omit()
gos <- emapper %>% dplyr::select(query, GOs) %>% na.omit()

• 構(gòu)建一個(gè)空的【gene2go】數(shù)據(jù)框，為后面填充數(shù)據(jù)用

gene2go = data.frame(GID = character(),
                     GO = character(),
                     EVIDENCE = character())

• 對(duì)【gene2go】空數(shù)據(jù)框進(jìn)行填充。

# library(stringr)
gos_list <- function(x){
  the_gos <- str_split(x[2], ",", simplify = FALSE)[[1]]
  df_temp <- data.frame(GID = rep(x[1], length(the_gos)),
                        GO = the_gos,
                        EVIDENCE = rep("IEA", length(the_gos)))
  return(df_temp)
}
gene2gol <- apply(as.matrix(gos),1,gos_list)
gene2gol_df <- do.call(rbind.data.frame, gene2gol)
gene2go <- gene2gol_df
gene2go$GO[gene2go$GO=="-"]<-NA
gene2go<-na.omit(gene2go)

• 循環(huán)的作用簡(jiǎn)單來(lái)說(shuō)就單行變多行，舉例如下：

#                                               query1  GO:XXXXXX1
#query1 GO:XXXXXX1,GO:XXXXXX2,GO:XXXXXX3   =>   query1  GO:XXXXXX2
#                                               query1  GO:XXXXXX3

2.4. 同理，提取KEGG信息

• 把emapper中的query列、KEGG_ko列提取出來(lái)

gene2ko <- emapper %>% dplyr::select(GID = query, Ko = KEGG_ko)
gene2ko$Ko[gene2ko$Ko=="-"]<-NA
gene2ko<-na.omit(gene2ko)
gene2kol <- apply(as.matrix(gene2ko),1,gos_list)
gene2kol_df <- do.call(rbind.data.frame, gene2kol)
gene2ko <- gene2kol_df[,1:2]
colnames(gene2ko) <- c("GID","Ko")
gene2ko$Ko <- gsub("ko:","",gene2ko$Ko)

• 下載KO的json文件并【放到當(dāng)前文件夾下】，下載地址：https://www.genome.jp/kegg-bin/get_htext?ko00001
  
• 定義函數(shù)，不用管是啥直接跑，其中只用到了下載的【ko00001.json】文件，所以可以放心跑

# library(jsonlite)
# 下面的json = "ko00001.json"，如果你下載到其他地方，記得加上路徑
update_kegg <- function(json = "ko00001.json") {
  pathway2name <- tibble(Pathway = character(), Name = character())
  ko2pathway <- tibble(Ko = character(), Pathway = character())
  kegg <- fromJSON(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))))}}}
  save(pathway2name, ko2pathway, file = "kegg_info.RData")}
# 調(diào)用函數(shù)后在本地創(chuàng)建kegg_info.RData文件，以后只需要載入 "kegg_info.RData"即可
update_kegg()
# 載入kegg_info.RData文件
load(file = "kegg_info.RData")

• 根據(jù)gene2ko中的Ko信息將ko2pathway中的Pathway列提取出來(lái)

gene2pathway <- gene2ko %>% left_join(ko2pathway, by = "Ko") %>% dplyr::select(GID, Pathway) %>% na.omit()

2.5. 獲取物種信息

• 首先在NCBI上查詢你物種的相關(guān)信息，網(wǎng)址：https://www.ncbi.nlm.nih.gov/taxonomy
  

• 比如說(shuō)小鼠Mus musculus

tax_id = "10090"
genus = "Mus" 
species = "musculus"

• 或者人Homo sapiens

# tax_id = "9606"
# genus = "Homo" 
# species = "sapiens"

• 又或者擬南芥Arabidopsis thaliana

# tax_id = "3702"
# genus = "Arabidopsis" 
# species = "thaliana"

• 其實(shí)，都可以隨便填寫，不會(huì)影響建庫(kù)結(jié)果，只會(huì)影響你輸出文件的名字，像個(gè)性化輸出名字就自定義

• 可以去重，以防萬(wàn)一，或者等下面報(bào)錯(cuò)提示去重再去跑代碼也行。原理是，makeOrgPackage不允許有重復(fù)的行，因此需要?jiǎng)h除，除了gene2pathway可能有重復(fù)的行，其他幾乎不可能有重復(fù)的情況。

# gene2go <- unique(gene2go)
# gene2go <- gene2go[!duplicated(gene2go),]
# gene2ko <- gene2ko[!duplicated(gene2ko),]
# gene2pathway <- gene2pathway[!duplicated(gene2pathway),]
# gene_info <- gene_info[!duplicated(gene_info),]

2.6. 構(gòu)建OrgDb

• 開始構(gòu)建

makeOrgPackage(gene_info=gene_info,
               go=gene2go,
               ko=gene2ko,
               pathway=gene2pathway,
               version="1.34.1",  #版本，使用？makeOrgPackage，拉到最下面查看
               maintainer = "your name <郵箱>",  #修改為你的名字和郵箱
               author = "your name <郵箱>",  #修改為你的名字和郵箱
               outputDir = ".",  #輸出文件位置
               tax_id=tax_id,  #你在NCBI上查并定義的id
               genus=genus,  #你在NCBI上查并定義的屬名
               species=species,  #你在NCBI上查并定義的種名
               goTable="go")

• 構(gòu)建成功后會(huì)在當(dāng)前路徑下生成一個(gè)【文件夾】，名字形如【org.XXX.eg.db】(其中XXX的構(gòu)成為你定義的genus的大寫首字母+species名，如小鼠是Mmusculus，人是Hsapiens，擬南芥是Athaliana)

3. 導(dǎo)入OrgDb庫(kù)

install.packages("./org.Mmusculus.eg.db", repos=NULL)
library(org.Mmusculus.eg.db)
# 查看所有列的信息
columns(org.Mmusculus.eg.db)
# 查看所有基因
keys(org.Mmusculus.eg.db)
# 查看特定基因的信息
# library(dplyr)
select(org.Mmusculus.eg.db, keys = "CW07G09620", columns = c("GO"))

4. 導(dǎo)出背景文件

• 因?yàn)槲覀冊(cè)诤竺孢M(jìn)行kegg分析的時(shí)候需要背景文件，因此為了方便即可導(dǎo)出pathway2name和pathway2gene

# 做一些常規(guī)格式化
pathway2name$Name <- gsub(" \\[BR:ko[0-9]{5}\\]", "",pathway2name$Name)
pathway2name<- na.omit(pathway2name)
pathway2gene <-gene2pathway[, c("Pathway","GID")]
# 輸出
write.table(pathway2name,file = "./pathway2name", sep = "\t", quote = F, row.names = F)
write.table(pathway2gene,file = "./pathway2gene", sep = "\t", quote = F, row.names = F)

5. 使用clusterProfiler富集分析

5.1 KEGG富集分析（不需要OrgDB）

library(clusterProfiler)
# 每次只需導(dǎo)入下面兩個(gè)文件即可
pathway2gene <- read.table("./pathway2gene",header = T,sep = "\t")
pathway2name <- read.table("./pathway2name",header = T,sep = "\t")
# 導(dǎo)入你鑒定到的差異基因列表，并轉(zhuǎn)化為向量
gene <- read.csv("/root/total_diff_gene.csv")
gene_list <- gene[,1]
# KEGG pathway 富集
ekp <- enricher(gene_list, 
                TERM2GENE = pathway2gene, 
                TERM2NAME = pathway2name, 
                pvalueCutoff = 1,  # 表示全部保留，可以設(shè)為0.05作為閾值
                qvalueCutoff = 1, # 表示全部保留，可以設(shè)為0.05作為閾值
                pAdjustMethod = "BH",
                minGSSize = 1)
dotplot(ekp)

5.2 GO富集分析（需要OrgDB）

library(org.Mmusculus.eg.db)
ego <- enrichGO(gene=gene_list,
                OrgDb=org.Mmusculus.eg.db,
                keyType="GID",
                ont="ALL",   #CC/BP/MF可選
                qvalueCutoff = 0.05,
                pvalueCutoff =0.05)
dotplot(ego)
# 無(wú)論是ekp還是ego都可以選擇將數(shù)據(jù)導(dǎo)出，然后可以自己使用ggplot畫
ego_results<-as.data.frame(ego)
write.table(ego_results, file = "./ego_results.txt", quote = F)

image.png