はじめに（おことわり）

この記事はかなり雑な記事です。しかし、bwaのインデックスについて記述している記事はあまり見ず、また公式の仕様書もないので、間違っているかもしれないけれども一応こういう情報もあるよという程度の認識で読んでいただけると幸です。

bwaとは、バイオインフォマティックスの領域で広く利用されているアライナーです。現在でもショートリードのDNAでは最もよく利用されているツールです。（後続の bwa-mem2というソフトもあります。）

bwaでは、リファレンスゲノムのインデックスを作成します。

bwa index reference.fa

インデックスを作成すると、.amb, .ann, .pac, .bwt, .sa の5つのファイルが作成されます。
これらのファイルの中身って一体何なんでしょう？

とういことを調べましたが .bwt, .sa は難しすぎて完全にはわかりませんでした。

この2つのファイルはアルゴリズムと密接な関連があり、FM-index,BWTといった高度なアルゴリズムは、私の知力ではとても理解が追いつきません。

しかし、だからといって、全く何も得るものがなかったわけではありません。それ以外の .amb, .ann, .pac に関してはかなり理解できた部分もあったので、理解できた部分に限って記録に残します。

5つつのファイルのおおまかな意味

こちらのサイトで、5つのファイルの大まかな意味が解説されています。それによると、

拡張子	種類	説明
.amb	テキスト	リファレンス fasta 内の N (またはその他の非 ATGC) の出現を記録します。
.ann	テキスト	リファレンスシーケンス、名前、長さなどを記録します。
.bwt	バイナリ	Burrows-Wheeler 変換されたシーケンスです。
.pac	バイナリ	パッケージ化されたシーケンスです (4 つの塩基対が 1 バイトをエンコードします)。
.sa	バイナリ	サフィックス配列インデックス

どうでしょう。これだけでもかなりのことが理解できました。

バイナリファイルを閲覧する方法 hexyl を利用する

さて、bwaのインデックス関連ファイルのうち、ambとannをのぞく3つのファイルはバイナリファイルです。これらのファイルについては今回十分に調査できたわけではないのですが、Rustでかかれた hexyl というコマンドを利用すると比較的簡単にバイナリファイルを閲覧することができます。おすすめです。

.amb は ambiguous を意味する

このambファイルは、リファレンスのACTG以外の塩基対が含まれているときに記録するものです。ACTG以外に含まれている塩基としては当然Nがあります。つい最近、ヒトゲノムのY染色体をのぞくほとんどすべての配列が解明されたという発表がありました。しかし、GRCh38などには多数のNが含まれています。

たとえば zcat GRCh38.p13.genome.fa.gz | head などとすると

>chr1 1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

こんな感じで初っ端からNがてんこ盛りです。人のリファレンスゲノムには通常、ACTGおよびN以外の文字は含まれていません。このことは、下記のような簡単なコマンドでもすぐわかります。

zcat GRCh38.p13.genome.fa.gz | grep -v '^>' | sed -e '/[ACTGN]/d' -e '/^$/d'

なんなら文字の数をカウントすることもできます。こちらは実行時間がそこそこかかるのでcrystal言語で書きました。

count.cr

count = {'A' => 0, 'C' => 0, 'T' => 0, 'G' => 0, 'N' => 0}
i = 0
while l = gets
  l.chomp.each_char {|c| count[c] += 1}
  i += 1
  p count if i % 1000000 == 0
end
puts "---"
p count

zcat GRCh38.p13.genome.fa.gz | grep -v '^>' | ./count

こんな感じで実行すると、

---
{'A' => 914265135, 'C' => 635937481, 'T' => 916993511, 'G' => 638590158, 'N' => 161331703}

こんな結果が得られます。contigのこともありますが、Nだけでみても5%程度は配列が不明な部分があることがわかります。幸いにして、ヒトのリファレンスゲノムにはA, C, G, T, N以外の不明な配列はありません。しかし、これは幸運なことで、実際にはACTGN以外の配列が入り込むことがあります。WikipediaのNucleic acid sequenceをご覧になると、ACTG意外にも非常に多くの記号が入り込む可能性があることがわかります。

これらの記号はめったに見ることはありませんが、ライブラリを作成するような場合には注意が必要です。

実際にためしてみる。

ここで、適当なリファレンスファイルをでっちあげてみます。poo.fa はA64文字からのみなるリファレンス配列です。

poo.fa

>chr1 1
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA

このリファレンスのインデックスを作成してみます。poo.fa.amb は

poo.fa.amb

64 1 0

のようになります。ちなみにリファレンス配列の染色体を増やすと、

poo.fa

>chr1 1
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
>chr2 2
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC

このようになります。

poo.fa.amb

128 2 0

それでは、適当にNを混ぜ込んでみましょう。

poo.fa

>chr1 1
NAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA

こうなります。

poo.fa.amb

64 1 1
0 1 N

