bwaのインデックスを一部だけながめる ( .amb .ann .pac ) 

## はじめに（おことわり）

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

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

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

```
bwa index reference.fa
```

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

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

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

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

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

http://seqanswers.com/forums/showthread.php?t=25553

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


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

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

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

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

https://github.com/sharkdp/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言語で書きました。

```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](https://en.wikipedia.org/wiki/Nucleic_acid_sequence)をご覧になると、ACTG意外にも非常に多くの記号が入り込む可能性があることがわかります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/144608/210b1981-cc2e-5aeb-225c-29f2a6362b00.png)

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

### 実際にためしてみる。

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

```fasta:poo.fa
>chr1 1
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
```

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

```poo.fa.amb
64 1 0
```

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


```fasta: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を混ぜ込んでみましょう。

```fasta:poo.fa
>chr1 1
NAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
AAAAAAAA
```

こうなります。

```poo.fa.amb
64 1 1
0 1 N
```

こうすると、

```fasta: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
```

のとき

```
128 2 11
0 chr1 1
0 64 1
0 chr2 2
64 64 1
```

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

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

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

```fasta: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の文字列検索アルゴリズムをある程度理解する必要がありそうです。
たぶん私の頭だとわからないんじゃないかという気がしますが、もしも薄っすらと理解できたらまた記事をかくかもしれません。

この記事は以上です。