localBLASTのblastdbcmdを使って配列エントリの特定の領域を取得する

新規なモデル生物では、ゲノム配列が公開されていても多くの場合染色体ごとに配列が一本につながっていない。どうなっているかというと、scaffoldという単位でゲノム断片配列が記述されている。 そういった場合でも、コマンドラインで使用する際のNCBI BLASTに含まれる makeblastdb というプログラムを使ってインデックスを作成、blastdbcmd を使ってエントリ名(-entry)と領域(-range)を指定することで、その領域の塩基配列を取得することが可能である。

#!/bin/sh

makeblastdb -in hoge.fa -dbtype nucl -hash_index -parse_seqids
blastdbcmd -db hoge.fa -entry scaffold001 -range 2000-2500

このコマンドでは、一行目でhoge.faファイル(FASTA形式)に対してインデックスを作成、二行目でscaffold001の2000塩基から2500塩基までの配列を抽出する例を示している。通常の makeblastdb のオプションに加えて、-parse_seqids が加わっている点に注意。


Written by bonohu in misc on 木 02 6月 2016.