FM index

シークエンサーのデータを解析するとき、リードをリファレンス配列にマッピングすることが第一歩である。一般にシークエンサーのリードは数百万個以上で、また、リファレンス配列の塩基数は数億にも及ぶ。これらの数百万のリードを、数億塩基からなるリファレンス配列にマッピングする際に、高速な文字列検索アルゴリズムが必要される。このようなアルゴリズムとして、FM index がある。実際に、現在のマッピングプログラムのほとんどは、FM index アルゴリズムを中心に使用してマッピングを行なっている。なお、FM index は完全一致でしか検索できない。そのため、ほとんどのプログラムでは、まずリードの一部分だけをリファレンスに対して完全一致で検索する。そして、完全一致でヒットした箇所で、リードの残った部分をリファレンスに対して Smith-Waterman アルゴリズムでアライメントを行う。

Burrows-Wheeler 変換

FM index は、Burrows-Wheeler 変換(BWT)後の文字列に対して、高速に検索を行うアルゴリズムである。ここで、BWT について紹介する。BWT は長さ L の文字列を受け取り、別の長さの L の文字列に変換する変換アルゴリズムである。変換前と変換後の文字列の長さは変化しないが、変換後の文字列は同じ文字が連続して並びやすくなり、文字列の高速検索アルゴリズムも適用しやすくなる。例えば、次に示した変換前と変換後の塩基配列を見比べると、変換後の配列では、同じ塩基が連続して出現していることがわかる。

before: TACGTAATCATGCGGCTAGCGCATGCATGCGTAGCGCATCGTAGCTC$
after:  CTTTTTACCCCTGGTGGGGAGTGGCCTAATTGACCCCG$GCGCAAAAA

BWT は、入力文字の後ろに終了文字 $ を付加してから、入力文字列を下図のように1 文字ずつ左にシフトして作られる文字列を縦に並べる。続けて、この文字列ブロックの最初 (frist) の列(F 列)をアルファベット順に並べ替える。その後、この文字列ブロックの最後 (last) 列(L 列)を構成する文字列を BWT 後の文字列として出力する。

BWT により文字列を変換する方法

BWT で変換された文字列から元の文字列を復元するために、F 列と L 列の文字列に着目して行う。まず、BWT 後の文字列は L 列であることがわかっている。そして、F 列は、L 列をアルファベット順に並べた文字列に等しい。これで、L 列と F 列の構成が復元できた。続けて、この文字列ブロックは入力文字列を 1 文字ずつシフトしながら作っていることを考慮して、その逆方向でシフトすることにより L 列を F 列の前に移動する。こうして得られる文字列ブロックの先頭の列に対して、並べ替えを行う。先頭列がソートされると、この文字列ブロックの最後の列は、BWT 後の文字列と同じ並び順になる。これで、文字列ブロックの最後の列を復元できる。次に、もう一度シフトして、並べ替えて、最後の列を復元する。この操作を文字列の長さと同じ回数だけ繰り返せば、BWT 前の文字列を復元できるようになる。

BWT により変換された文字列を復元する方法

この逆変換のアルゴリズムは簡単で実装しやすいものの、文字列が長くなると非常に効率が悪い。そこで、現在は、そこで効率のよい逆変換のアルゴリズムとして LF マッピングが考案され、使われている。

LF マッピング

LF マッピングは、L 列にある文字が、F 列のどこにあるのかを調べるアルゴリズムである。例えば、上の例では、F 列が $aaaccg であり、L 列が gc$aaac である。L 列にある g は F 列の 7 番目の g と同一文字であり、LF(g) = 7 とかける。しかし、L 列の 2 番目 c は、F 列の 5 番目の c に対応しているのか、あるいは 6 番目の c に対応しているのかは、明確ではない。この対応関係について調べていく。

この対応関係を調べる時に、BWT のアルゴリズムについてもう一度見ていく。BWT は、入力文字列 acaacg$ を 1 文字ずつ左にシフトして文字列ブロックを作成し、次に文字列ブロップの最初の列について並べ替えることによって行われる。

BWT 変換で L 列と F 列にある文字列の順序対応。

ここで、文字列ブロックのうち、c から終わる 2 行だけに着目する。この 2 行を 1 文字ずつ右にシフトする(BWT とは逆方向にシフトする)ことで、もとの文字列ブロックの 1 つ上の行の文字列を復元できる。これは、L 列ある c の順番は、F 列にある c の順番をそのまま反映していることを表している。

LF マッピングで L 列の文字と F 列の文字のランクが同じである。

このように、L 列にある文字と F 列にある文字の対応関係をわかると、これをコンピュータで扱いやすいように関数化しておく必要がある。ここで、F 列で文字 c が 最初に現れる位置を C[c] とし、F 列の 0 から i 番目の文字の中で c が現れる回数を occ(c, 0, i) とし、さらに L 列の i 番目にある文字を L[i] とすると、LF マッピングは、次のような式で書き表すことができる。

\[ LF(i) = C[L[i]] + ooc(L[i], 0, i) \]

BWT のときに行なった 1 文字シフトと LF マッピングの関係を合わせると、次のような図を作成できる。これらの F 列にある $ から初めて、矢印を辿っていけば文字列 $gcaaca が作られ、これを反転させると acaacg$ が得られる。このように、LF mapping を利用することで、高速に BWT 後の文字列を復元できる。

FM index

FM index は、入力文字列(クエリー)と同じ文字列構成の部分が検索対象文字列(リファレンス)のどこにあるのかを調べるアルゴリズムである。リファレンス配列 acaacg に対して、クエリー配列 caac がどこにあるのかを検索する例を示す。検索のとき、クエリー配列の後ろから前に検索する。まず、クエリー配列の最後の c が、BWT 後の文字列のどこにあるのかを調べる。BWT 後の文字列の c の範囲が絞られると、その範囲のなかで c の後に a が現れる行を調べる。これで、クエリー配列の最後の 2 文字の検索が完了したことになる。続けて、ca の検索で絞られた a に着目して、その a の範囲でさらに ac に対応するものを検索していく。FM index では基本的にこの操作を繰り返すことによって、クエリー配列と完全一致の文字列をリファレンスから見つける。

FM index により文字列検索を行う方法。

完全一致の文字列が見つかると、これを suffix array とよばれる配列と照合して、その位置を出力する。suffix array は BWT を行う途中プロセスで生成できる。suffix array の作り方として、リファレンス配列の後ろに $ をつけて、1 文字ずつ左にシフトして、文字列ブロックスを作る。次に文字列ブロックの最初の列の文字に対してアルファベット順で並べ替える(BWT と同じ)。並べ替えによってできた各行が、並べ替え前の行の何行目であるのかを記述している行列が suffix array である。

FM index により文字列検索を行い、suffix array でヒット位置を調べる。

suffix array のすべての情報を保存することで膨大な容量を占用してしまう。そこで、suffix array 全体を保存せずに、間引きしてから保存する方法が取られている。以下は、suffix array を 1 行ごとに間引きして保存しているれいである。この際に、存在しない行に検索がヒットした場合は、その直前の文字の情報をしらべて、その位置番号から現在のヒット位置の番号を計算することができる。

FM index により文字列検索を行い、suffix array でヒット位置を調べる。

一般に FM index は完全一致でしか検索できない。そのため、ほとんどのプログラムでは、まずリードの一部分だけをリファレンスに対して完全一致で検索する。そして、完全一致でヒットした箇所で、リードの残った部分をリファレンスに対してアライメントを行う。この際にギャップやミスマッチなどを可能にする Smith-Waterman アルゴリズムが使われることが多い。