ナノポアアセンブリツールについて
気が付くとブログを書くのも随分久しぶりですが、年末のAdvent カレンダーという締め切りを持つことで、バラバラになっている知識をまとめて、また理解が足りない点を補足する自身の勉強のためにも書きました。今回は特にナノポアのシークエンスデータに対してよく使われているいくつかのツールについて、そのアルゴリズムの背景を書いてみました。まだ知識不足や書き足りない個所もありますが、誰かの理解の助けになれば幸いです。
ナノポアのデータの特徴は非常に長いリードを得ることが可能であること（ただし入力したデータが短いDNAである場合は、短いリードしか出てきません）、しかしながらエラーがショートリードと比べると高いことなどから、今回は特にこの長いリードを使って、De novo アセンブリのアルゴリズムがどのような工夫を凝らしているのかを各種ツールの工夫とともに見ていきたいと思います。
私は様々なツールのHow it works? (どうやって機能しているのか)、アルゴリズムの中身にいつも興味が向いてしまうので、主に文献やWebなどからわかる範囲でツールの内容の詳細を記載しました。利用方法については、それぞれのツールの開発者のサイトや日本語でもすでに大変良いサイトが用意されているので、今回は特に利用方法については記載していません。それらを期待している方々には大変申し訳ないのですが、ご理解いただき興味を持っていただければ読み進めてください。とても長いです（汗）
De Novo アセンブリは、ゲノムを再構築する手法です。現在市場には、ショートリードとロングリードが存在しており、ロングリードでのDe Novo アセンブリツールはショートリードのアセンブリをそのまま利用することが難しく（計算量が違う、エラーがショートリードよりも多いなど）、少し異なるアルゴリズムでの解析が必要となり、様々なツールが開発されています。多くは、OLC (Overlap Layout Consensus)というリード同士を比較しながらオーバーラップを見つけつなげていく手法が多く採用されていますが、ショートリードで利用されていたde Bruijnグラフを応用したものなどもあり、とても興味深いです。現在、ナノポアでよく使われるDe Novo アセンブリツールとしては、様々ありますが、今回は特に以下の4つについてみていきます。Shastaは比較的新しいアセンブリツールですが、それ以外は割とよく利用される選択肢の一つになっているかなと思います。
- Canu
- Flye
- Redbeans (wtdbg2)
- Shasta
Canu: https://github.com/marbl/canu
Canu (Koren et al., 2017) はヒトゲノムプロジェクトなどで用いられていたCelera Assembler の後継としての開発され（Celera Assemblerの開発は現在停止）、ロングロード用のアセンブリツールとして開発が進んでいます。CanuのアルゴリズムはCorrect, Trim, Assembly という３つのパートからなっている。これらは、一部だけ実施することも可能となっています。論文に以下のようなアルゴリズムの流れの図があります。
それぞれのステップにあるgkpStore がリードを保存し、ovlStoreがオーバーラップした情報を保存しています。それぞれのステップで使うオーバーラップを見つける方法は同じですが、パラメータの設定がそれぞれのステップごとに少し違っています。Supplementに記載のある以下の図は、gkpStoreとovlStoreの関係を理解するのに役立つと思います。
Canuでは、著者の一人にもなっているBerlinが提案したMinHash Alignment Program (MHAP) を使い高速にオーバーラップを検出できる方法を提案しています (Berlin et al., 2015)。このMHAPは、Webのテキスト検索などで用いられるMinHash (Broder, 1997)という方法をシークエンスデータに応用したものになります。MHAPの論文に良い図があるので、その図を使ってみていきましょう。
まず最初のステップ（A）ではリードがk-mersの長さに分割される。この例では3塩基（実際はもっと長いです）を示している。これにより、14塩基のリードから12種類のk-merが作成された。次の（B）ステップでは、これら12種類のk-mersを4種類のハッシュ関数で数値化しています。Wikipediaによるとハッシュ関数とは、ある値が与えられた時にそれを代表する数値を得る操作、またはその値を得るための関数とあります。おそらく身近なハッシュ関数による値は、ファイルをダウンロードした時に与えられるMD5ではないかと思います。ファイルのサイズだけではなく、このMD5が一致することを確認して、初めてそのファイルが確実に同じだと確認するように、ハッシュ関数によりある情報を要約した数値列を得ていると考えます。長くなりましたが、MHAPでは、複数のハッシュ関数を使い、k-mersから要約したハッシュ値を得ます。（C）つぎのステップCでは、それぞれのハッシュ関数で作成されたハッシュ値が一番小さいものを集めたSketchがS1,S2それぞれに作成されます。（D）そして、このうち、最小のハッシュ値が小さいものが、S1とS2で同じものがいくつあるかを調べます。この一致する最小のハッシュ値を持つ数が２つの配列を比較したときの一致率を計算すると、これは集合の類似度を計算するJaccard係数に一致するというのが、MinHashの重要な性質になります。これを応用してリードの配列比較アルゴリズムMHAPは作成され、それを組み込んだのがCanuという訳です。このMinHashとJaccard関数の関係については、すぐにすんなり理解できないかもしれないので、興味のある人はこれら（https://tech.preferred.jp/ja/blog/minhash/ 、
https://www.youtube.com/watch?v=bQAYY8INBxg）を参考にしてもらうといいと思います。
最後のステップの（E）では、一致したk-mers をもとの配列S1,S2に戻し、その配列での類似度を計算し、オーバーラップのOffsetを計算するというステップになります。
Canuではさらに、リピート配列に対して柔軟に対応できるように、ゲノム中に頻繁に表れるk-merは繰り返し現れるため、重みを付けて値を低くし、逆に特徴的なk-merは重要に扱われるように重みを付けています。これを実現するために画像データの検出に使われたMinHashの拡張、tf-idf を使った重み付きのMinHash (Chum, Philbin, & Zisserman, 2008)が採用されています。すごいですね。
Canuの性質としてSparse Graph の構築がありましたが、ちょっと時間切れです（ごめんなさい）。冬休みに追記出来たら、追記します！
Redbean (wtdbg2) https://github.com/ruanjue/wtdbg2
もともとの名前がwtdbg2という非常に覚えにくい名前であったため、名前がRedbeanと変わったことに安堵したバイオインフォマティシャンは多くいたと思います。このツールは、中国のAgricultural Genomics InstituteチーフサイエンティストをしているJue Ruanが開発を行い、共著にはHeng Liが入っています。Redbeanの原理はショートリードのDe Novo アセンブリで利用されていたDe Bruijnグラフを応用したFuzzy De Bruijnグラフとなり、k-mer単位で作るde Bruijnグラフをk-bin単位というもっと大きなサイズでグラフを作成し、高速にアセンブリが行えるようになっています。Binの長さは256bpとなっていて、オーバーラップさせずに、設定します。オーストラリアのクイーンズランド大学で毎年行われているWinter Schoolというバイオインフォの冬の学校で、Jue Ruan が招待され、その時のスライドに詳細が記載されているので、こちらを使ってみていきましょう。
まずk-binの作られ方ですが、1本のリードを256bpに分割します。通常、de Bruijnグラフでは1塩基ずつずらしたk-merを作成するため、ここがde Bruijnグラフとの大きな違いになります。理由としては多くのk-merは重要ではなく、計算量とメモリを使っていること、またエラーが多いため短くカットしたk-merにエラーが多く含まれるものが出てくるから、という事のようです。
全てのリードを256bpに分割した後は、256bpの単位でDynamic Programming を行いオーバーラップしている箇所を決めていきます。
RedbeanのLimitationとして、リードの最大長が256kbとなっているようで、それより長いリードは現在はカットして利用するというようにマニュアルには記載があるので、この点は非常に長いリードを入れて、それによるScaffoldを期待している人は要注意かもしれませんね。
https://github.com/ruanjue/wtdbg2/blob/master/README-ori.md
Flye: https://github.com/fenderglass/Flye
それではFlye (Kolmogorov, Yuan, Lin, & Pevzner, 2019)を見ていきましょう。これはアメリカのUCSDで開発された手法になります。Flyeを理解するには、GithubのUsageのところにAlgorithm概要があり、これが一番シンプルで分かりやすかったので、その内容に沿って記載したいと思います。
Flyeはまず最初にドラフトのコンティグを作ります。まずリードをk-merへ分割し、明らかにエラーと思われるようなk-merは利用せずにオーバーラップを見つけドラフトのコンティグを作ります。ここでは大雑把なコンティグで間違いなどは後で修正することを想定しているようです。
そして次にリピートの解析を行います。まずリードがリピートの両端をまたぐようなBridged リピートを解消し、ドラフトのコンティグからリピートグラフを作成します。Githubにリピートグラフの図があります。De Bruijnグラフよりも、もう少し要約しているようなグラフですが色のついているエッジがリピートです。
そして、Bridgedのリピートが解消された後、Bridgedされていないシンプルなリピート（繰り返しが2回）について、Trestle モジュールがリピートを解決できるかリードの多様性などから解析するようですが、Trestle モジュールについてはあまり記載がありませんでした。
そしてFlyeは最後のステップとして、残りのエラーを補正するためにminimap2を使い、すべてのリードをコンティグに張り付けてエラーを補正していくという作業を数回繰り返します。
SHASTA: https://github.com/chanzuckerberg/shasta
Shastaは11人のヒトゲノムを9日間で解析したという論文とともにアナウンスされたde Novo assemblyツールになります (Shafin et al., 2019)。アルゴリズムの基本はOverlap Layout Consensus となりますが、ホモポリマー領域のエラーをうまく回避するためにRun Length Encoding (RLE) によるホモポリマー領域の記述、そのほかのエラーに対してもオーバーラップしているリードの個所を見つけるために、一部のk-merのみをマーカーとして利用し、そのマーカーの一致からリードのオーラ―ラップを見つけていきます。類似度の計算はMinHashを使っています。
オンラインComputational Methodsに詳細が記載されているので、その中の図を使いながら見ていきましょう。
https://chanzuckerberg.github.io/shasta/ComputationalMethods.html
RLEはホモポリマーのエラー領域を改善するために利用されているEncoding方法です。たとえば、AAAAとなっている箇所をA4と記載するような手法です。リード配列の一部で例を作ると以下のようになります。
現在のナノポアのデータで最も多いエラーのタイプはホモポリマー領域に見られるので、その対策ですね。
RLEを使うことで、ホモポリマー領域のエラー対策は行われていますが、それ以外の領域のエラーについても、Shastaは工夫をしています。オーバーラップの検出にCanuでも使われていたMinHashを使っていますが、リード全体の情報から作成するすべてのk-merを使うのではなく、Marker というランダムに選んだk-merで部分的な配列情報をオーバーラップの有り無しに判定するために利用しています。
つまりリードにそのk-merを照らし合わせると、以下のような感じです。
ちょっとびっくりですが、ある程度の長さを想定しているのであれば、確かにすべての配列情報を使わず、一部だけ使ってオーバーラップを見つけるのは、なるほどなと思ったりしました。図では3塩基ですが、Defaultは10塩基で、これはパラメータで変更可能です。
このRLEとMarkerをリードに当てはめた図をみてみましょう。これを見るとかなりすっきりですね。黒い塩基配列の上に表示されているのがRLEですね。
Shastaの利用にあたっては、想定されているデータについても知っておくほうが良いでしょう。筆者らのチームは非常にナノポアのシークエンスデータ取得に慣れているため、リード長の多くが50kbであること、そして10kbのリードは10％以下で、それよりも多い100kb以上のリードがあること、などを想定しているデータとして挙げています（https://chanzuckerberg.github.io/shasta/ComputationalMethods.html）。これらのデータを得るためには、論文で記載されているように、シークエンスにかける前に短い断片をCirculomicsのShort Read Eliminatorで取り除き、さらに長いDNA断片をピペッティングなどで切ってしまわないように注意を払ってライブラリ調整を行う事が必要となります。Shastaを使う際には、自分のデータがこの筆者らのデータとどの程度似ているのかを知ったうえで使うほうが結果の解釈の役に立つかもしれませんね。
さてさて、一体誰得なの？という内容のわりに随分長くなってしまいましたが、この記事を書くまでの調べ物や文献読みは大変楽しかったです。もう少しグラフ構築については勉強しないとなという点など浮き彫りになりました。ツールの中身に興味があるけど、なかなかじっくり理解する時間が取れないな、というような誰かの役に立っていれば幸いです。
最後にこのような良い機会を与えていただき、企画くださったKurokiさんに感謝いたします。
引用文献
Berlin, K., Koren, S., Chin, C.-S., Drake, J. P., Landolin, J. M., & Phillippy, A. M. (2015a). Assembling large genomes with single-molecule sequencing and locality-sensitive hashing. Nature Biotechnology, 33(6). https://doi.org/10.1038/nbt.3238
Berlin, K., Koren, S., Chin, C. S., Drake, J. P., Landolin, J. M., & Phillippy, A. M. (2015b). Assembling large genomes with single-molecule sequencing and locality-sensitive hashing. Nature Biotechnology, 33(6), 623–630. https://doi.org/10.1038/nbt.3238
Broder, A. Z. (1997). On the resemblance and containment of documents. Proceedings of the International Conference on Compression and Complexity of Sequences, 21–29. https://doi.org/10.1109/sequen.1997.666900
Chum, O., Philbin, J., & Zisserman, A. (2008). Near Duplicate Image Detection: min-Hash and tf-idf Weighting. Procedings of the British Machine Vision Conference 2008, 27(5), 50.1-50.10. https://doi.org/10.5244/C.22.50
Kolmogorov, M., Yuan, J., Lin, Y., & Pevzner, P. A. (2019). Assembly of long, error-prone reads using repeat graphs. Nature Biotechnology, 37(5), 540–546. https://doi.org/10.1038/s41587-019-0072-8
Koren, S., Walenz, B. P., Berlin, K., Miller, J. R., Bergman, N. H., & Phillippy, A. M. (2017). Canu: scalable and accurate long-read assembly via adaptive k -mer weighting and repeat separation. Genome Research, 27(5), 722–736. https://doi.org/10.1101/gr.215087.116
Shafin, K., Pesout, T., Lorig-Roach, R., Haukness, M., Olsen, H. E., Bosworth, C., … Paten, B. (2019). Efficient de novo assembly of eleven human genomes using PromethION sequencing and a novel nanopore toolkit A Preprint. https://doi.org/10.1101/715722
