Introduction
TaxonKit是一個用于處理生物分類學(xué)數(shù)據(jù)的命令行工具悯恍。
它的主要功能是處理NCBI的生物分類學(xué)數(shù)據(jù)锐帜，包括對分類單元（如物種、屬验毡、科等）的查找、分類單元的上下位關(guān)系查詢典蜕、分類單元名稱的標準化等。
為了方便R社區(qū)用戶（自己）使用和流程整合，我把Taxonkit工具整合進了R包pctax，也開發(fā)了一些配套的系統(tǒng)發(fā)育分析和可視化方法斋射。
R調(diào)用Taxonkit
準備工作
- 安裝
pctax
pctax穩(wěn)定版本可在CRAN上獲得:
install.packages("pctax")
或者你可以通過以下方式從GitHub安裝pctax的開發(fā)版本:
# install.packages("devtools")
devtools::install_github("Asa12138/pctax")
- 安裝taxonkit：
library(pctax)
pctax::install_taxonkit(make_sure = TRUE)
#成功后taxonkit會安裝在下面這個目錄??
tools::R_user_dir("pctax")
- 下載NCBI Taxonomy數(shù)據(jù)文件:
pctax::download_taxonkit_dataset(make_sure = TRUE)
#成功后Taxonomy數(shù)據(jù)文件會在下面這個目錄??
file.path(Sys.getenv("HOME"), ".taxonkit")
該函數(shù)會下載官網(wǎng)最新版本的Taxonomy數(shù)據(jù)庫，如果需要制定版本的數(shù)據(jù)庫但荤，可以自己在官網(wǎng)下載：https://ftp.ncbi.nih.gov/pub/taxonomy/绩鸣，然后指定位置：
pctax::download_taxonkit_dataset(make_sure = TRUE,taxdump_tar_gz = "~/Downloads/taxdump.tar.gz")
使用
# 下列命令不報錯說明可以正常使用
check_taxonkit(print = FALSE)
主要功能與taxonkit一致：
| 函數(shù) | 功能 |
|---|---|
taxonkit_list |
列出指定TaxId下所有子單元的的TaxID |
taxonkit_lineage |
根據(jù)TaxID獲取完整譜系（lineage） |
taxonkit_reformat |
將完整譜系轉(zhuǎn)化為“界門綱目科屬種株”的自定義格式 |
taxonkit_name2taxid |
將分類單元名稱轉(zhuǎn)化為TaxID |
taxonkit_filter |
按分類學(xué)水平范圍過濾TaxIDs |
taxonkit_lca |
計算最低公共祖先(LCA) |
并且help(taxonkit_*)可查看詳細使用說明。
# 列出[genus] Homo下的所有子單元
taxonkit_list(ids = c(9605), indent = "-", show_name = TRUE, show_rank = TRUE)
## [1] "9605 [genus] Homo"
## [2] "-9606 [species] Homo sapiens"
## [3] "--63221 [subspecies] Homo sapiens neanderthalensis"
## [4] "--741158 [subspecies] Homo sapiens subsp. 'Denisova'"
## [5] "-1425170 [species] Homo heidelbergensis"
## [6] "-2665952 [no rank] environmental samples"
## [7] "--2665953 [species] Homo sapiens environmental sample"
## [8] "-2813598 [no rank] unclassified Homo"
## [9] "--2813599 [species] Homo sp."
## [10] ""
taxonkit_lineage, taxonkit_reformat, taxonkit_name2taxid, taxonkit_filter 與 taxonkit_lca 默認從文件中讀取數(shù)據(jù)纱兑，也可通過指定text = TRUE從字符串輸入讀取輸入數(shù)據(jù)：
# 查詢9606和63221的完整譜系
taxonkit_lineage("9606\n63221", show_name = TRUE, show_rank = TRUE, text = TRUE)%>%
pcutils::strsplit2(split = "\t",colnames = c("taxid","lineage","name","level"))
## taxid
## 1 9606
## 2 63221
## lineage
## 1 cellular organisms;Eukaryota;Opisthokonta;Metazoa;Eumetazoa;Bilateria;Deuterostomia;Chordata;Craniata;Vertebrata;Gnathostomata;Teleostomi;Euteleostomi;Sarcopterygii;Dipnotetrapodomorpha;Tetrapoda;Amniota;Mammalia;Theria;Eutheria;Boreoeutheria;Euarchontoglires;Primates;Haplorrhini;Simiiformes;Catarrhini;Hominoidea;Hominidae;Homininae;Homo;Homo sapiens
## 2 cellular organisms;Eukaryota;Opisthokonta;Metazoa;Eumetazoa;Bilateria;Deuterostomia;Chordata;Craniata;Vertebrata;Gnathostomata;Teleostomi;Euteleostomi;Sarcopterygii;Dipnotetrapodomorpha;Tetrapoda;Amniota;Mammalia;Theria;Eutheria;Boreoeutheria;Euarchontoglires;Primates;Haplorrhini;Simiiformes;Catarrhini;Hominoidea;Hominidae;Homininae;Homo;Homo sapiens;Homo sapiens neanderthalensis
## name level
## 1 Homo sapiens species
## 2 Homo sapiens neanderthalensis subspecies
從文件中讀取數(shù)據(jù)：
names <- system.file("extdata/name.txt", package = "pctax")
taxonkit_name2taxid(names, name_field = 1, sci_name = FALSE, show_rank = FALSE)%>%
pcutils::strsplit2(split = "\t",colnames = c("name","taxid"))
## name taxid
## 1 Homo sapiens 9606
## 2 Akkermansia muciniphila ATCC BAA-835 349741
## 3 Akkermansia muciniphila 239935
## 4 Mouse Intracisternal A-particle 11932
## 5 Wei Shen
## 6 uncultured murine large bowel bacterium BAC 54B 314101
## 7 Croceibacter phage P2559Y 1327037
系統(tǒng)發(fā)育樹
如果是做16S測序的話，在分析過程中就會得到一個帶距離的系統(tǒng)發(fā)育樹化借。宏基因組分析如果組裝MAG后用GTDB-Tk比對數(shù)據(jù)庫后也可以獲得有距離的系統(tǒng)發(fā)育樹潜慎。
但有時候我們想要從物種名或taxid獲取整齊的譜系信息，用來一個構(gòu)建系統(tǒng)發(fā)育樹（層級樹蓖康，沒有真實的距離铐炫，只展示包含關(guān)系）。這是一個常見的需求蒜焊，很多文章都會畫一個這樣的樹圖來展示自己的數(shù)據(jù)倒信。
可以實現(xiàn)這個需求的工具有一些：
- iPhylo：https://iphylo.net/，免費泳梆，快速鳖悠，支持NCBI taxonomy和一些化學(xué)物質(zhì)分類樹，贊
- R包
taxtree优妙，很慢 - PhyloT：https://phylot.biobyte.de/乘综，收費
當然可以使用pctax包快速完成，對于分析流程都在R里做的人來說非常方便:
names <- system.file("extdata/name.txt", package = "pctax")%>%readLines()
# 首先通過`name_or_id2df`獲取整齊的系統(tǒng)發(fā)育分類：
tax_df=name_or_id2df(names,mode = "name")
# 去除部分NA套硼，原因可能是學(xué)名不標準卡辰，或者在新數(shù)據(jù)庫里刪除了，因為taxonomy數(shù)據(jù)庫是不斷變化的
tax_df=na.omit(tax_df)
#用`df2tree`將分類層級表轉(zhuǎn)化為樹對象
tax_tree=pctax::df2tree(tax_df[,3:9])
# tax_tree是phylo對象邪意，可以用ape包直接簡單繪圖
ape:::plot.phylo(tax_tree)
可視化
pctax還提供了一些系統(tǒng)發(fā)育信息展示方法：
- 系統(tǒng)發(fā)育樹
data(otutab, package = "pcutils")
#otutab是豐度數(shù)據(jù)九妈，taxonomy是分類層級表（可通過name_or_id2df獲得）
ann_tree(taxonomy, otutab) -> tree
easy_tree(tree, add_abundance = TRUE) -> p
p
添加主要Phylum的strip：
easy_tree(tree, add_abundance = TRUE,add_tiplab = FALSE) -> p
some_tax <- table(taxonomy$Phylum) %>%
sort(decreasing = TRUE) %>%
head(5) %>%
names()
add_strip(p, some_tax)
當然，更多系統(tǒng)發(fā)育樹的繪制可以參考我之前寫的R繪制優(yōu)美的進化樹（基礎(chǔ)）和R繪制優(yōu)美的進化樹（進階）雾鬼，或者使用iPhylo網(wǎng)站來交互式繪圖：iPhylo 生成并繪制優(yōu)美的分類樹
- 擅戎欤基圖：
sangji_plot(tree)
3.旭日圖
sunburst(tree)
TaxonKit 使用
TaxonKit是采用Go語言編寫的命令行工具， 提供Linux, Windows, macOS操作系統(tǒng)不同架構(gòu)（x86-64/arm64）的靜態(tài)編譯的可執(zhí)行二進制文件呆贿。
發(fā)布的壓縮包不足3Mb嚷兔，除了Github托管外森渐，還提供國內(nèi)鏡像供下載，同時還支持conda和homebrew安裝冒晰。
用戶只需要下載同衣、解壓，開箱即用壶运，無需配置耐齐，僅需下載解壓NCBI Taxonomy數(shù)據(jù)文件解壓到指定目錄即可。
- 源代碼 https://github.com/shenwei356/taxonkit 蒋情，
- 文檔 http://bioinf.shenwei.me/taxonkit （介紹埠况、使用說明、例子棵癣、教程）
選擇系統(tǒng)對應(yīng)的版本下載最新版 https://github.com/shenwei356/taxonkit/releases 辕翰，解壓后添加環(huán)境變量即可使用”芬辏或可選conda安裝
conda install taxonkit -c bioconda -y
表格數(shù)據(jù)處理喜命，推薦使用 csvtk 更高效：
conda install csvtk -c bioconda -y
測試數(shù)據(jù)下載可直接 https://github.com/shenwei356/taxonkit 下載項目壓縮包，或使用git clone下載項目文件夾河劝，其中的example為測試數(shù)據(jù)
git clone https://github.com/shenwei356/taxonkit
TaxonKit為命令行工具壁榕，采用子命令的方式來執(zhí)行不同功能， 大多數(shù)子命令支持標準輸入/輸出赎瞎，便于使用命令行管道進行流水作業(yè)牌里， 輕松整合進分析流程中。
- 輸出：
- 所有命令輸出中包含輸入數(shù)據(jù)內(nèi)容务甥，在此基礎(chǔ)上增加列牡辽。
- 所有命令默認輸出到標準輸出（stdout），可通過重定向（
>）寫入文件缓呛。 - 或通過全局參數(shù)
-o或--out-file指定輸出文件催享，且可自動識別輸出文件后綴（.gz）輸出gzip格式。
- 輸入：
- 除了
list與taxid-changelog之外哟绊，lineage,reformat,name2taxid,filter與lca均可從標準輸入（stdin）讀取輸入數(shù)據(jù)因妙，也可通過位置參數(shù)（positional arguments）輸入，即命令后面不帶 任何flag的參數(shù)票髓，如taxonkit lineage taxids.txt - 輸入格式為單列攀涵，或者制表符分隔的格式，輸入數(shù)據(jù)所在列用
-i或--taxid-field指定洽沟。
- 除了
TaxonKit直接解析NCBI Taxonomy數(shù)據(jù)文件（2秒左右）以故，配置更容易，也便于更新數(shù)據(jù)裆操，占用內(nèi)存在500Mb-1.5G左右怒详。 數(shù)據(jù)下載：
# 有時下載失敗炉媒，可多試幾次；或嘗試瀏覽器下載此鏈接
wget -c https://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz
tar -zxvf taxdump.tar.gz
# 解壓文件存于家目錄中.taxonkit/昆烁，程序默認數(shù)據(jù)庫默認目錄
mkdir -p $HOME/.taxonkit
cp names.dmp nodes.dmp delnodes.dmp merged.dmp $HOME/.taxonkit
Taxonkit的作者大大貼心地提供了中文文檔：https://bioinf.shenwei.me/taxonkit/chinese/吊骤，非常詳細，大家可以參考使用静尼。
