ナノポアアセンブリツールについて

気が付くとブログを書くのも随分久しぶりですが、年末の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

Canuhttps://github.com/marbl/canu

Canu (Koren et al., 2017) はヒトゲノムプロジェクトなどで用いられていたCelera Assembler の後継としての開発され(Celera Assemblerの開発は現在停止)、ロングロード用のアセンブリツールとして開発が進んでいます。CanuのアルゴリズムはCorrect, Trim, Assembly という3つのパートからなっている。これらは、一部だけ実施することも可能となっています。論文に以下のようなアルゴリズムの流れの図があります。

f:id:marimiya_analysis:20191223065041p:plain

論文 (Koren et al., 2017)より抜粋

それぞれのステップにあるgkpStore がリードを保存し、ovlStoreがオーバーラップした情報を保存しています。それぞれのステップで使うオーバーラップを見つける方法は同じですが、パラメータの設定がそれぞれのステップごとに少し違っています。Supplementに記載のある以下の図は、gkpStoreとovlStoreの関係を理解するのに役立つと思います。

f:id:marimiya_analysis:20191223070535p:plain

論文 (Koren et al., 2017)より抜粋

Canuでは、著者の一人にもなっているBerlinが提案したMinHash Alignment Program (MHAP) を使い高速にオーバーラップを検出できる方法を提案しています (Berlin et al., 2015)。このMHAPは、Webのテキスト検索などで用いられるMinHash (Broder, 1997)という方法をシークエンスデータに応用したものになります。MHAPの論文に良い図があるので、その図を使ってみていきましょう。

f:id:marimiya_analysis:20191223070605p:plain

論文 (Berlin et al., 2015b)より抜粋

