fastacmdの後継者: blastdbcmd
遺伝子上流配列など、ゲノム配列の任意の場所を切り出すのには、BLATパッケージに含まれている nibFrag をこれまで勧めてきました。FASTAフォーマットの配列ファイルに対して別個にnibインデックスを作らないといけなくて管理が面倒なのと、やはりBLATはfor profitには有償という大きな難点がありました(ライセンス料は決してお安くない)。
そこでpublic domainなソフトウェアのNCBI BLASTパッケージにも似たプログラムがきっとあるはずだろうということで探してみると… fastacmd コマンドがそれだったようです。しかしながらBLASTパッケージのバージョンアップでBLAST+になってからは fastacmd が含まれなくなっているようです。更に調べると blastdbcmd というプログラムがその後継という位置づけのようで、BLAST+をインストールする時に blastdbcmd は同時にインストールされるプログラムです。その使い方を調べてみました。
まず、makeblastdb コマンドでBLAST用のindexを作成します(かつての formatdb コマンド相当だと思います)。-dbtype nucl で核酸配列であることを指定し、-hash_index でインデックスの作成も指定しておきます(追記: -parse_seqids を足さないと個別エントリの切り出しには対応しない模様。別のblastdbcmd関係エントリ参照。ここで -in で指定している 1.fa はゲノム配列(1番染色体)が1本だけ入ったFASTAフォーマットのファイルです。
#!/bin/sh
makeblastdb -dbtype nucl -hash_index -in 1.fa
そして blastdbcmd の指定の領域を -range オプションで指定すれば良いようです。この例の場合は5000000番目から5000100番目ということになります。
#!/bin/sh
blastdbcmd -db 1.fa -entry all -range 5000000-5000100
そうすると以下の様な結果が返ってきて取得できました!
>gnl|BL_ORD_ID|0:5000000-5000100 1 dna:chromosome chromosome:NCBIM37:1:1:197195432:1
GGTGGGCCTTGGTGATGCAGCTGTCTTTGAGTCATTTCTGCTCCTGTAAGCAACCTCTCACCCATATTCTTATAAGTAAC
CCTAGTAAAACTCATGGTTAA
追記：より新しい「localBLASTのblastdbcmdを使って配列エントリの特定の領域を取得する」というエントリも参照。