有時(shí)候不得不感慨档叔，時(shí)間已經(jīng)來(lái)到了2024年的尾聲，但是對(duì)于許多非模式的物種來(lái)說(shuō)芋膘，各種工具鳞青、數(shù)據(jù)集、參考等等依然還在寒冬为朋。今天帶大家手動(dòng)的構(gòu)建屬于小眾物種的富集分析全攻略臂拓。
1. 簡(jiǎn)介
相信每一個(gè)做過(guò)生信的人來(lái)說(shuō)，對(duì)于GO/KEGG等富集分析并不陌生潜腻。有越來(lái)越多的工具和云平臺(tái)可以提供這些分析埃儿。然而，對(duì)于許多對(duì)于許多農(nóng)口或者微生物研究人員來(lái)說(shuō)融涣，由于低質(zhì)量的參考基因組童番，和無(wú)人構(gòu)建的數(shù)據(jù)集都使得這一“簡(jiǎn)單”的需求莫名的走了許多彎路精钮。
其實(shí)，簡(jiǎn)單來(lái)說(shuō)剃斧，對(duì)于任何一個(gè)物種都需要對(duì)應(yīng)的基因->條目Term這樣的關(guān)系表轨香，有了這樣的關(guān)系就可以構(gòu)建我們的抽獎(jiǎng)模型“超分布幾何檢驗(yàn)”。只不過(guò)對(duì)于大多數(shù)的模式物種來(lái)說(shuō)幼东，有大量的前輩已經(jīng)把這些整理好了臂容；前人栽樹，后人乘涼即可根蟹。所以脓杉，理論上當(dāng)我們沒(méi)有現(xiàn)成的數(shù)據(jù)庫(kù)時(shí)，就需要我們構(gòu)建自己的物種的注釋數(shù)據(jù)庫(kù)简逮，就可以輕松的完成今后的富集分析球散。
2.注釋庫(kù)構(gòu)建
注釋庫(kù)的構(gòu)建我們使用eggNOG，eggNOG（evolutionary genealogy of genes: Non-supervised Orthologous Groups）是一種用于基因組功能注釋的數(shù)據(jù)庫(kù)和工具散庶，廣泛應(yīng)用于基因組學(xué)和生物信息學(xué)研究蕉堰。目前已經(jīng)更新到了eggNOG v5版本，他可以一次性的完整GO/COG/KEGG/PFAM等一系列的注釋悲龟。
2.1 準(zhǔn)備工作
1. 全轉(zhuǎn)錄本/基因的核酸或蛋白序列： 通常為fasta格式屋讶，里面包含了我們物種的所有基因/轉(zhuǎn)錄本。當(dāng)然如果是組裝程度極低或者甚至無(wú)參考的须教，可能是一些contig/scaffold皿渗。
2.eggNOG-maaper：該工具可以在github上下載https://github.com/eggnogdb/eggnog-mapper/releases/
3.eggNOG數(shù)據(jù)庫(kù)：這部份建議手動(dòng)下載，目前最新版本的地址如下：
http://eggnog6.embl.de/download/emapperdb-5.0.2/eggnog_proteins.dmnd.gz
http://eggnog6.embl.de/download/emapperdb-5.0.2/mmseqs.tar.gz
http://eggnog6.embl.de/download/emapperdb-5.0.2/eggnog.db.gz
http://eggnog6.embl.de/download/emapperdb-5.0.2/eggnog.taxa.tar.gz
http://eggnog6.embl.de/download/emapperdb-5.0.2/pfam.tar.gz
2.2 比對(duì)數(shù)據(jù)庫(kù)
安裝好eggNOG-maaper没卸，和對(duì)應(yīng)的python依賴后將我們下載的數(shù)據(jù)庫(kù)解壓到eggNOG-mapper的data目錄下羹奉，結(jié)構(gòu)大致如下：
data/
├── eggnog.db
├── eggnog_proteins.dmnd
├── eggnog.taxa.db
├── eggnog.taxa.db.traverse.pkl
├── mmseqs
│ ├── mmseqs.db
│ ├── mmseqs.db.dbtype
│ ├── mmseqs.db_h
│ ├── mmseqs.db_h.dbtype
│ ├── mmseqs.db_h.index
│ ├── mmseqs.db.index
│ ├── mmseqs.db.lookup
│ └── mmseqs.db.source
└── pfam
├── Pfam-A.clans.tsv.gz
├── Pfam-A.hmm
├── Pfam-A.hmm.h3f
├── Pfam-A.hmm.h3i
├── Pfam-A.hmm.h3m
├── Pfam-A.hmm.h3m.ssi
├── Pfam-A.hmm.h3p
└── Pfam-A.hmm.idmap
如果我們的輸入是蛋白序列則使用nohup emapper.py -i proteins.fa -o myeggnog --cpu 12 --dbmem -m diamond &
如果我們的是核酸序列則需要先翻譯，此時(shí)使用命令nohup emapper.py -i gene.fa --itype genome -o myeggnog --translate --cpu 12 --dbmem -m diamond &
具體的其他參數(shù)可以通過(guò)-h進(jìn)行查看约计，cpu數(shù)可以根據(jù)實(shí)際調(diào)整诀拭。
2.3 整理結(jié)果
當(dāng)emapper比對(duì)完畢后會(huì)得到 myeggnog.emapper.annotations文件，該文件就是我們后面用于構(gòu)建我們自己GO/KEGG數(shù)據(jù)庫(kù)的關(guān)鍵煤蚌，我們需要對(duì)結(jié)果進(jìn)行進(jìn)一步的整理耕挨。
library(tidyverse)
egg <- read_tsv('myeggnog.emapper.annotations', comment = '##') %>% mutate( `#query` = str_replace( `#query`, "_\\d+$", "")) %>% group_by(`#query`) %>% top_n( n=1, wt = score)
egg_df <- egg %>% select(GID = `#query`, GOs, KEGG_ko) %>% filter(GOs != '-', KEGG_ko != '-') %>% separate_rows(GOs, sep = ",") %>% separate_rows(KEGG_ko, sep = ",")
egg_df <- egg_df %>% mutate(KEGG_ko = str_replace_all(KEGG_ko,'ko:','')) %>% rename(KO = KEGG_ko) %>% distinct()
library(KEGGREST)
kegg_pathways <- keggList("pathway", "ko") # "ko" 表示 KO 專用的 pathway
kegg_ko_pathways <- keggLink("ko", "pathway")
kegg_db <- tibble(Pathway = names(kegg_pathways), Name = kegg_pathways)
kegg_db <- tibble(Ko = sub("ko:","", kegg_ko_pathways), Pathway = sub("path:", "", names(kegg_ko_pathways))) %>% filter(str_detect(Pathway, '^ko')) %>% inner_join(kegg_db, by = 'Pathway')
gene_info <- egg %>% select(GID = `#query`, GENENAME = seed_ortholog)
gene2go <- egg_df %>% select(GID, GOs) %>% mutate(EVIDENCE = 'IEA') %>% distinct() %>% rename(GO = GOs)
koterms <- egg_df %>% select(GID, KO) %>% distinct()
gene2pathway <- kegg_db %>% inner_join(koterms, by = c('Ko' = 'KO')) %>% select(GID, Pathway) %>% distinct()
library(AnnotationForge)
makeOrgPackage(gene_info=gene_info, go=gene2go, ko=koterms, pathway=gene2pathway, version="1.0.0", maintainer='Terry Xizzy <xxx@gmail.com>', author='Terry Xizzy <xxx@gmail.com>',outputDir=".", tax_id="12345", genus="Genus", species="species",goTable="go")
kegg_df <- koterms %>% inner_join(kegg_db, by = c('KO' = 'Ko')) %>% write_tsv('MyKEGG.xls')
在實(shí)際運(yùn)行中只需要修改自己的輸入文件名即可。這里的結(jié)果是基于5.0.2構(gòu)建尉桩，后續(xù)可能出現(xiàn)列名變化的情況筒占，需要根據(jù)實(shí)際進(jìn)行調(diào)整。
代碼運(yùn)行完畢后會(huì)在目錄下生成用于GO富集的Org會(huì)根據(jù)genus和species生成蜘犁，這里就是org.Gspecies.eg.db翰苫，使用install.packages("org.Gspecies.eg.db",repos = NULL, type = "source")安裝我們的包，后續(xù)分析時(shí)候只需要參考正常的clusterProfiler方法，詳細(xì)的方法可以參考我以前的文章奏窑，核心步驟如下：
library(org.Gspecies.eg.db))
Db <- org.Gspecies.eg.db
enrichGO(gene_list, OrgDb = Db, ont='ALL',pAdjustMethod = 'BH',pvalueCutoff = 1,qvalueCutoff = 1,keyType = 'SYMBOL')
這里有一點(diǎn)需要注意的是导披，我們構(gòu)建數(shù)據(jù)使用的GID應(yīng)該和我們輸入的gene_list的使用同一種。
而KEGG的分析則需要我們手動(dòng)使用 enricher進(jìn)行埃唯，關(guān)鍵代碼如下：
library(tidyverse)
library(clusterProfiler)
kegg_df <- read_tsv('MyKEGG.xls')
#這里用于測(cè)試我們?nèi)∏?00個(gè)GID作為gene_list進(jìn)行分析
gene_list <- kegg_df$GID[1:100]
kegg_res <- enricher(gene_list, TERM2GENE = select(kegg_df, Pathway, GID), TERM2NAME = select(kegg_df, Pathway, Name), pvalueCutoff = 1,qvalueCutoff = 1)
kegg_res %>% as_tibble
輸出如下：
kegg_res %>% as_tibble
# A tibble: 134 × 9
ID Description GeneRatio BgRatio pvalue p.adjust qvalue geneID Count
<chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <chr> <int>
1 ko01120 Microbial … 70/81 162/46… 3.75e-97 5.02e-95 3.04e-95 unass… 70
2 ko00010 Glycolysis… 46/81 46/4637 1.58e-88 1.06e-86 6.41e-87 unass… 46
3 ko01110 Biosynthes… 75/81 441/46… 1.05e-71 4.69e-70 2.84e-70 unass… 75
4 ko01200 Carbon met… 51/81 113/46… 9.27e-67 3.10e-65 1.88e-65 unass… 51
至此撩匕，我們就可以對(duì)任意的物種輕松的進(jìn)行GO/KEGG分析，一起動(dòng)手試試吧墨叛！