まず最初のステップ(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で同じものがいくつあるかを調べます。この一致する最小のハッシュ値を持つ数が2つの配列を比較したときの一致率を計算すると、これは集合の類似度を計算する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 が招待され、その時のスライドに詳細が記載されているので、こちらを使ってみていきましょう。

http://bioinformatics.org.au/winterschool/speaker-items/long-read-bioinformatics-and-applications-ruan/

 まずk-binの作られ方ですが、1本のリードを256bpに分割します。通常、de Bruijnグラフでは1塩基ずつずらしたk-merを作成するため、ここがde Bruijnグラフとの大きな違いになります。理由としては多くのk-merは重要ではなく、計算量とメモリを使っていること、またエラーが多いため短くカットしたk-merにエラーが多く含まれるものが出てくるから、という事のようです。

f:id:marimiya_analysis:20191223222136p:plain

スライドより抜粋

全てのリードを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グラフよりも、もう少し要約しているようなグラフですが色のついているエッジがリピートです。

f:id:marimiya_analysis:20191223222417p:plain

Githubより抜粋

そして、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と記載するような手法です。リード配列の一部で例を作ると以下のようになります。

f:id:marimiya_analysis:20191223222759p:plain

Computational Methodsより抜粋

現在のナノポアのデータで最も多いエラーのタイプはホモポリマー領域に見られるので、その対策ですね。

RLEを使うことで、ホモポリマー領域のエラー対策は行われていますが、それ以外の領域のエラーについても、Shastaは工夫をしています。オーバーラップの検出にCanuでも使われていたMinHashを使っていますが、リード全体の情報から作成するすべてのk-merを使うのではなく、Marker というランダムに選んだk-merで部分的な配列情報をオーバーラップの有り無しに判定するために利用しています。

つまりリードにそのk-merを照らし合わせると、以下のような感じです。

f:id:marimiya_analysis:20191223222831p:plain

Computational Methodsより抜粋

ちょっとびっくりですが、ある程度の長さを想定しているのであれば、確かにすべての配列情報を使わず、一部だけ使ってオーバーラップを見つけるのは、なるほどなと思ったりしました。図では3塩基ですが、Defaultは10塩基で、これはパラメータで変更可能です。

このRLEとMarkerをリードに当てはめた図をみてみましょう。これを見るとかなりすっきりですね。黒い塩基配列の上に表示されているのがRLEですね。

f:id:marimiya_analysis:20191223222933p:plain

Computational Methodsより抜粋

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

ナノポア解析ワークフロー

鶴岡の学会で発表した内容が早口でメモできなかったというお声を聴きまして、内容を掲載します。諸般の事情により、スライドをそのままアップすることが出来ないため、ワークフローの部分のみの掲載となる点、ご了承ください。

 

以下、ナノポアのDe novo アセンブリの解析ワークフローと構造変異のワークフローとなります。ご注意いただきたい点は、実際の解析での評価ではなく、いろいろ見聞きした情報を元に集めている点、特にDe novo アセンブリは、扱うデータの難しさ(配列の複雑さや倍数性、得られるデータとリード数)によりプロセスは大きく変わります。そのため、以下のフローはあくまでも、シンプルな方法と理解いただき、個々のデータに合わせて修正が必要となる点もご理解ください。

 

では、まずde novo アセンブリから!スライドのスクショのコピーという点はご勘弁を・・

 

       f:id:marimiya_analysis:20180925001748p:plain

 

まずはQCやフィルタリングですね~ なんと、ONTのデータにはアダプターが残ったりしているので、それはporechopで取り除いてくださいね。ここに記載のあるツールは特別な記載がない限りオープンソースです。nanostatは数値での統計値情報が得られるのが個人的には好きですが、nanoplotを使っている人が多いなと思ったりします。フィルタリングはnaonfilt で。

 

                            f:id:marimiya_analysis:20180925002032p:plain

そしてロングリードのアセンブリツールとしておそらく一番有名どころなのは、canu。時間がかかる、というような声も聞きますが、そのあたりはエラー補正のツール選択にもよるらしいですが、精度では非常に高いものがでるようです。

        f:id:marimiya_analysis:20180925002203p:plain

エラー補正単独、特にDe novo アセンブリの前に使える補正(ナノポア自身のデータで補正する)では、nanopolishがあげられるかと思います。naonpolishは、補正以外にも1塩基の変異検出や、最近はRNAのPolyAの長さの検出など、様々なオプションというかツールが備わっています。いつか詳細を記載できればと思います。作成しているのは、Jared Simpson というABySSを作った人ですね。これで「おー!」と思う人ってどれぐらい居るんだろうなと思ったりしますが、とりあえず次にいきます。

 

 

 f:id:marimiya_analysis:20180925002558p:plain

次は、De novo アセンブリ、ポリッシング、それぞれのパートを行うツール群です。原理などがもちろん違うわけですが、とりあえずこのポストでは、ツール名の列記で失礼します。あとでリンクの調べ先を記載しますね。

このブログを書かなくてはと思っている間に新しくアップデートされたのが、MaSuRCAアセンブリ!なんとヒトゲノムのデータで、論文では6.4MbのNG50だったものを8.4Mbまでに伸長させるなど、Contiguinityが改良されているようです。スライドに入れられず、残念!
masurca.blogspot.com

ポリッシングのステップは、主にHybrid(イルミナリード)を組み合わせることが多いかと思います。ナノポアのデータだけで、高いコンセンサス精度を得ようとすると、解析が大変になったり、カバレッジをあげないといけないため、コスト面からすでに手持ちのイルミナと組み合わせたり、ちょっとだけイルミナを足したりすることが多いようです。

f:id:marimiya_analysis:20180925002750p:plain

そして、De novo アセンブリとポリッシングを一体にしたツールも存在します。UnicyclerとPomoxisですね。Unicyclerは環状化(環状ゲノムを環状に出来る場合、環状にする操作)もできるそうで(名前からすぐにピンとくればよかったのですが)、それを気に入っている方もいるようです。

canuは環状化のオプションがあります(が、つないではくれないそうです(荒川さん、ご指摘ありがとうございます!))。miniasm はありません 訂正:デフォルトで環状化できる場合してくれるようです。オプションになかったのでないのかと思っていました(鈴木さん、ご指摘ありがとうございます)。

PomoxisはONT社が作成したパイプラインでONT社のGithubで公開されています。ONT社のツールはほとんどお魚の名前なんですよね。ほとんどが・・・

そして、以下が構造変異(SV)検出ツールです。QCなどは前述と同じです。

      f:id:marimiya_analysis:20180925003237p:plain

マッピングツールは本当にいろいろありますが、その後の解析ツールとの相性やスピード、パラメータの選択が自分の解析データにあっているか(BAMファイルで必要な項目がちゃんと記載されているかetc)、などが選ぶところでしょうか。マップ率の高い、低いはパラメータでどうとでも変わってしまうので、マップされるべきものがちゃんとマップあされているかなど、見れるといいですね。

     f:id:marimiya_analysis:20180925003439p:plain

そして、SV検出。発表の時は、NanoSVを書き忘れていました、ごめんなさい、、とりあえず追記しています。

 

そしてマッピングとSV検出をまとめているツールがPickyです。

   f:id:marimiya_analysis:20180925003532p:plain

これらのツールは、おそらくググるとすぐに出てくると思いますが、これ以外にもいろいろと新しく出てきているものなど、オーストラリア国立大学のベンジャミンさんが以下にまとめてくださっているので、是非定期的にのぞいてみてください。

ナノポアツールまとめ

この記事を書かなくては…と思っている間にも、またporeTally というツールがリリースされ、以下のパイプラインが実装されているようです・・・明日の新幹線で論文を読んでみます。。これだけやれるといいですね。。Flyeはde Bruijnグラフではなく、アライメントベースのA Bruijnグラフを使ったアセンブリですね。あまり使ってる人を聞いた事がないのですが、以下に含まれているようですね。

f:id:marimiya_analysis:20180925004527p:plain

 

www.biorxiv.org

 

と、ナノポアのツールは次から次へと出てくるので、大変ですが、こちらが何かの参考になれば幸いです。

NGMLR とSniffles

すっごく久しぶりになりますが、ロングリード用のマッピングツール(NGMLR)とその構造変異のツール(sniffles)に関する論文を少し細かく読んだので、ブログに書いておくことにしました。

Accurate detection of complex structural variations using single molecule sequencing

Nature Methodsに2018年4月に掲載されています(OAではありません)。

https://www.nature.com/articles/s41592-018-0001-7

プレプリントで同様の内容が読めます。

https://www.biorxiv.org/content/early/2017/07/28/169557

サプリメントがとても充実しているのと、方法の詳細は、サプリメントに記載があるため、詳細はサプリメントを中心に記載します。

https://static-content.springer.com/esm/art%3A10.1038%2Fs41592-018-0001-7/MediaObjects/41592_2018_1_MOESM1_ESM.pdf

 

この論文は、 Baylor College of MedicineのFritz Sedlazeck とUniversity of ViennaのPhillip Reschenederが開発しており、Fritz とPhillipは昔University of Viennaで研究をしていたようで、NGMLRを出す前にNexGenMap として論文をだしています。NGMLRはその発展版のような位置づけになるようです。

 

それぞれのツールの特徴を大きくまとめると以下。あ、このブログはおそらく長くなるので、ご容赦ください。

NGMLR

  • Long Read 用に開発されたMapping ツール
  • NextGenMap というMappingツールのLong Read 版
  • NGMLRでは、アライメントのGap コストにConvex gap を利用している

Sniffles

  • Deletion, Duplications, Insertions, Inversions, Translocations, Nested Events といった様々なタイプ(論文ではAll typesと記載)の構造変異の検出が可能
  • Sniffles では、リード1本のアライメントに見つけることの出来る変異とSplit reads と呼ばれる場合の変異も検出ができる。
  • Snifflesの新しい点は、新しいSVスコアリング方法。これによりエラーが多い、PacBioやONTのリードを使った構造変異検出に利用できるとしている。

ではNGMLRから

f:id:marimiya_analysis:20180915230355p:plain

Fig. 1の図に大まかな2つのツールの流れが出ています。左がNGMLR、右がSnifflesですが、右のSnifflesは、さらに詳細な図が後述されるので、左図で、まずはNGMLRの説明をします。

Subsegment alignment とLMPの検出

まず、最初のステップがSubsegment alignment。図ではSub-readとなっていますが、のちに修正されているようなので、Subsegmentで行きます。

  • Read を256bpにカット(Subsegment)し、Reference 配列へアライメント。
  • Referenceスコアへ高いアライメントスコアでマップされた領域をSubsegementのマッピング領域候補とする。
  • アライメントスコアが75%以下のものは捨てる
  • 1000個以上マップされる場合(Repetitive領域など)は、
    そのSubreadは捨てる。
  • 残った候補から、Linear Mapping Pair (LMP) を定義.
  • ゲノムに対して同方向で近い位置でマッピングされているペアを見つける
  • このステップでは、SVを含んだ領域ではなく、まず精度よくマッピングされているサブセグメントの領域を見つけるというステージ。

 サプリメントのFig 1.1が良いので入れておきます。

f:id:marimiya_analysis:20180915231329p:plain

 

a)では、オーバーラップがないようにカットされ、 その候補が スコア/ポジション として記載されています。そうすると、1本のリードが複数個所にマップされ、そのうちいくつかポジションが近いものをペアとすることができます。b)の黄色の個所と青の個所です。これらはLMPとされます。そしてさらに、d)の2つのサブセグメントは、距離は離れてはいるものの離れた距離がReference上もリード上も同じ距離の場合、そのセグメントを結合してLMPとしています。

 