こうすると、

poo.fa

>chr1 1
NNAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
>chr2 2
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCN

こうなります。

128 2 2
0 2 N
127 1 N

これで、なんとなくルールがわかったと思います。最初の1列目は {塩基の数、染色体の数、N等の数｝2列目以降は {Nの開始位置、終了位置、N等}となっているのだと思います。違うかもしれません。違う場合は教えてください。

.ann

こちらはもう少し単純なファイルの様子で、

>chr1 1
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
>chr2 2
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC
CCCCCCCC

のとき

みたいになります。完全にはわかりませんが、配列名や、塩基長、オフセットなどを記録していると思われます。

.pac は配列をバイナリでパックしている

さて、pacファイルですが、実は、bwa index でつくられる pac ファイルと、 bwa fa2pac コマンドでつくられるファイルは少し違います。

poo.fa

>chr1 1
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA

があったときに、bwa index poo.fa で得られるpacファイルは heyxl poo.fa.pac とすると

┌────────┬─────────────────────────┬─────────────────────────┬────────┬────────┐
│00000000│ 00 00 00 00 00 00 00 00 ┊ 00 00 00 00 00 00 00 00 │00000000┊00000000│
│00000010│ 00 00                   ┊                         │00      ┊        │
└────────┴─────────────────────────┴─────────────────────────┴────────┴────────┘

こんな感じです。一方で、hexyl poo.fa.pac する場合は、

┌────────┬─────────────────────────┬─────────────────────────┬────────┬────────┐
│00000000│ 00 00 00 00 00 00 00 00 ┊ 00 00 00 00 00 00 00 00 │00000000┊00000000│
│00000010│ ff ff ff ff ff ff ff ff ┊ ff ff ff ff ff ff ff ff │××××××××┊××××××××│
│00000020│ 00 00                   ┊                         │00      ┊        │
└────────┴─────────────────────────┴─────────────────────────┴────────┴────────┘

こんなふうになります。この違いはなんでしょう？その秘密は -f オプションにあって、どうやら、pacには配列全体のリバースコンプリメントを詰め込んでいるようなのです。考えてみると、bwaなどのアライナーは器用に相補的な塩基配列に対してもマッピングしてくれます。その秘密はこんなところにあるようなのです。（間違っているかもしれません）

さて、00 00 00 の部分は実際には 00000000　00000000　00000000 を意味します。塩基配列はACGTの4種類(実際にはN等もあるが無視するとして)なので、AAAAが00であろうという推測がつきます。すると、AAACは01、AAAGが02, AAATが03, AACAが04, AACCが05, AACGが06, AACTが07 などという推測がつきます。本当にそうか試してみます。

>chr1 1
AAAA
AAAC
AAAG
AAAT
AACA
AACC
AACG
AACT
AAGA
AAGC
AAGG
AAGT
AATA
AATC
AATG
AATT

とすると

┌────────┬─────────────────────────┬─────────────────────────┬────────┬────────┐
│00000000│ 00 01 02 03 04 05 06 07 ┊ 08 09 0a 0b 0c 0d 0e 0f │0•••••••┊•__•__••│
│00000010│ 00 00                   ┊                         │00      ┊        │
└────────┴─────────────────────────┴─────────────────────────┴────────┴────────┘

を得ることができました。どんぴしゃです。間違いなさそうです。リバースコンプリメントも含めてみます。面倒くさいので seqkit seq -vt dna -w4 -rp poo.fa > poo.revc.fa として、

>chr1 1
AATT
CATT
GATT
TATT
ACTT
CCTT
GCTT
TCTT
AGTT
CGTT
GGTT
TGTT
ATTT
CTTT
GTTT
TTTT

をつくりました。これをパックすると bwa fa2pac poo.revc.fa

┌────────┬─────────────────────────┬─────────────────────────┬────────┬────────┐
│00000000│ 0f 4f 8f cf 1f 5f 9f df ┊ 2f 6f af ef 3f 7f bf ff │•O××•_××┊/o××?•××│
│00000010│ 00 01 02 03 04 05 06 07 ┊ 08 09 0a 0b 0c 0d 0e 0f │0•••••••┊•__•__••│
│00000020│ 00 00                   ┊                         │00      ┊        │
└────────┴─────────────────────────┴─────────────────────────┴────────┴────────┘

２列目がしっかりときれいな数列になっていることに注目してください。やはりここにはリバースコンプリメントがつまっているのでしょう。.bwtファイルは、こちらのリバースコンプリメントが詰まっている方のpacファイルからつくられます。

さて、今回の調査はここまでです。

.bwt ファイルや .sa についても少し調査をしたのですが、理解するためには、BWAの文字列検索アルゴリズムをある程度理解する必要がありそうです。
たぶん私の頭だとわからないんじゃないかという気がしますが、もしも薄っすらと理解できたらまた記事をかくかもしれません。

この記事は以上です。