Local BLASTの結果からLCA (Lowest common ancestor)で由来生物種を推定する

きっかけ

最近仕事でLocal BLASTの結果をパースして配列が由来する生物種を推定する機会があったのだが、なかなか思うような機能を実装したツールが存在していなかった*1ので作ってみた。

準備

以下のインストール・ダウンロードを行う

• NCBI BLAST+
• BLAST database files
  • 以下のDBはBLAST用のフォーマットに変換されたものがNCBIによって提供されているのでそれをダウンロード
    • RefSeq protein
    • NCBI nr
• Python 3
  • sqlite

使い方

スクリプト

https://github.com/tmaruy/blast2taxonomy

1. SQLiteでDBを作成する

git clone https://github.com/tmaruy/blast2taxonomy.git
cd blast2taxonomy
chmod +x preprare_taxdb.sh
./prepare_taxdb.sh

上記のコマンドで次の2つのテーブルを含むDBファイル (taxonomy.db) が作られる。トータルで30GB強のサイズとなるので注意。

• accession2taxid: NCBIのアクセッション番号とTaxonomy IDの対応表
  • ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/から作成
• taxonomy: Taxonomy IDと生物種名の対応表
  • ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdmp.zipから作成

accession2taxid

SELECT * FROM accession2taxid LIMIT 5;

accession|taxid
A0A009IHW8|470    
A0A023GPI8|232300 
A0A023GPJ0|716541 
A0A023GS28|37565  
A0A023GS29|37565  

taxonomy

SELECT * FROM taxonomy LIMIT 5;

taxid|category|name|t_domain|t_phylum|t_class|t_order|t_family|t_genus|t_species
1|no rank|root|||||||
2|superkingdom|Bacteria|Bacteria||||||
6|genus|Azorhizobium|Bacteria|Proteobacteria|Alphaproteobacteria|Rhizobiales|Xanthobacteraceae|Azorhizobium|
7|species|Azorhizobium caulinodans|Bacteria|Proteobacteria|Alphaproteobacteria|Rhizobiales|Xanthobacteraceae|Azorhizobium|Azorhizobium caulinodans
9|species|Buchnera aphidicola|Bacteria|Proteobacteria|Gammaproteobacteria|Enterobacterales|Erwiniaceae|Buchnera|Buchnera aphidicola

2. BLAST+

タブ形式で結果を出力 (-outfmt 6)

blastx -db /path/to/refseq_protein \
       -query test/input.fasta \
       -outfmt 6 \
       -out test/blast.txt

3. 由来生物種推定

• -i : Input file, result of blast search (made at step 2)
• -db : DB file taxonomy.db (made at step 1)
• -o : Output file
• -evalue : Threshold of evalue (Ignore hits if their evalues are above this threshold)
• -lca : Top N hits used for calculation of lowest common ancestors (default N=10)
• -lca_threshold : Returns classification if shared by X % of top hits (default X=50)

#git clone https://github.com/tmaruy/blast2taxonomy.git
python blast2taxonomy_lca.py -i test/blast.txt \
                             -db taxonomy.db \
                             -o test/tax_lca.txt \
                             -evalue 0.01 \
                             -lca 10 \
                             -lca_threshold 50

*1:すでにサポートが終了したものは見つかった https://github.com/frederikseersholm/blast_getLCA