Convex Gap Alignment

そしてスコアの計算はConvex Gap スコア。ここでアライメントスコアのおさらいです(参照する配列と検索をかける配列の類似度を測るスコア)。

 

一般的なアライメントのGapスコア。上の配列を参照配列、下を検索をかけるクエリ配列とします。Wikipediaを参考にして記載しています ( Gap penalty - Wikipedia)。

Linear

      f:id:marimiya_analysis:20180915232017p:plain

  • Match score = 1, Gap cost = -1 の場合、上記アライメントは、7-(1)*3 = 4 となる

Affine

  • Gap open cost, Gap extension cost を採用してスコアを計算。Gap open cost = -3, Gap extension = -1とすると前述のアライメントの場合は、7+(-3)+(-1)*3 = 1 となる

Convex

Affineと同様にGap Open cost を設定するが、Extension に定数×自然対数(Gapの長さ)というような設定をし、Gapが伸びた際のコストをマイルドにする。この論文では、もう少し違った式でConvex Gap を実装しています。

 

 Gapスコアに関しては、構造変異が起きているような箇所は、小さくばらばらにGapが起こるのではなく、ある程度まとまったGap として観察されるという事を前提にGapのコスト計算がされています。

NGMLRでは、Convex を以下のように定義しています。

f:id:marimiya_analysis:20180915233037p:plain

 g_o は、Gap Open コストでGap を開くときのコスト、g_e はExtension コストで、gapを伸長するときにかかるコスト、g_MはGapペナルティが大きく0以下にならないようにするMinimumの値、そしておそらくこのDecay を入れているのが新しいのかなと思いますが、g_D として定義しているGap Decay コスト。

