BWA マッピングツール

BWA，Burrows-Wheeler Alignment Tool　を、ひとから勧められたので試してみた。
このアライメント（マッピング）ツールは、ギャップを許してアライメントすることができる。
また、リード配列は全部が同じ長さである必要が無い。
ということで、SNPなどのバリアントを見つけるにはいいと、そのひとは言っていたのだが、ここは議論の余地が残るところだ。

さて、そんなことはさておき、さっそく実行ファイルをダウンロード
http://sourceforge.net/projects/bio-bwa/
http://sourceforge.net/projects/bio-bwa/files/

インストールするマシンはLinuxの64ビット。
ダウンロードしたら、以下のコマンドを投げていこう！
$ bunzip2 bwa-0.5.9rc1.tar.bz2　これは解凍ね
$ tar xvf bwa-0.5.9rc1.tar
$ cd bwa-0.5.9rc1/　ディレクトリに移動して、
$ make　メークして、
$ sudo apt-get install bwa　スーパーユーザー権限でインストール！

だ、だ、だ、、、、ってインストールが完了したら、準備OK

失敗したら、近くのLinuxマニアに助けを呼ぼう。

ここからが本題。

BWAはBurrows-Wheeler Transform (BWT)を基本としているアライメントツールで、ショートリード向けのものと、ロングリード向けのものの2種類がある。
ショートリード向けは200塩基未満のリード長を3%未満のエラー率で、ギャップを許すグローバルアライメントを行い、ペアエンドにも対応している。
ロングリード向けのものはBWA-SWといい、Smith-Watermanアルゴリズムに似たヒューリステッィクなアルゴリズムで高いスコアを見つけながらローカルアライメントを行う。　

ショートリードでもこのBWA-SWを実行できなくはないが、エラーは多くなる。

・・・なんて、いかにも前から知っていたような書き方をしたが、ここに書いてあるのを訳しただけ。

もっと詳しく知りたいひとは文献までさかのぼってみよう。

さて、ランの流れだ。　こんな感じで進む。

1．データを用意する
リファレンス配列はFASTAフォーマット
リード配列はFASTQフォーマット　今回は例としてイルミナ・ペアエンドのショートリードを使う。

ほかのフォーマットに対応しているかは？

2．リファレンスデータベースの作成
リファレンス配列のデータベースを先に作る必要があり、FASTAファイルに、indexコマンドを使ってインデックスを貼る。

3．アライメントコマンド実行
アライメントはalnコマンドで実行する。Suffix Arrayという検索アルゴリズムの原理を使ってひとつひとつのリードのベストヒットを検索するらしい。

4．SAMファイルへ出力
alnコマンドで作成したsaiファイルを基に、マッピングファイルSAMを作成する。
samse/sampeコマンドはSuffix Arrayのインデックスを染色体の位置情報に変換する。
シングルリードならsamse、ペアエンドリードなら、sampeを実行。

sampeだよ。sampleと打ってしまいそうなので注意！

それでは実例を

まず、リファレンスに対し、インデックスを貼る（ヒト染色体1番のみを例に）
$ bwa index -p ./work/humanCh1 -a is ./Path_to_directory/hs_ref_GRCh37_chr1.fa

重要パラメータ
  -p: アウトプットするリファレンスデータベースの名前
  -a: インデックスの種類を指定　-a is または -a bwtsw
isは簡単で速い。ただしリファレンス塩基数の5.37倍のメモリを必要とし、2GB以上のサイズのデータベースは作れない。
bwtswは、BWT-SWのため。2GBの制限が無いので、ヒトの全ゲノムでデータベースを作ることができる。ただし10MB未満のデータベースは作れないし、速度もISに比べて遅い。

こんなインデックスファイルができたかな？



リファレンスデータセットの名前をhumanCh1にしたので、以下、alnの後のデータベース名はhumanCh1 になる。

リードファイルは、SRR027863_1.fastqとSRR027863_2.fastqで、これはペアエンドのファイルだ。

ひとつずつ、アラインさせる。

$ bwa aln ./work/humanCh1 ./work/SRR027863_1.fastq > ./work/SRR027863_1.sai
$ bwa aln ./work/humanCh1 ./work/SRR027863_2.fastq > ./work/SRR027863_2.sai

アラインして、拡張子saiのファイルができたら、いよいよマッピングファイルSAMの作成。
これはペアエンドなので、sampe

$ bwa sampe ./work/humanCh1 ./work/SRR027863_1.sai ./work/SRR027863_2.sai ./work/SRR027863_1.fastq ./work/SRR027863_2.fastq > ./work/SRR027863_Ch1.sam

ちゃんとパスを通しておけば、もっとすっきりするけどね。
基本は、
$ bwa sampe 「インデックスデータベース」 「アラインファイル1.sai」 「アラインファイル2.sai」 「リードファイル1.fastq」 「リードファイル2.fastq」 > 「マッピングファイル.sam」

これでSAMファイルは出来上がる。1.26Gくらいのでかいファイル。

今回は、シンプルにデフォルトで行った。

全部通しでも1時間くらいで終わった。

気をつけるとしたら、ショートリード用のリファレンスインデックスをつくるには、2GBのサイズ制限があるので、そのままヒト全ゲノムは使えないということ。

染色体ごとに分けるか、目的に応じて、いらない配列を削るかする必要があるかな。