自家製BLAST用DBから必要な配列エントリ取得

BLASTでヒットがあったエントリの配列データを切り出してきたい時はままあるかと。公共データベースのエントリならtogowsのRESTで取得するなり、ググるなりNCBIで検索するなりで配列データにありつけるかと思います。ところが、自分で作成したDBに対してBLASTした場合にはそうもいかないことがあるわけで。その場合どうするか、なのですが、とくに新規のコマンドをインストールしなくてもできるいい方法があります。

それにはまず、BLASTのDBを作成するコマンドmakeblastdbを実行するときに-parse_seqidsというオプションも足して実行します。例えば、自分の手持ちのDBがoreno.aaというFASTA形式のファイルだとすると

#!/bin/sh

makeblastdb -in oreno.aa -dbtype prot -hash_index -parse_seqids

な感じで。 そして、BLASTヒットがorenogene1だったとすると

#!/bin/sh

blastdbcmd -db oreno.aa -entry orenogene1

とすればFASTA形式でそのエントリの配列が取得できます。ただ、FASTAのヘッダ行のラベルが|(パイプ)などで長くなっていると、行頭からスペースが現れるまでの文字列を指定する必要があるようです。例えば、 >gi|5524211|gb|AAD44166.1| cytochrome b のようになっている場合です。その場合は以下のように指定することになります。

#!/bin/sh

blastdbcmd -db oreno.aa -entry 'gi|5524211|gb|AAD44166.1|'

Written by bonohu in misc on 金 08 8月 2014.