少し細かくなりましたが、サプリメントのFig 1.2 に良い図があります。

  f:id:marimiya_analysis:20180915233531p:plain

 著者らが実装したNGMLRのConvexでは、大きなGapを受け入れやすくなっています。Fig 1.2 でも、Affineコストでは、大きなGapと複数のGapが同じスコアになるのに対して(a)、採用したConvexでは、大きなGapがより大きいスコアとなっていることからもわかります(b)。

 

Small InDel Detection

次のステップは、小さなInDelの検出です。

小さなSmall Inversionの検出のために以下を行っています:

  • LMPの3bpごとのWindowsでアライメントスコアを計算
  • スコアがある閾値(65%)より小さいものがある場合、その領域を抜き出す
  • Reverse Compliment を作成し、それがReferenceにマッチし、アライメントスコアが上がるか調べる

これもサプリメントに良い図があるので、記載しておきます。

 f:id:marimiya_analysis:20180915234127p:plain

 

Alignment Selection

NGMLRで最後のステップは、Alignment Selection です。LMPでは、1つのリードが複数個所に分断されてマップされる個所を見つけていますが、ゲノムビューアで閲覧できるように、1つのリードをどのポジションにマップさせるか、決めなくてはいけません。そのために以下のようにしてマップさせるポジションとマッピングクオリティを決めています。

  • Linear Mapping Pair (LMP)の組み合わせで、オーバーラップなく、もっともスコアが高くなる組み合わせを選択
  • Mapping Quality は、SubreadsのMapping Qualityの平均

 

