bwaとは
bwaは高速シーケンサーによって出力される大量の塩基配列(リード)を参照配列(リファレンス)のどこと一致するのかをすばやく検索するプログラムです。通常こうしたソフトウェアをマッパー(mapper)と呼びます。
bwaはBurrows-Wheeler Alignment Toolを意味します。BurrowsとWheelerは1994年にブロックソートと呼ばれるデータ圧縮の前処理に使われる文字列の可逆変換のアルゴリズムを開発しました。bwaはこのアルゴリズムを使って参照配列の索引をつくって高速な検索を可能にしていると思われます。
bwa用のインデックスの作成
$ bwa index reference.fasta
これによりカレントディレクトリにあるreference.fastaからreference.amb、reference.ann、reference.bwt、reference.pac、reference.saの5つのファイルが作成されます。およそ360MBの塩基配列からインデックスを作成するのに5分ほどかかりました。
FASTAファイルとは別のディレクトリにインデックスを作成するには以下のようにします。
$ bwa index dir/reference.fasta -p dir/reference_name
とりあえずやってみる
デフォルトのパラメーターで動かしてみます。
$ bwa mem reference fastq.gz |grep -v '^@' >out.sam $ head -n 10 out.sam 55KC:06608:08873 16 ch01 1044 0 55S52M * 0 0 TTCAGCCGAGCTATAGCTCATGAAATTTTGGCTGATAGCNAAAAACCATGAGTTATAAACCCTAAACCCTAACCCTAAACCCTAACCCTAAACCCTAAACCCTAAAC 779;:38,,,,,888<<<<?==5<:4<<<6@<<<==<<;*;;;;;4;;<;<:67785;;1;87+;<5<<=8<7<<<8<=8>==9=8===8@=7<=<4?:4:46.::4 I55KC:07929:13415 16 ch01 1044 1 112S67M * 0 0 ATTTGGGCCAGGTNAAAAAATCACGTGTTTTTCAAAGNTTAACACTTAGAGGCAGTTTCAGCCGANGCTATAGCTCATGAAATTTTGGCTGATAGCAAAAAACCATGAGTTATAAACCCTAAACCCTAACCCTAAACCCTAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACC <7>8.781::6<8'8888888??<<<;1;;;;;1;;8*880==<<8===<8<<==8=?=<:39,,)---668;:9:=C<5>>2<<<9=@>?@@=<<,??@?<7>=<<<=9<=D6C909997>>8>>>9>8?;;6<;5<=;6;6<=<8><8==<6<=8===7=<6<<<8?=8=<<6<<4? I55KC:06385:03280 16 ch01 1056 0 5S57M * 0 0 CTNANCCTAACCCTAAACCCTAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCT //)/)///392::;4<<5><;5;4=;;3<;6;<;6?:3:;:6<<6=<<7<<7><<6><6<<< I55KC:04516:06048 16 ch01 1079 60 26S84M * 0 0 ACCTNACCGTNANCTAACCTAANCCTAANCCTNAANCCTAAACCCTAAACCCTAAACCCTAAACAGCTGACAGTACGATAGATCCACGCGAGAGGAACCGGAGAGACAAC 62--(42---)-*--59495/5*77708*;;8*88*=<=5=;6<<<7=<7@><6<;7@@A8@====<=====<>==<<<>===9==?BCCE?>>?9===9=>====>9=9 I55KC:01768:02245 16 ch01 1100 60 18M1D166M * 0 0 AACCCTAAACCCTAAACACTGACAGTACGATAGATCCACGCGAGAGGAACCGGAGAGACAACGGGATCCAGGCGCCAGCGACGGATCCGGGCGAGAGGGGAGTGGAGATCATGGATCCGTGCGGGAGGGGAAGAAGTCGCCGAATCCGACCCTCCCATCGCCATCGACAGTAGGTCTGCGCGAG 3;1899.883;..,?...9:99:89678<<<<<<<8<<=<<<;;;7;6;661577<;77299499::6:<8>>=7><<>8886///7>6=>=<<<<2<<<;;;5;;<<<<=<8<<=9;;;><6<<<4@@83888<<<;>9A=9=>9A<909986<<<<<<7=<=>=<==@>?=>>========= I55KC:02103:07442 16 ch01 1115 60 82M1D91M * 0 0 ACAGCTGACAGTACGATAGATCCACGCGAGAGGAACCGGAGAGACAACGGGATCCAGGCGCCAGCGACGGATCCGGGCGAGAGGGAGTGGAGATCATGGATCCGTGCNGGANGGGAAGAAGTCGCCGAATCCGACCCTCCCATCGCCATCGACAGTAGGTCTGCGCGAGAGGG >>==>>757<<<<=>>==<<<8<<><<<===9>9?9>9==?===<8<<8===>9=A9>==9?@=====9=;;7=7==<;::7+888:;4<<=<<<=<8;;;6;;;;8*775*8880;;8=<><=;DD9=<8<<<6<<<8>?>=>>9>=>@>=<?>@<9==<;::<=<<<;4;< I55KC:09422:11780 16 ch01 1128 60 76M1I23M2I7M1I20M1I4M * 0 0 CGATAGATCCACGCGAGAGGAACCGGAGAGACAACGGGATCCAGGCGCCAGCGACGGATCCNGGCGAGANGGGAGTGGGAGATCATGGATCCGTGCGGGANGGGGGAAGNAAGTCGCCGAATCCGANCCTNCCCA <<<<<::84///:::<=<8<8<8<8::;<:994982:7755:84::;8<;55595/55835*557===5*9998<9+999;;;<<<8<<;599995+877'::::;6;5*588<<<8<<8995884)444(8888 I55KC:05782:12911 0 ch01 1129 60 88M1D65M1I3M1D40M * 0 0 GATAGATCCACGCGAGAGGAACCGGAGAGACAACGGGATCCAGGCGCCAGCGACGGATCCGGGCGAGAGGGGAGTGGAGATCATGGATCGTGCGGNAGGGNAAGAAGTCGCCGAATCCGACCCTCCNATCGCCATCGACAGTAGGTCTGCGCGNAGAGGNCACCGGCGCTGGCTCTGAGCGAGATGCAACGCCGGCC 8=<===<<9=<<<<<<;;6;7<9@=>>=>=><9=<:475548;6;9:6:;;<:<;7;;<5;;49<<=>?>>7<<<<9<=<==<::5:66399666(3333'3-3305534427:06647:77.133*1778848<::;;::7999=::9::--)-4---(-77,60677384888:576:<=:<9>==8===8>8?7
277万のリードを360Mbpの参照配列にあてるのに約23分かかりました。リード数は2,770,697であるのに対して出力は2,804,375行ありました。
bwaの出力
bwaの出力はSAM形式のテキストファイルとなっています。SAM形式の各列の意味は次の通りです。
| 1列目 | リードのID |
|---|---|
| 2列目 | FLAG: リードのマップ状況を示す数値。2進数にして各桁ごとに意味を読みとる。 |
| 3列目 | 参照配列の名前 |
| 4列目 | リードと一致した参照配列の位置。1から始まることに注意。 |
| 5列目 | MAPQ: リードと参照配列の関係が1:1に近いと高い数字になる(0-60)。 |
| 6列目 | CIGARというアライメントを表す文字列。 |
| 7列目 | ペアリードが与えられたとき相手方のリードの名前 |
| 8列目 | ペアのリードがあたった位置 |
| 9列目 | ペアであたったときリファレンスから推定されるインサートの長さ。 |
| 10列目 | リード本体 |
| 11列目 | リードの塩基配列の質 |
| 12列目 | その他あれば |
FLAGの数値は2進数に変換して各桁ごとに意味を読みとります。
| 0x0001 | このリードはペアの方割れである |
| 0x0002 | ペアが適切にマップされた |
| 0x0004 | このリードはマップされなかった |
| 0x0008 | ペアの相手はマップされなかった |
| 0x0010 | 参照配列とは逆向きである |
| 0x0020 | ペアの相手は参照配列と逆向きである |
| 0x0040 | ペアの最初のほうである |
| 0x0080 | ペアの後のほうである |
| 0x0100 | このアライメントは第一候補ではない |
| 0x0200 | このリードは品質検査に不合格であった |
| 0x0400 | 同じ塩基配列が存在する(PCRで過剰に増やしすぎたか、画像処理で同じ点を二回読んでいる可能性がある) |
| 0x0800 | このアライメントは予備である(samtoolsのマニュアルにはのっていない)。参照 |
FLAGの数値解釈ツール
FLAGの数値を入力してください: (上記の表の一致する項目がハイライトされます)
bwaのマッピング結果の解釈
2,770,697リードをゲノムにあてたところ2,804,375行の出力がありました。
染色体ごとのリード数を出したものが以下の表です。
| 3列目 | 出現回数 |
|---|---|
| ch01 | 2605354 |
| ch02 | 18406 |
| ch03 | 11840 |
| ch04 | 15623 |
| ch05 | 8430 |
| ch06 | 9272 |
| ch07 | 6587 |
| ch08 | 6519 |
| ch09 | 8404 |
| ch10 | 10341 |
| ch11 | 9032 |
| ch12 | 15424 |
| * | 79143 |
1番染色体にほとんどがあたっているのは使用したfastqファイル(Torrent Serverよりダウンロード)にもともとある偏りによるものです。
3列目が*になっているものはマップされなかったものを示します。これらは除く必要があるでしょう。
次にMAPQ値の分布を調べてみました。
60%強のリードがマッピングスコア60で参照配列にあたっていました。 MAPQ値が高いとそのリードはゲノムのあたった領域を読んだものである可能性が高いと考えられます。 一方MAPQ値が0のものが約20%含まれていました。これらはゲノムの様々な場所に一致し、どこを読んだものであるのか判別できないことを示します。 試しにMAPQ値が0であったリードをbwaの-aオプションをつけて再度検索を行うと48箇所もゲノム上の場所がでてきました。 MAPQ値が低いものはSNPの検索には使わない方がよいでしょう。 bwaの-Tオプションを使うと指定した数値以下のMAPQ値のものを除くことができます(実際にはFLAGが4になる)。
重複してあたっているリードを調べる
各リードあたり何回出力されているのか調べたのが下の表です。
| 1回 | 2,738,176 |
| 2回 | 31,365 |
| 3回 | 1155 |
| 4回 | 1 |
| 計 | 2,770,697 |
ほとんどのリード(98.8%)が1回だけ出力されています。 同じリードが複数個所に当たっているときの状況を抜き出してみました。
I55KC:10195:05899 0 ch12 13408042 60 142S55M * 0 0 I55KC:10195:05899 2048 ch09 9343681 48 106H53M38H * 0 0 I55KC:10195:05899 2048 ch01 24842224 36 35M162H * 0 0
Read AGCCNAANGGGACCGGAACAGTCTAGTCTCTAGCTGCCTAA(...)CCCTACTGTCCTAATTCTGCTTTAGGACAAAGGAAAGGGNAAGACTTCCTCTAATAAAGGAGGAAAGGAGAGATTTTCATTAAACCAAAAGCCAACCA ch012 ACTTCCTCTAATAAAGGAGGAAAGGAGAGATTTTCATTAAACCAAAAGCCAACCA ch09 TCCTAATTCTTCTTTAGGACAAAGGAAAGGGGAAGACTACCTCTAATAAAGG ch01 TGCCCAAAGGGACCGGAACAGTCTAGTCTCTAGCT
先頭のアライメントが最もMAPQ値が高く55bpにわたって完全一致でした。 リードは全体で197bpの長さがあり、残りの142bpは由来が不明です。 次のアライメントはMAPQ値が48で、52塩基の一致と1塩基の不一致を含みます。 最後のアライメントは35bpの範囲にわたっていますが、リードに二つのNが含まれる領域で1塩基のミスマッチが存在します。
複数のアライメントが出力されたとき二つ目以降のアライメントのFLAGには0x0800がつくようです。 おそらくミスマッチやギャップの数を考慮してより一致度の高いものを優先していると考えられます。 次はMAPQがともに60である二つのアライメントの例です。
I55KC:04859:11456 0 ch11 26799599 60 54S138M * 0 0 I55KC:04859:11456 2048 ch01 10652610 60 54M138H * 0 0
Read CATTGCCATCCTAAGTTGAAATGGTGTACATAAAGCCCTCCACCCAAATACTGTCAAAACTGGGGTCATTCGCCTGTCACTGCAGCGGAGCAAAACTGGGGTCATTCACCTGTCAGTGAAGCAGGTAGCATGAACTACAATGGGGCTGTCAATGAAGCAGGCAGTTGGAGCCACAGGGGACTTCCTCCTTCA ch11 CAAAACTGGGGTCATTCGCCTGTCACTGCAGCGGAGCAAAACTGGGGTCATTCACCTGTCAGTGAAGCAGGTAGCATGAACTACAATGGGGCTGTCAATGAAGCAGGCAGTTGGAGCCACAGGGGACTTCCTCCTTCA ch01 CATTGCCATCCTAAGTTGAAATGGTGTACATAAAGCCCTCCACCCAAATACTGT
FLAGの0x0800を見るとより一致度の高いアライメントを採用できそうです。 おそらくこの例はライブラリの作成過程でキメラになったものと考えられます。
bwaをうまく使う方法
- -T オプションに適切な数値を与えて不正確なアライメントを除く必要があると考えられます。
- リファレンスに一致しなかったもの(3列目が*またはFLAGの0x0004が1のもの)を除去するとデータ量が減らせるでしょう。
- CIGARを見てリードをトリミングする必要があるでしょう。
- 挿入と欠失のあるリードをどう扱うべきか検討が必要です。欠失(リードのほうが短い)ときは-を入れて、挿入のときは前の塩基を数字に置き換えるというのはどうかな。