レトリバのCTOの武井です。
https://twitter.com/goth_wrist_cut
新型コロナウィルスが世界で猛威を振るっていますが、皆様安全に過ごせておりますでしょうか。 レトリバでは フルリモート化 や、 交流なども オンライン飲み会 にするなど、工夫して過ごしています。
さて、今回はそんな新型コロナウィルス、COVID-19の遺伝子配列をターゲットに、 Longest Common Substring(最長共通部分文字列)を求めたり、そのアルゴリズムの解説をしてみようと思います。
アルゴリズムの説明自体は William Fiset さんという方が動画で説明しており、 非常に分かりやすい図示がありますので、是非ご覧になって下さい。
https://www.youtube.com/watch?v=Ic80xQFWevc
https://www.youtube.com/watch?v=DTLjHSToxmo
このブログでは日本語での解説と、実際の実装まで踏み込んで紹介しています。
はじめに
さて、世界で猛威を振るっているCOVID-19ですが、 その遺伝子情報がNCBI(アメリカ国立生物工学情報センター)に登録されており、ダウンロードできるようになっています。
2020-05-27現在、COVID-19( SARS-CoV-2 ) の配列のうち、概要にcompleteが含まれるもの( Complete )で 2883 種登録されています。
最初に試しにデータを触ってみた、2020-04-13時点では700ぐらいでしたので、ここ一ヶ月半で登録が四倍以上になっているようです。 二本鎖のDNAに比べ、一本鎖のRNAしかもたないRNAウィルスは変異しやすい、というのもありそうですが、今まさに多くの人が研究をしていることが分かります。
さて、そんな変異しやすいRNAですが、逆に変異していない部分はあるのでしょうか？
参考として、HIVウィルスの一つ、HIV-1 ウィルスでは、ほぼ全域に渡って変異が起きているようです。
http://data.dbcls.jp/~meso/meme/phd/p3/ (DBCLS 内藤雄樹 理学博士による研究)
変異していない部分、つまり全ての変異において保存されている部分が見つかれば、 たとえば、その領域を対象にした遺伝子創薬などができるかもしれません。
そこで、世界を救う文字列アルゴリズム、というのはいささか大げさな表現ですが、 バイオインフォマティックスの分野で活躍している文字列アルゴリズムの中から、 今回は複数の遺伝子に対し、変異していない最長領域を見つけることのできる、 Longest Common Substrings を求めるアルゴリズムを紹介したいと思います。
ちなみにではありますが、塩基配列のデータは A C G T の塩基を表す記号のほか、
「C または G」を表す S や、 「A 以外」を表す B など、あいまいな塩基を表す記号が含まれていることがあります。
漏れの無い検索などでは記号の包含関係を含めた比較が必要となるのですが、今回のLCSのプログラムでは文字の一致のみをみています。
今回利用したデータでは、あいまいな塩基を表す記号が含まれているデータは、2883件中896件。
あいまいな塩基を表す記号がもっとも多いデータは、29871塩基中1451塩基が N (任意の塩基を表す)となっているデータでした。
Longest Common Substrings(最長共通部分文字列)とは
読んで字の通り、複数のデータに対し、共通して出現する部分文字列のうち、最長のものを指します。
実例で見ていきましょう。 COVID-19の配列は1データあたり3万塩基ぐらいあって、そのまま例に使うとキビしいので、短いサンプルを使って説明します。
データ1: CATTTACG データ2: ACACACATTT データ3: GCATATTT
この三つのデータに共通して出てくる部分文字列を探してみます。
たとえば、データ1の先頭三文字の CAT は他のデータにも含まれているため、共通部分文字列になっていそうです。
データ1: CAT TTACG データ2: ACACA CAT TT データ3: G CAT ATTT
CAT の前後は各データそれぞれ別の文字になっており、
CAT を伸ばしてそれより長い共通部分文字列を見つけることはできなさそうです。
ところが、この CAT から C をとって、 TT を足した ATTT は、
データ3では異なる位置に含まれているため、 ATTT も共通部分文字列となります。
データ1: C ATTT ACG データ2: ACACAC ATTT データ3: GCAT ATTT
今回のデータの場合、 ATTT が共通部分文字列のうち最長となり、最長共通部分文字列(Longest Common Substrings)となります。
さきほどの例から、一度見つけた共通部分文字列を延長しても局所解となってしまうことがわかりましたので、 最長共通部分文字列を求めるためには適切なアルゴリズムが必要そうです。
COVID-19の遺伝子配列で実際に求めた場合、2020-05-27現在 (2883種)での最長共通部分文字列は、
TCAGCTGGTTTTCCATTTAATAAATGGGGTAAGGCTAGACTTTATTATGATTCAATGAGTTATGAGGAT CAAGATGCACTTTT
の83文字(83塩基)でした。
100文字近いパターンですので「全パターン作って検索」みたいなアプローチでは、 文字種4つしかない遺伝子データでも、4の100乗で60桁ぐらいですから、一般的なケースだともはや話にならなさそうです。
全パターン作るのではなく、どれか一個のデータを m 文字ずつ取り出して(m-gram)、クエリにして検索する、
みたいな方針ですと、一回の検索がデータの平均長 N に対して、BM方法で O(N/m) 、データ個数が M 個なら O(N*M/m) 、
パターンをずらして作っていくので N-m+1 回クエリを試すことになり、 O(N*M/m*N) 。
検索したあとは、全データに共通して出現する隣り合う m 文字のクエリをまとめあげれば答えを出せそうですが、
計算量としては平均データ長の二乗のオーダーとなりそうです。
(ちなみに m は大きくするほど速度が速くなりますが m より短いLCSを求められなくなります)
COVID-19の場合、平均データ長は3万ぐらいなのでまだなんとかなりそうでしょうか。
なお、LCSの問題のバリエーションとして、データ全てに出現するのではなく、 少なくともK個のデータに出現する最長共通部分文字列を探す K-Longest Common Substrings という問題設定もあります。 変異が大きすぎるデータをいくつか無視しつつも、大多数での共通部分を探す、という感じですね。
こうなると、先ほどのアルゴリズムだと、クエリ生成に選んだデータが除外されるデータの可能性がありますので、
先の操作をデータ数M個に対して少なくとも M-K+1 回繰り返す必要がでてきます。
あるいは先に、出現する全ての m-gram を列挙してしまう、という手もありそうです。
今回のブログエントリでは、検索でおなじみのSuffix Arrayを使った、 Longest Common Substrings(および K-Longest Common Subsrings)を求めるアルゴリズムを紹介します。
なお、以下 Longest Common Substrings と打つのが大変なので、LCSと略そうと思いますが、 一般的にLCSと書いた場合、よく似た違う問題、Longest Common Subsequence(最長共通部分列)を指すことが多いようです。
Longest Common Subsequenceは、同様に全てのデータに出現するパターンを見つける問題ですが、 Substrings(部分文字列)と違い、連続していなくてもよく、前後関係だけ保った「Subsequence(部分列)」を見つける問題になります。
たとえば、冒頭の例題ですと CATTT が最長共通部分列となります。
C が ATT と離れていますが、各データ、この順で出現することが分かります。
データ1: C ATTT ACG データ2: ACA C AC ATTT データ3: G C AT ATTT
(上図の位置合わせは一例で、他の位置の合わせかたもできます。Subsequenceは「出現順」だけを気にします)
Longest Common Subsequenceもまた、バイオインフォマティックスでよく出てくる問題で、複数の遺伝子配列のアライン、 つまり突然変異(遺伝子の置換挿入欠損)があった部分を除いた箇所を合致させたりする時などに利用されます。
いずれLongest Common Subsequenceの話もできれば、とは思いますが、 今回はSubstringsについてお話しようと思いますので、 LCSをLongest Common Substringsの意で用います。 同様に、少なくともK個のデータに共通な最長部分文字列、K-Longest Common Subsringsも K-LCSと表記します。
Suffix Array と LCS
さて、部分文字列を探す話はSuffix Arrayと相性がいい問題です。 Suffix Arrayは全ての接尾辞をソートして並べたものになりますので、同じ部分文字列は同じ箇所に集まります。
冒頭のデータをサンプルに、実際にSuffix Arrayを作成してみるのですが、 Suffix Arrayを作るためには、一個のデータになっている必要があるため、 データとデータの間に区切り文字をいれて連結してからSuffix Arrayを作ります。
区切り文字は、通常のSuffix Arrayを作る際の終端文字 $ と同様に、他の文字より小さい文字として扱います。
一般的にこういう場合、区切り文字は終端文字と区別し # をつかうことが多いですが、
ちょっと面倒なので終端文字と同じ $ を使いました。
(# は結合する、というニュアンスで使われることが多いです)
データ1: CATTTACG データ2: ACACACATTT データ3: GCATATTT
上の元データから、一つの文字列 CATTTACG$ACACACATTT$GCATATTT$ を作り、Suffix Arrayを作ります
28 $ 8 $ACACACATTT$GCATATTT$ 19 $GCATATTT$ 9 ACACACATTT$GCATATTT$ 11 ACACATTT$GCATATTT$ 13 ACATTT$GCATATTT$ 5 ACG$ACACACATTT$GCATATTT$ 22 ATATTT$ 24 ATTT$ 15 ATTT$GCATATTT$ 1 ATTTACG$ACACACATTT$GCATATTT$ 10 CACACATTT$GCATATTT$ 12 CACATTT$GCATATTT$ 21 CATATTT$ 14 CATTT$GCATATTT$ 0 CATTTACG$ACACACATTT$GCATATTT$ 6 CG$ACACACATTT$GCATATTT$ 7 G$ACACACATTT$GCATATTT$ 20 GCATATTT$ 27 T$ 18 T$GCATATTT$ 4 TACG$ACACACATTT$GCATATTT$ 23 TATTT$ 26 TT$ 17 TT$GCATATTT$ 3 TTACG$ACACACATTT$GCATATTT$ 25 TTT$ 16 TTT$GCATATTT$ 2 TTTACG$ACACACATTT$GCATATTT$
さて、LCSの ATTT や最長候補だった CAT などを見ると、やはり同じ箇所にまとまっていることが分かります。
このSuffix Arrayを走査して、横幅が最大となるような「同じ文字列が並んだ長方形」部分、を見つければよいのですが、 ただ単に並んでいればいいわけではありません。
LCSの制約条件
先ほどのSuffix Arrayの図に、長方形を足しました。
下二つの水色の枠は、これまでにでてきた CAT や ATTT でLCSの候補です。
一方、上の黄色で囲った ACA も CAT や ATTT と同様に、同じ文字列が3つ並んでいます。
これもLCSの候補でしょうか？
これは、元のデータを見直すと、いずれもデータ2 ACACACATTT の前半部分から出てきていることが分かります。
共通部分文字列は異なるデータで共通な文字列であって、同じデータに複数発生する文字列はカウントしたくありません。
そこで、Suffix Arrayの図に、データ番号(docid)と、部分文字列がどのデータ由来なのかを分かり易くするため色を付けてみました。
区切り文字と終端文字はデータ番号なし、にしています。
ACA はデータ番号が 2 の一種類しかなく、複数のデータをカバーしていないことがわかります。
一方で、 CAT や ATTT は 1 2 3 全てのデータ番号を含んでいることが分かります。
共通部分文字列になるのは、データ番号を「全種類含む」長方形、となりそうです。 K個以上のデータに共通して出現するK-LCSの場合には、全種類ではなく「K種類以上含む」という条件を使えばよさそうです。
この制約を満たしつつ、長方形の横幅、つまり共通部分文字列の長さを最大化したい、となります。
共通部分文字列の長さ
共通部分文字列の長さは範囲内の各文字列の共通の接頭辞(prefix)の長さとなります。 文字列が二つ三つぐらいであれば文字列を横断していけばよいですが、データ数が100や1000になると大変です。
そこで、隣り合う二つの文字列で先に共通接頭辞の長さを計算しておいて、必要な範囲で結果をまとめるようにします。 文字列が三つ以上の場合は、隣り合う二つの文字列の共通接頭辞の長さの最小値が、全体での共通接頭辞の長さとなります。
共通接頭時の長さとして、この図にLCPというカラムを足します。 LCPは、Longest Common Prefixの略で使っており、一つ手前の文字列とその文字列との共通接頭辞の長さをいれています。
例えば、 Suffix Arrayのテーブルの ATTTACG$ACACACATTT$GCATATTT$ の行場合、
一つ手前の行は ATTT$GCATATTT$ ですので、共通接頭辞は ATTT で、LCPは4になります。
ATTT の場合、LCPのカラムは 5 と 4 で、そのうち小さい方の 4 が、
CAT の場合は 3 と 5 のうち小さい方の 3 が、長方形の幅 = 共通部分文字列の長さとなることが分かります。
一つ注意点としては、LCPは、今見ているところと手前、の二つを見ているため、先頭のLCPは使いません。
データ数 M に対して、隣り合う二つずつの「あいだ」は M-1 個、いわゆる「植木算」なことに注意です。
前述の制約を満たしつつ、この「共通部分文字列の長さ」=「LCPの最小値」を最大にすることができれば、LCSを求めることができそうです。
「LCPの最小値」は、「最小値」ですからLCPの個数が減るほど、値が大きくなるはずです。 LCPの個数は範囲の中間部分を飛ばすわけにはいかないですので、「長方形」の高さを縮めるしか、ありません。
つまり、「制約を満たしつつ、最小の連続した範囲を得る」 言い換えれば「なるべく少ない連続したデータで制約条件を満たす」ことができればLCSの候補を求めることができます。
Slide Window
制約を満たしつつ、最小の連続した範囲を探したい時には、制約付きSliding Windowが使えます。 競技プログラミングでは「しゃくとり法」(two pointers)と呼ばれているようです。
制約付きSliding Windowは、単純には「制約条件を満たすまではWindowの拡大する」 「制約条件を満たしたら制約条件を満たさなくなるまで縮小する」を繰り返す可変幅のSliding Windowです。
「拡張して制約条件を満たした時」ではなく「縮小して制約条件を満たさなくなる 直前 」が、 求めたい「制約を満たしつつ、最小の連続した範囲」となります。
実際にやってみましょう。
スタートは区切り文字や終端文字を除いた部分から開始します。
この時点では、データ番号は 2 のみですから条件をみたしません。
制約条件、docidが全種類(1,2,3)揃うまで範囲を拡張します。
この時点では制約条件は満たしていますが、最小になっていません。 条件を満たさなくなる直前まで範囲を狭めます。狭めるときは上端を縮めます。
制約条件を満たす、最初の最小セットが手に入りました。
LCPは前述の通り、一件目を無視するので、 2 と 1 のうち小さい方、 1 がLCSの長さとなり、 LCSの候補 A を得ます。
拡張をする前に、もう一つ進めて制約条件を満たさなくなるまで縮めます。
制約条件が成り立たなくなった状態から、 また「制約条件を満たすまで拡張」、制約条件を見たらしたら「条件を満たさなくなる直前まで縮小」を繰り返します。
今回は制約条件を満たすまで拡張した時点で、一つでも縮めるとデータ番号 3 が無くなってしまうため、
範囲は既に最小となっています。
データ番号 1 のデータが二つあって無駄なように見えますが、これ以上縮小するとデータ番号 3 が欠損してしまいます。
LCPは、 1 2 5 のうち最小の 1 で、再びLCSの候補は A です。
もう一つ縮めて、また繰り返しです。
制約条件を満たさなくなるまで縮めたら、また制約条件を満たすまで拡張します。
拡大して制約条件をみたしたら、制約条件を満たさなくなる直前まで縮小。
LCPは 5 と 4 のうち最小値は 4 で、LCSの候補 ATTT を得られました。
以下同様に繰り返すことで、 CAT や ATTT の部分文字列 TTT なども同様に得ることができます。
画像が大きくて全体の流れが分かりづらかったので動画も作りました。
道具立て
さて、アルゴリズムの説明としては以上なのですが、このアルゴリズムを実装するにあたって、必要な道具がいくつかあります。
一つ目は制約条件をチェックするための、異なり数をカウントするキュー。 Sliding Windowが伸び縮みしている間、 その範囲にデータ番号が全種類含んでいるか？(K-LCSでは K種類以上含んでいるか？)をチェックできるキューが必要となります。
二つ目は最小値キュー。 同様に、Sliding Windowが条件を満たしたときにLCPの最小値を取り出したいわけですが、 都度都度さかのぼって最小値を計算しなおすのは面倒です。 Sliding Windowが後戻りをしないアルゴリズムなので、最小値も後戻りをしない方法でとれると嬉しそうです。
最後にそれらを利用した、SlidingWindowを部品として作っていきます。
CardinarityQueue(異なり数カウントキュー)
キューの中の要素の種類数(異なり数)を数えるキューの実装です。
難しいことはなにもなく、ただのキューに、連想配列と連想配列のうち非ゼロな要素の個数を数えるカウンターを足しただけです。
要素の追加(Push)によって、ある要素の個数が非ゼロになったら異なり数を +1 、
要素の削除(Pop)によって、ある要素の個数がゼロになったら異なり数を -1 、としているだけです。
今回は連想配列で実装していますが、数えたい対象がデータ番号なので、配列でも実装可能です。
配列で実装すれば要素アクセスは定数時間ですので、 Push Pop Cardinality いずれも O(1) となります。
MinQueue(最小値キュー)
キューの中の要素のうち最小値を返すキューの実装です。
プログラミング面接なんかでもおなじみの課題ですね。 Min付きStackを実装して、Stack二つでキューを実装する方針もありますが、今回は両端キューを使った実装をしてみています。
最小値を保持するキュー min_queue_ は、通常のキューと同じように末尾追加、先頭削除をしますが、
制約として内部のデータが常に広義単調増加となるようになっています。
つまり、末尾追加するさい、追加する新しい要素より小さい既存要素を全て末尾から削除し、新しい要素を追加します。
また、先頭削除するさいに、通常のキューの先頭と最小値キューの先頭が同じ値の場合には、最小値キューから先頭削除をします。
これで最小値キューの先頭が、キューの要素の最小値となります。
イメージ的には、野球部の新一年生にエースが入ってきたら、 それより年上の選手は先に卒業してしまうので、卒業までトップになれないので監督のエースメモから除名、 新一年生は、先輩にエースがいたとしても、その先輩が卒業すればトップなるかもしれないので、 とりあえず監督メモには残留(ただし次の新一年生に抜かされたら除名)、みたいなイメージでしょうか。
……ちょっといい例えではなかったですね。
計算量としては Pop Front Min は自明に O(1) 、
Push は時々両端キューを後ろから掘り返すのでコストが掛かりそうに見えますが、
各要素に対して高々一回ずつしか発生しませんので、 N件のデータに対して O(N) になりますので、
一件あたりにすると O(1) となります(いわゆる、ならし計算量ですね)。
SlidingWindow
最後はSlidingWindowの実装です。
CardinarityQueueとMinQueueを使っています。
コンストラクタの least_docnum が K-LCSの K の値、
つまり少なくとも least_docnum 個のデータに共通して出てくる部分文字列を探します。
あとはアルゴリズムの説明通り、拡大と縮小を繰り返すのですが、
説明では「制約条件を満たさなくなる直前まで縮める」という説明をしていましたが、
LCPが先頭一個を省いて計算するため、制約条件を満たさなくなるまで縮めてから(= 一回多く Shrink を呼んでから)、
LCPの最小値を(先頭一個を省かず)求めています。
計算量は各キューへの追加は前述の通り O(1) です。
途中 Shrink をループで呼んでいますが、これは高々追加した要素数分しか呼ばれませんので、これもならして O(1) となります。
コードと解説
さて、必要な道具が揃ったため、全体の実装です。
実験では社内に転がっている大規模検索向けのSuffix Arrayライブラリを利用しましたが、 サンプルコードでは手元で実行可能なように、愚直にsortを使ったSuffix Array構築をしています。
ある程度の文量なら耐えられますが、文章量が多くなると単純なソートだと太刀打ちできなくなってきます。 SA-ISなどのアルゴリズムを使えば線形でSuffix Arrayが構築可能ですが、今回の本筋とずれてしまいますので省略します。
また、文書数があまり多くないため GetDocID は線形探索していますが、文書量によっては二分探索にするなどの工夫が考えられます。
コードとしては単純にSuffix Arrayを順に走査して、隣り合うテキスト同士のLCPを求めてSlidingWindowに足してくだけです。
Suffix Arrayの構築は、SA-ISを使うと文章量に線形の O(N * M) (N: 平均文書長、 M: 文書数、つまり総文字長オーダー)で構築することができます。
ループはSuffix Arrayのサイズ分回って、中の操作は、
SlidingWindowの操作はが O(1) 、 GetDocID は二分探査なら O(log M) となりますので、全体で O(N * M * log M)
最終的な計算量は O(N * M * log M) となります。
一般的にはデータ長 N が大きく、 M はそこそこ、となりますので N の二乗よりはよくなっていそうです。
実際のCOVID-19の配列(84MB / 2883種)での実行時間は、Suffix Arrayの構築で4分程度、それを使った探索が15分程度となりました。
まとめ
今回は、Suffix Array から Longest Common Substrings (およびK-Longest Common Substrings)を求める方法を紹介しました。
COVID-19のSuffix Arrayは、もともと、 レトリバがエンジン提供をしている 超絶高速ゲノム配列検索 GGGenome で検索用に作られていたデータでした。 今回、GGGenomeにCOVID-19の遺伝子情報が検索対象として加わったことをきっかけに、 文字列アルゴリズムでなにか役に立てないだろうか、と思い試しに実装してみました。
Suffix Arrayで共通な文字列を探すアルゴリズムは自分でも試行錯誤してみたことがあるのですが、 今回、SlidingWindowでの制約条件つき探索など、興味深いアルゴリズムに触れることができました。
COVID-19は、まだ日が浅いためか、変異を起こしている部分が少なく、 求めたLCSもまだ直接役に立つようなフェーズでは無さそう、という感じのようです。 私個人はバイオ関係の知識は薄く、治療薬の開発やウィルスの解析などのような直接的な貢献はできませんが、 いち文字列アルゴリズム好き人間として、何か役に立てればと思ってます。
直近では、先ほど名前をだしました、 超絶高速ゲノム配列検索 GGGenome のパッケージ版などをリリースしております。
COVID-19の配列なども、カスタムインデックスとして追加可能ですので、もしご興味ありましたら是非お問い合わせ下さい！
新型コロナのいち早い収束を心より願っております。