そしてようやく次がSnifflesです。

 

 これもサプリメントに詳細なプロセス図があるので、まずはそちらを掲載します。

  f:id:marimiya_analysis:20180915234634p:plain

Prameter Estimation 

最初のステップは、パラメータの推定です。Snifflesでは、3つのパラメータを推定します。このパラメータ推定は、ランダムに選択した10,000本のインプットリードから行います。SVは稀に観察される事象であるという仮定で、この10,000本がSVにオーバーラップせずにベースの値として利用できることを想定しています。

  1. アライメント中のIndelやMismatchによる違い。
  2. 100塩基中にみられるのミスマッチとIndelの個数の95%分位点。
  3. 各リードがマップされたBestと2番目に良かった場所のアライメントスコアの比

 

Detection of SV

Filtering : 以下の条件に当てはまるものは捨てる

  • Mapping Quality < 20
  • Bestと2番目のアライメントスコアの比が2以下 (つまり似ている箇所のアライメントがある)
  • 7か所以上にSplitしたリードがマップされる場合

そして次のステップでSVの領域を検出します。

  1. まず最初に小さな(<1kb)Insertion, Deletion とNoisyな領域(ミスマッチが多く、1-5bpのInDelが見つかる領域)の検出。
  2. 1つのアライメントでは説明できない大きなGap(Split reads)によるSVを検出

    Split reads の情報は、マッピング結果のSAM/BAMファイルのSAタグを利用

  3. ステップ1と2で検出されたSVを支持するリードは、それらが同じSVに由来するか分類
  4. 疑わしいSVは捨てる

サプリメントのFig 2.4 にSnifflesが検出できるSVと、Split reads のマップのされ方が記載されています。

f:id:marimiya_analysis:20180916001712p:plain

SVをそれぞれ上記分類に分けるプロセスは以下。

1.異なるアライメントが同じ染色体で同方向の場合、かつ長さが指定したSVの長さ(default = 30bp)っよりも長い場合、Insertion とコール。Deletion も同様にコール。Duplicationは、Reference上の距離がユーザーが指定するSVの距離よりも大きく、またリード側の距離がユーザーが指定するSVの距離より小さい場合にDuplicationとコール。

2.2つのセグメントが同じ染色体上にあるが、異なるStrandにある場合、Inversion とコール。2つのセグメントが少なくとも200bpオーバーラップし、長さが40%以上すくないものはInversion Duplicatin (U-Turn)とコール

3.2つのセグメントが異なる染色体にある場合は、Translocation とコール

Break pointの検出

  • 異なるRead が支持する同じSVは、Breakpointが一致またはそのばらつきは小さな範囲に収まるはずと仮定。
  • Sniffles では、このばらつきは正規分布に従い、その標準偏差を以下に定義

f:id:marimiya_analysis:20180916002633p:plain

 

  • 1st、4th quantile (5%, 95%)のばらつきを持つリードを除去することで、より正確なBreak point を検出できるとしている。

   f:id:marimiya_analysis:20180916002522p:plain

と、こんな感じが、SVの検出のところまでです。Snifflesはこのほか、ジェノタイプの検出やフェージングも行えるようになっています。

このツールを使った結果は、シミュレーションと実データでの検出が行わわれていますが、その辺りは論文に分かりやすく記載されていると思います。

補足としてパラメータの記載がサプリメントに会ったので、入れておきます。

NGMLRのパラメータ

f:id:marimiya_analysis:20180916003044p:plain

Snifflesのパラメータ

f:id:marimiya_analysis:20180916003210p:plain

以上になります。また時々こういうブログが書けたらいいなと思います。お気づきの点がありましたら、ご連絡ください。

アンガーマネージメント

このブログには基本データ解析に関する話を書こうと思っていたのですが、自分の備忘録的に久しぶりに投稿します。

