2020 9/9,9/10 コード修正
タイトルの通り。
インストール
condaでpython3.7の仮想環境を作ってテストした（macos10.14）。
依存
- ncbitax2lin requires python-3.7
conda create -n ncbitax2lin -y python=3.7
conda activate ncbitax2lin
pip install -U ncbitax2lin
> ncbitax2lin
NAME
ncbitax2lin - Converts NCBI taxomony dump into lineages.
SYNOPSIS
ncbitax2lin NODES_FILE NAMES_FILE <flags>
DESCRIPTION
NCBI taxonomy dump can be downloaded from
ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz
POSITIONAL ARGUMENTS
NODES_FILE
path/to/taxdump/nodes.dmp from NCBI taxonomy
NAMES_FILE
path/to/taxdump/names.dmp from NCBI taxonomy
FLAGS
--output=OUTPUT
NOTES
You can also use flags syntax for POSITIONAL ARGUMENTS
データベースの準備
NCBI taxonomyからtaxdumpファイルをダウンロードする。
ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/
wget -N ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz
mkdir -p taxdump && tar zxf taxdump.tar.gz -C ./taxdump
実行方法
ncbitax2lin taxdump/nodes.dmp taxdump/names.dmp
元のファイル
> head names.dmp
> head nodes.dmp
出力
先頭だけexcelで開いた。
補足
localでnucleotide collection(Nr/Nt)データベースに対してblastサーチを行い(実行例)、得られたNCBIのGenbank accession numbersからlineageファイルに変換したい場合、まずNCBI FTPサーバにあるnucl_*accession2taxid.gzをルックアップテーブルとして使いaccession IDからtaxIDに変換する（ちなみにGI numberはすでに廃止されている） 。
nucl_*accession2taxid.gzをダウンロードする。
ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/
それから統合してデータベースにする。
#accession2taxidを結合、”sort -u”で重複を除きつつソート
cat nucl_gb.accession2taxid nucl_wgs.accession2taxid dead_nucl.accession2taxid dead_wgs.accession2taxid |sort -u > concatenated_accession2taxid
#1列目とGI numberは不要（２列目はバージョンつき。".1"ならバージョン１）
cut -f 2,3 concatenated_accession2taxid > genbank_AC_db
#高速なripgrepでaccesion IDを検索（*1）（無いならbrew install ripgrep）、taxidを取り出す。"-m 1"で１ヒットのみ取り出す。
#例えばID AAAA02000036を結合したデータベースgenbank_AC_dbから探す
rg -m 1 AAAA02000036 genbank_AC_db |cut -f 2 > taxid_list
それから、NCBItax2linで作成したリストファイルを参照するか、taxonkitなどを使ってlineageファイルに変換する。taxonkit を使うなら
taxonkit lineage taxid_list > output_lineage
#ワンライナーで
rg -m 1 <genbank_accesion_ID> genbank_AC_db |cut -f 2 | taxonkit lineage > output_lineage
#複数IDをループ処理.echoで１列目にクエリのID、２列目に見つかったIDを出力
cat genbank_accesion_ID_list | while read line
do
echo -n -e "$line\t" >> taxid_list ;rg -m 1 $line genbank_AC_db |cut -f 2 >> taxid_list
done
#最後にtaxid => lineage
taxonkit lineage taxid_list >> lineage_list
出力
lineage_list
（taxid1つを処理）
引用
GitHub - zyxue/ncbitax2lin: 🐞 Convert NCBI taxonomy dump into lineages
Used in
Mahmoudabadi, G., & Phillips, R. (2018). A comprehensive and quantitative exploration of thousands of viral genomes. ELife, 7. https://doi.org/10.7554/eLife.31955
参考
RefSeq と GenBank の違い
Question: Accession number to taxonomy id after blasting
https://www.biostars.org/p/222183/
ripgrep
*1
index検索してないので低速だが、急ぎでないなら問題なし。macだとgrepがとんでもなく遅い場合があるのでlinux 推奨。
関連
