kraken2のレポートをkrona plotで視覚化する - macでインフォマティクス

2024/02/14 誤字修正

メタゲノムデータ解析レシピ（ISBN 978-4-7581-2255-9）3章のWEB年度更新で、kraken2のunclassifiledの割合には注意しましょうという説明をしました。その中で、unclassifiledがkrona plotには反映されないと書いたのですが、これはKrakenの出力をBrackenに読み込み、Brackenで推定したレポートをkrona plotにプロットした時の話でした。Krakenの生のレポートをそのまま使った時にはunclassifiledもkronaでプロットされます。ここで修正させていただきます。

ではなぜBrackenをKrakenと共に使うのでしょうか。今日はこれについて書いてみます。まずkrakenのアルゴリズムについて簡潔に説明します。年度更新でも書いてますが、krakenはk-merの完全一致による超高速なシークエンシングリードの分類学的分類のプロファイラーです（Wood DE et al, 2014）。krakenは各リード配列内の全てのk-merを取り出し、それらを全ゲノムのLCA分類群にマッピングしてk-merの数で重み付けします。ルートからリーフ（種）まで、パスに沿ったすべての重みの合計を計算することによってノードの得点を計算します。最大得点のパスが正しい分類とみなし（通常は一意）、リードそれぞれに対して対応する分類学的ラベルを割り当てます。krakenは塩基またはアミノ酸（kraken2のみ）の完全一致を利用しているため、同種であっても、データベースから少し距離のある株のリードの場合、どの分類群とも一致せず、unclassifiedなリードが多数になってしまう可能性があります。数万以上のリファレンスゲノムにスケールできるものの、データベースが充実していないと正しい分類が難しくなるという欠点があるわけです。2019年にメモリ使用量を数分の一と大幅に削減したkraken2の論文も出ていますが、DB依存性が高い点は変わらずです（Wood DE et al, 2019）。

krakenのアルゴリズムのもう１つの問題点として、種間や属間で高度に保存されたゲノム領域があると、最大得点のパスが種や属レベルで2つ以上出てきてしまい、krakenは種や属のレベルまでユニークにリードを分類することができなくなる点が挙げられます。そのような時、krakenはlast common ancestor (LCA)までリードに分類をアサインします。種レベルでユニークな分類が不可能なら、ルートから属まで、属も無理ならルートから科まで、ということです。しかし、このようなリードが多数になると、種レベルや属レベルの定量をしたい時、誤差が無視できないほど大きくなります。2017年に発表されたBracken（Jennifer Lu et al, 2017）はこの点に注目した実装で、リードを確率的に再分布させることによってデータ中の種の存在量を推定します。再分布させたいランクはユーザーが指定します。例えば種を指定したとしたら、Brackenは、論文の式１にあるベイズ推定によって、種レベルより上のノードにアサインされたリードを種ノードの下に分散し、株レベルでアサインされたリードをそれらの親の種上に再分散して確率的なカウント値を算出しようとします。詳しくはそれぞれの原著論文をご確認下さい。

Brackenがこのような性質を持つため、Brackenをkrakenのあとに実行することで、種や属レベル定量結果を改善できます。

1、Kraken2

この流れを実際に確認してみましょう。ここではKraken v2を使います。Kraken2とBrackenは前もってインストールされ、データベースも既にローカル環境にダウンロード・ビルドされているとします（詳しくは年度更新を参照して下さい）。

#way1 single sample
kraken2 --db <path>/<to>/kraken2_DB/ --threads 20 --gzip-compressed --paired --report kraken2_report.txt --use-names --classified-out seqs#.fq R1.fq.gz R2.fq.gz > output.txt

#way2 forループ
mkdir report
for file in `\find *R1.fq.gz -maxdepth 1 -type f`; do
 fastq1=$file
 fastq2=${file%_R1.fq.gz}_R2.fq.gz
 kraken2 --db <path>/<to>/kraken2_DB/ --threads 20 --gzip-compressed --paired --report report/${file%_R1.fq.gz}_kraken2_report.txt --use-names --classified-out seqs#.fq $fastq1 $fastq2 > report/${file%_R1.fq.gz}_output.txt
 rm seqs_1.fq seqs_2.fq
done

#way3 たくさんfastqがあるなら、DB読み込みを1度で実行できるフォークを使用(解説)
mkdir report
kraken2 --use-names --report-zero-counts --threads 20 -db <path>/<to>/kraken2_DB/ --output report/ --report report.txt *fastq

-reportで指定するレポートテキスト：kraken2_report.txtの先頭が下になります。

unclassifiedがまずプリントされて、2行目以降にルートから種までの全リードの分類結果の％がプリントされます。レポートは、検出された分類ごとに繰り返しこの結果をプリントします。

２，Bracken

続いてBrackenを使い、指定した分類階級のカウント数を推定します。ここでは種を　指定します（-l S）。krakenのレポートテキスト（-i kraken2_report.txt）とリード長（-r 150）を指定します。"-d"でkraken DBを指定します。

#way1 single sample
bracken -d <path>/<to>/kraken2_DB/ -i kraken2_report.txt -r 150 -l S -t 5 -o  bracken_report.txt

#way2 forループ
for file in `\find *R1.fq.gz -maxdepth 1 -type f`; do
 bracken -d <path>/<to>/kraken2_DB/ -i ${file%_R1.fq.gz}_kraken2_report.txt -r 150 -l S -t 5 -o ${file%_QT1.fq.gz}_bracken_report.txt
done

通常のレポートなら数秒以内に結果は出力されます。Brackenのレポートテキスト：bracken_report.txtの先頭が下になります。

種レベルを指定したので、先ほどのKrakenレポートのようにルートから種までの全リードの分類がプリントされるのではなく、1行ずつ種レベルの確率的分類がプリントされているのが特徴です。種レベルなので、アサインされていないunclassifiedも消えています。

また、これとは別にkraken形式のレポートファイルも再作成されて保存されます。こちらも未分類のリードはレポートに含まれません。

kraken形式のレポートファイルですが、Brackenのレポートには推定レベル以下のレベルはプリントされず、カウント値が閾値を下回ったレベルは含まれていません。また、Brackenでは扱えない未分類のリードもレポートに含まれていません（特にここがkrakenレポートと異なる部分）。

3、krona

最後にkronaを使ってBrackenの結果を階層的な円グラフ形式で視覚化します。Brackenレポート（画像3枚目のkrakenスタイルのレポートの方）をKraken toolsを使ってkrona互換のテキストファイルに変換し、それからkronaでプロットします。

#step2 brackenからのkraken形式brackenレポートをkrona互換形式に変換
git clone https://github.com/jenniferlu717/KrakenTools.git
python KrakenTools/kreport2krona.py -r kraken2_report.txt -o bracken_convert

#step2
#korna plot。複数レポートがあるならワイルドカード指定
ktImportText *bracken_convert -o krona.html

krona plot

このようにプロットできますが、年度更新で説明した通り、unclassifiedはKrakenの精度を測るための重要な指標になります。参考資料（*1）のようにレポート先頭を別に保存しておくか、Brackenを介さずにkrona plotを描画しておきましょう。つまり、説明の1のステップ、３のステップの順番です。

#step1のkrakenからのレポートをkrona互換形式に変換
git clone https://github.com/jenniferlu717/KrakenTools.git
python KrakenTools/kreport2krona.py -r kraken2_report.txt -o kraken2_convert

#step2
#korna plot。複数レポートがあるならワイルドカード指定
ktImportText *kraken2_convert -o krona.html

krona plot