今日は久しぶりにfuriousな事案が発生しました。furiousは日本語だと怒り狂うみたいな感じですが、なんというか、怒った(angry)のレベルではない、furiousという感じ。

私も最初は今ひとつピンと来ませんでしたが、アメリカで仕事を数ヶ月していたときに友達が、今、僕はfuriousなんだよ、分かる?angryじゃなくてfuriousなんだよ、と彼の仕事に対して文句を言っていた時、ああ、このぐらい怒ってるときにfuriousを使うんだなーと思ったのを覚えています。そう、あのレベル。

怒った理由を書くのも大人げないのと、平常心に戻って考えれば、それほど怒らなくても良かった案件だなとも思ったりするので、詳細はスルーしますが、久しぶりにこのレベルの怒りが発生して、平常心に戻るまでにいくつかのステップがあったので、ここに書いておこうと思います。

外国だと、Anger Management はマネージメントのレベルになった人が受けさせられるようなクラスだったりするのですが、日本ではあまり聞きません。外国では部下を怒鳴る、という事はあまりなく、そういう人は、Anger Managementが出来てない、Education level の低い人という烙印を押されるような場合もあるので、頭ごなしに怒鳴ったりというのは、ほとんど見かけません。まぁ居るとは思いますが、どちらかというと良い評価にはなりません。なので、みんなAnger Managementが出来ている感じです。 私はAnger Managementのクラスを受けた事がないので、自分なりにどうやって今日の怒りをコントロールしたか、書いてみました。

1.まず、何に対して怒ったのか、自分を分析

これ、何かの本に書いてあったと思うのですが、自分の感情を分析する立場に立つ事で、客観的に物事を見る事が出来、冷静になれると。確かに少し落ち着くために自分の感情を分析してみました。なんでも分析大好きなので、自分が怒っている理由を考えると、それは自分が思ったようになってないから。そしてそれは、自分の中では、約束したら守るもの、という定理のようなものが守られなかったから、怒っているのだなと。

2.場所を変える

私は在宅勤務のため、あまり話す相手もおらず、これは良くないと、とりあえず晩御飯の買い物も兼ねてスーパーへ買い物へ行きました。少し外に出ると気分も落ち着きました。

3.ぜんぜん関係ない人と、くだらない話をする

この事案に全く絡んでいない人と、Lineでくだらないスタンプを送りあったりして、今度会う予定を決めていたりしたら、面白いスタンプを見てつい笑ってしまい、笑えた事が嬉しくなりました。

4.別の事で頭を集中させる

これ、よく言われることですが、実際に実行するのは難しいと思っていました。今回はスーパーで「あれ、切れてた調味料なんだっけ?なんで切れてる事に気づいたんだっけ???」と昨日や数日前のおかずを思い出していたら、自然と頭はその事だけ考えるようになり、「あー、お酢がないんだ!」と思い出すと、それはそれで嬉しくなりました。これはいいなという事で、他にも何かあったよな、と家中の洗剤やら、お風呂場で使うものを順番に思い出しながら、買ったほうが良かったもの、ストックが切れてるものを考えると、自然と怒った内容は考えなくてすむようになりました。これは大きな発見でした。次回もこの手で頭を怒りの事案から頭を遠ざけようと思います。

5.プランBまたはCを考える

スーパーの中をぐるぐる歩きながら、自分のゴールは怒って相手に仕返しするとか意地悪するとかじゃなくて、今やってる仕事を成功させることなんだから、そのために何が今足りなくなって、どうやって埋めるのかというのを考えていたら、いくつかアイデアが出てきました。そして、もしかすると、代わりのプランの方がいいかもしれないとも思うようになりました。次のプランが出てきて、それの成果の可能性が見えると、怒りなんてどこかに飛んでいく感じでした。

6.好きな音楽を聴く

私はEDMやハウスのようなハードな音楽と、クラシックが好きです。今回はEDM系で、いくつか音楽を聴きながらスーパーをぐるぐるまわって、頭をトランス系に切り替えると、なんでも出来そうな気がしてきました。EDMやハウスミュージックのいいところは、こういうトランス状態が健全に作り出せるところだと個人的には思っています。

