有時候不得不感慨，時間已經(jīng)來到了2024年的尾聲，但是對于許多非模式的物種來說，各種工具、數(shù)據(jù)集、參考等等依然還在寒冬。今天帶大家手動的構(gòu)建屬于小眾物種的富集分析全攻略。
1. 簡介
相信每一個做過生信的人來說，對于GO/KEGG等富集分析并不陌生。有越來越多的工具和云平臺可以提供這些分析。然而，對于許多對于許多農(nóng)口或者微生物研究人員來說，由于低質(zhì)量的參考基因組，和無人構(gòu)建的數(shù)據(jù)集都使得這一“簡單”的需求莫名的走了許多彎路。
其實(shí)，簡單來說，對于任何一個物種都需要對應(yīng)的基因->條目Term這樣的關(guān)系表，有了這樣的關(guān)系就可以構(gòu)建我們的抽獎模型“超分布幾何檢驗(yàn)”。只不過對于大多數(shù)的模式物種來說，有大量的前輩已經(jīng)把這些整理好了；前人栽樹，后人乘涼即可。所以，理論上當(dāng)我們沒有現(xiàn)成的數(shù)據(jù)庫時，就需要我們構(gòu)建自己的物種的注釋數(shù)據(jù)庫，就可以輕松的完成今后的富集分析。
2.注釋庫構(gòu)建
注釋庫的構(gòu)建我們使用eggNOG，eggNOG（evolutionary genealogy of genes: Non-supervised Orthologous Groups）是一種用于基因組功能注釋的數(shù)據(jù)庫和工具，廣泛應(yīng)用于基因組學(xué)和生物信息學(xué)研究。目前已經(jīng)更新到了eggNOG v5版本，他可以一次性的完整GO/COG/KEGG/PFAM等一系列的注釋。
2.1 準(zhǔn)備工作
1. 全轉(zhuǎn)錄本/基因的核酸或蛋白序列： 通常為fasta格式，里面包含了我們物種的所有基因/轉(zhuǎn)錄本。當(dāng)然如果是組裝程度極低或者甚至無參考的，可能是一些contig/scaffold。
2.eggNOG-maaper：該工具可以在github上下載https://github.com/eggnogdb/eggnog-mapper/releases/
3.eggNOG數(shù)據(jù)庫：這部份建議手動下載，目前最新版本的地址如下：
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 比對數(shù)據(jù)庫
安裝好eggNOG-maaper，和對應(yīng)的python依賴后將我們下載的數(shù)據(jù)庫解壓到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 &
如果我們的是核酸序列則需要先翻譯，此時使用命令nohup emapper.py -i gene.fa --itype genome -o myeggnog --translate --cpu 12 --dbmem -m diamond &
具體的其他參數(shù)可以通過-h進(jìn)行查看，cpu數(shù)可以根據(jù)實(shí)際調(diào)整。
2.3 整理結(jié)果
當(dāng)emapper比對完畢后會得到 myeggnog.emapper.annotations文件，該文件就是我們后面用于構(gòu)建我們自己GO/KEGG數(shù)據(jù)庫的關(guān)鍵，我們需要對結(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)行完畢后會在目錄下生成用于GO富集的Org會根據(jù)genus和species生成，這里就是org.Gspecies.eg.db，使用install.packages("org.Gspecies.eg.db",repos = NULL, type = "source")安裝我們的包，后續(xù)分析時候只需要參考正常的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的分析則需要我們手動使用 enricher進(jìn)行，關(guān)鍵代碼如下：
library(tidyverse)
library(clusterProfiler)
kegg_df <- read_tsv('MyKEGG.xls')
#這里用于測試我們?nèi)∏?00個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
至此，我們就可以對任意的物種輕松的進(jìn)行GO/KEGG分析，一起動手試試吧！
