krakenに代表されるメタゲノムの分類学的プロファイリングツールは、fastqのそれぞれのReadに対してダイレクトに分類学的分類を行う。そのために、kraken1ではJellyfishを使ってリファレンスゲノムからk-merが取り出され、ゲノムの分類学的情報と共にデータベース化されて利用される。データベース作成後、kraken実行時は論文の図1(引用)にあるように各シークエンシングリードから複数のk-mer配列が取り出され、メモリにロードされた巨大なデータベースとk-mer配列の完全マッチが探索される。それぞれのk-merの完全マッチのLCAが決定され、そのヒット結果を総合して1つのリードの分類学的情報が得られる。
kraken2の解析例。単離シアノバクテリアを読んだfastqのたった1つの配列だけを分析した。一本のリード(2列目)はrootから種レベルまでアサインされている。
kraken2ではk-merの代替としてminimizerを使うことでメモリ使用量が大きく削減されている(それもでもスタンダードDBでは数十GB使い、GTDBなど巨大なDBでは300GBほど使う(古いバージョンにも関わらず))。これをインメモリデータベースとして使うこともあって、krakenは超高速に動作し、例えば細菌ゲノム数万個からなるデータベースであっても、マッピングベースの手法(*1)よりも非常に速い速度で分類学的プロファイリングを行うことができる。例えば浅く読まれたメタゲノム(shallow metagenomicsと呼ばれる)なら巨大なデータベースを使っているにもかかわらず、1時間で数百〜千サンプルほどのプロファイリングを行うことができる(紹介)。これは驚くべきスピードと言える。一方で明確な制限もあり、krakenはk-merの完全マッチに基づくため、リードのエラー率が高いと、シークエンシングリードのどの部分から取り出されたk-merも完全マッチがなくなってしまう。環境サンプルであればリファレンスとの完全マッチがなくなってしまう事が多くなるため、エラー率が高いと結果はさらに悪化してしまう。例えばilluminaの高精度なリードを使っていても、メタゲノムサンプルであればkrakenの分析で数十%がアンマッチとなることはそこまで珍しくはない。また、k-merの完全マッチは複数の分類群にまたがって起こりやすいため、後のKrakenUniq(リコールと精度を改善するアルゴリズム)(論文)やbracken(krakenの結果のベイズ推定)(論文)などの派生版が生まれている。特にkraken1/2の結果をbrackenと組み合わせることはよく行われている。それでも環境サンプルは、研究で使用されるDBが地球の全ての生命のゲノムのスーパーセットではなく、ある時点でサンプリングされたゲノムのサブセットでしかないため、種レベルの推定には誤差が大きくなってしまう。誤判定を減らすため、結果を属より上の分類まで刈り上げて評価することもよく行われている。
このkraken2はロングリードでも動作はするが、どのくらい精度が出るかはシークエンシングリードの精度に依存していて、よく分からないことが多い。例えば100塩基の配列にランダムに10%のエラーがあれば(Phredクオリティスコア=10)、1-bpずつずらしながら30mer前後のk-mer配列を取り出した時のDBとの完全マッチはほぼゼロになると考えられる。シークエンシングエラーによって誤判定(偶然のマッチ)がどのくらい出るかも予想が難しい。評価する研究も存在するが、評価のためにground truthが分かる模擬メタゲノム(単離された配列特徴量が異なる菌を指定の量混ぜ合わせ、そこからDNAを取る)が使われていて、"外"の環境サンプルへの適用性は不明な点が多い(*2)。
k-merベースの手法とは異なり、mOTUs3は、リードを固有のマーカーデータベースにマッピングする。これはmetagenomic OTUのシングルコピー系統マーカー遺伝子で構成されるデータベースであり、このデータベースへのリードのマッピングに基づき、mOTUsは分類群の相対的な存在量を算出する(カバレッジスコアに基づく)。ハウスキーピングなマーカー遺伝子を使うため、メタゲノムだけでなくメタトランスクリプトームの定量も現実的に行うことが可能とされる。mOTUsは最近ロングリード向けのオプションが追加された。これはロングリードを複数のショートリードに分割して定量するだけの単純なものだが、実際にどのように動作するかどうかは気になるところである。mOTUsのロングリード向けコマンドを確認しておく。
Profile long reads
https://github.com/motu-tool/mOTUs/wiki/Profile-long-reads
インストール手順は以前紹介しました。
> motus prep_long
$ motus prep_long
Usage: motus prep_long -i <long_read_file> -o <converted_file> [options]
Input options:
-i FILE long read file to convert, can be fasta(.gz) or fastq(.gz)
Output options:
-o FILE converted file (gzipped), ready to be used by motus profile
-no_gz save the output file without gzipping it
Algorithm options:
-sl INT splitting length for the long reads [300]
-ml INT minimum read length, shorter are discarded [50]
-v INT verbosity level: 1=error, 2=warning, 3=message, 4+=debugging [3]
(mOTUs 3.0.3以降が必要)
実行方法
1、各ロングリードを複数のショートリードセグメントに分割する。
motus prep_long -i long_reads.fasta -o converted_long_reads.fasta.gz
converted_long_reads.fasta.gzができる。
2、mOTUsのラン
motus profile -s converted_long_reads.fasta.gz -t 20 -n sample1 -o taxonomy_profile_sample1.txt
出力例(単離E.coliのONTリード使用)
taxonomy_profile_sample1.txt
問題なく動作している。
引用
Cultivation-independent genomes greatly expand taxonomic-profiling capabilities of mOTUs across various environments
Hans-Joachim Ruscheweyh, Alessio Milanese, Lucas Paoli, Nicolai Karcher, Quentin Clayssen, Marisa Isabell Keller, Jakob Wirbel, Peer Bork, Daniel R. Mende, Georg Zeller & Shinichi Sunagawa
Microbiome volume 10, Article number: 212 (2022)
関連
*1
例えば細菌ゲノム50万個単純にを連結するとヒトゲノムの500倍の大きさになってしまう(ゲノムサイズ2.9MB換算)。マッピングベースの手法でそれほど大きいリファレンスは想定されておらず、例えばHTSlibの32-bit制限で4GB以上のゲノムでトラブルがあったり(samtools #1117)、大きな配列のindexingでエラーが起きるなどする(参考)。
*2
データベースにない種では挙動の予測が難しい場合がある。例えばあるDBをつかった時に属レベルのカウントがゼロだったものが、DBを変えると、同じfastqの分析でも全体の1%くらいその属へのカウントが出ることがある。DBにその属のゲノムが1つでも含まれば確率的にマッチするため、カウントは減ってもゼロにならない気がするが、必ずしもそうはならないということ。曖昧さを許容するマッピングベースの手法ではこのようなピーキーな結果にはならないと思うが、明確な原因は不明。
参考にした論文