そして家に戻って、ご飯を作ってから、あらためてAnger Managementでネットを見ていると、Mayo ClinicからAnger Management: 10 tips to tame your templer というのが出ていて、そのなかで I から始まるセンテンスでコミュニケーションするというのが6番目にありました。 For example, say, "I'm upset that you left the table without offering to help with the dishes," instead of, "You never do any housework." と、こんな感じ。なるほどねぇ~ と思い、早速さっき怒っていた相手にメールで I was upset when I head about the news but let's work together for an alternative plan. みたいな感じでメールしたら、相手もさすがに悪いと思ったのかSkypeで電話がかかってきて、まぁ一件落着。一緒になんとか代替案を考えるという事で、雨降って地固まる、みたいな感じとなりました。

社内の相手ですが、多くの場合悪意を持って約束を守らないというわけではなく、事情があっての事ですから、そんなに怒るよりは、貸しを作ったと思った方がいいですね。 という事で、比較的うまく怒りをmanageできたと思います。

平常心に戻るまでには、今回50分かかりました。次回は30分程度で気持ちを切り替えられるようになりたいものです。

EdgeRと一般化線形モデルまわり

ちょっと仕事も含めて、EdgeRの一般化線形モデルをじっくり見る必要があったので、そのあたりのメモ。

EdgeRは以前見たときは、1因子はExact Test をしていて、多因子の時が一般化線形モデルという話だと思ったけれど、最近は1因子でも一般化線形モデルという事らしい。DESeqもそうなっているので、まぁそういうことなのだろう。

一般化線形モデルはリンク関数を使って、応答変数を変換し、線形予測子でモデルを作るようにしたもので、応答変数が正規分布をしてなくてもいいというところで応用の幅が広い(ただし指数分布族である必要はある)。

で、いきなりEdgeRでのモデルを見てみる。マニュアル18ページに出てる。ちなみに2015年10月8日にReviseされたバージョン。

 \begin{align} \log\mu_{gi} = {\bf x}^T_i \beta_g + \log N_i \end{align}

ここで \mu_{gi} g番目の遺伝子、 i番目のサンプルの観測値の期待値、{\bf x}^T_i はデザイン行列、 \beta_gは係数。そして \log N_iはサンプル i番目のライブラリサイズ、つまり総リード数。ここで、ふと思ったわけです、この \log N_iってなんでそこに居るの、と。普通教科書に出てくる一般化線形モデルは、以下

 \displaystyle \log\mu_{gi} = {\bf x}^T_i \beta_g

みたいな感じ。で、とりあえず分かりやすくするために左辺をもとの値に戻してみるかと、式変形してみると


\begin{align}
\log\mu_{gi}&={\bf x}^T_i \beta_g + \log N_i \\\
\exp(\log\mu_{gi})&=\exp({\bf x}^T_i \beta_g + \log N_i) \\\
\mu_{gi}&=\exp({\bf x}^T_i \beta_g)N_i \\\
\end{align}

で、マニュアルの14ページあたりに、以下の記載があるわけです。


\begin{align}
E(y_{gi})=\mu_{gi}=N_i \pi_{gi}
\end{align}

ここで、 y_{gi}が実際の観測されたリード数、その期待値が、総リード数に割合 \pi_{gi}をかけたものであらわされていて( \Sigma^G_{g=1}\pi_{gi}=1)なんだ、そういうことだったんですか、Robinsonさんと著者の人を思っていました。まぁ自明だから説明省くよ、って事なのかもしれませんが、なるほど、総リード数の割合のところを推定するようなモデルなのだなーとなるほどなるほどと思っていました。
バイオインフォのアルゴリズムは、当たり前といえば当たり前なのでしょうが、教科書どおりとはならないので、有る意味解読的に読まないといけませんが、おや?と思ったところが分かると、パズルが解けたような面白さですね。と、教科書には載ってないと書いて、よくよくいろいろな本を見てみると、一般化線形モデルの教科書的な久保先生のみどり本には、47ページ脚注に説明があり、138ページで同様な形のGLMが紹介されていました。あはは。。。GLM奥深し。。大学院で習った記憶はあるものの、あまり専門ではないのですが、これを機に勉強しなおそうかなとか思ったり。

