2010年12月28日火曜日

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のサイズ制限があるので、そのままヒト全ゲノムは使えないということ。
染色体ごとに分けるか、目的に応じて、いらない配列を削るかする必要があるかな。






0 件のコメント:

コメントを投稿