まぁSlow and Steady でいきませう。

肌に付着した菌叢からの個人特定

 

これはアドベントカレンダーの企画の記事です。これを機にブログなどはじめてみようかなと。

仕事柄、比較的よく論文を読み、ほほーと思っているので、今年一番面白かった論文というのを決めるのは難しいのですが、クライムサスペンスドラマ好きな私にとっては、肌に付着した微生物の菌叢から個人を特定しようとした論文が、個人的興味と重なり印象に残りました。

論文はコロラド大学の環境系の研究室からPNASに発表された論文、

Forensic identification using skin bacterial communities.です。

 

では見ていきましょう。

 

まず最初のデータは、キーボードに付着した菌叢と指先の菌叢の比較。20歳から35歳、過去6ヶ月の間に抗生物質をとっていない健康な人。3人のうち二人は同じオフィスで勤務。キーボードは個人が専用で利用しているもので、サンプル採取30分前は誰も触れていないもの。16S rRNAを454でシークエンスしています。

 

 

論文中Fig. 1は以下のとおり。PCoAで距離はAがUnweighted UniFrac、BがWeighted UniFrac。

 

f:id:marimiya_analysis:20151219214119j:plain 

 

 

 おおお、結構似てる~ とちょっとびっくりしました。

 

そして次は、どのぐらい肌の菌叢の状態が保存されるのか。サンプル保存期間(DNA抽出までの期間)を3日と14日、そしてそれぞれ20度と-20度で保存。サンプルは健康な大人2人のわきの下から菌叢の採取。肌が選ばれているのは、2009年に発表された論文で肌が腸内や他の場所に比べて個人間の多様性が高い場所であるからで、さらにわきの下は、ある程度の菌の量が期待できるからという事でした。 たしかにわきの下にはいっぱい居そうですね。。。

 

そして結果が以下、論文のFig 2. こちらもPCoAでAがUnweighted, BがWeighted のUniFracが距離。日が変わってもPCoAは割と綺麗に分かれている。CはGenus (gen), Family (Fam), Order レベルでのrelative abundance. レベルにもよりますが、個人内の変化は個人間の変化よりも小さいのだなという印象。 

 

f:id:marimiya_analysis:20151219222431j:plain

 

 

これまでの結果は数人の結果なので、次のデータでは9人からその人が使っているコンピュータマウス表面から菌叢採取、そしてマウスを使っていた手のひらからも菌叢採取。マウス所有者の手のひらとマウス表面の距離を計算し、さらにその9人のマウスとデータベースに登録されている270人分のデータを比較して距離を計算。距離は、Unweighted UniFracとweighted UniFrac。赤がUnweighted、青がWeighted。横軸がサンプルの9人。色がついてないシンボルは、マウス所有者の手のひらと、その人が使っていたマウス。色が付いているシンボルは、マウス表面と270人の手のひらのデータ。

f:id:marimiya_analysis:20151219230531j:plain

つまりマウス所有者の手のひらと、マウス表面の菌叢は、マウスとデータベース270人分の類似度よりも似ているという事を示している。F8、M2のWeightedはちょっと似ているけれど、エラーバーにかかっていないので、ぎりぎり区別はできそうな雰囲気。

 

この論文は2010年の論文で、たまたま今年読んで個人的にほほーとなった論文でした。この論文を引用している論文は2015年12月19日段階で236本(Google Scholarしらべ)。いろいろ面白そうなものがあったので、またこのブログで紹介したいなと思います。

 

こういったメタゲノム解析以外にも、NGSを使った個人鑑定に関連するものとしては、STR(Short Tandem Repeat)を見つけるものや、ミトコンドリアのHyper-variable Region を読んで個人を特定しようとするものなどもNGSで行われているので、機会を見つけて書いていければと思います。

 

と、最新の知見でなくて恐縮ですが、こういった技術がいずれは法医学の現場で使われるのかな~など思いをはせていました。真実は小説より奇なり。論文は小説よりも面白いなと思うことがあるので、冬休みにかけていろいろ読んでみたいですね。

 

ではまた。