macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

マルチサンプルに対応したkraken2のフォーク

2023/12/20 追記、12/21 インストール手順修正

 

Kraken 2は、k-merの完全一致を利用したシークエンシングリードの分類学的プロファイリングツールで、メタゲノムやメタアンプリコンの分類や汚染のチエックなどに幅広く使用されている。データベースは自分で作成することもできるが(custom databaseの項を参照)、おそそらくは公開されているスタンダードデータベースがよく使われている(standard databseの項を参照)。このスタンダードデータベースはバクテリア古細菌、ウイルスドメインのRefSeqの全ゲノム、ヒトゲノム、既知のベクターのコレクションなどが含まれている(構築後のサイズは2022年8月でhash.k2dは57GBくらい)。また、より大きなGTDB互換kraken2データベースも有志の手によって公開されている(構築後のhash.k2dのサイズは308GBくらい)。kraken2の実行時は、このデータベースがRAMにロードされる。スタンダードデータベースではこのデータベースのロードにより50GB以上のメモリを必要とする(v1よりかなり抑えられている)。kraken2 GTDB-D.Bだと6倍ほど多くメモリを使う。

kraken2を大量のサンプルに対して実行する場合、このD.Bのロード時間が一番のボトルネックとなる。例えば自分の計算環境ではシーケンシャルリード6500MB/sくらいのエンタープライズ向けPCI接続SSDにスタンダードデータベースを配置してD.Bのロードを行っているが、それでもおよそ30秒かかる。それに対して、ロード後の100万リードの分類は素晴らしく高速で5秒程度で終わる。SATA接続SSDだとデータベースのロードにはより長い時間がかかり、HDDから読み出すとだといつまで経っても終わらないような感覚になる。サンプル数が増えるにつれて、データベースのロード時間が全ランタイムに対して支配的になる(1,000サンプルではD.Bロードに30,000秒(8.3-h)、分析に5,000秒(1.4-h))*1。。データベースロード時間の問題を抜本的に解決するには、複数サンプルを扱う時もデータベースのロードが1回であることが望ましい (すなわち、1,000サンプルではD.Bロードに30秒、分析に5,000秒)。しかし2022年12月現在では本家kraken2は複数サンプル指定に対応しておらず、ワイルドカードでまとめて指定してもスキップされる。

そこで、ここではマルチサンプルに対応したkraken2のフォークを試す。これを使う事で、データベースのロードは一度のみで大量のfastqの分類学的プロファイリングを行うことができる。

 

kraken2 D.B

https://benlangmead.github.io/aws-indexes/k2

 

インストール

kraken2.0.9-betaを試した。

Github

mamba create -n kraken2folk -y
conda activate kraken2folk
git clone https://github.com/daydream-boost/kraken2.git
cd kraken2/
./install_kraken2.sh <path>/<to>/pre-build_kraken_database
#メッセージを読んで、シンボリックリンクを貼る。ここでは仮想環境のbin/に
ln -s <path>/<to>/pre-build_kraken_database/kraken2* ~/mambaforge/envs/kraken2folk/bin/

>kraken2 -v

 

スタンダードD.B

kraken2-build --standard --db kraken2_DB --threads 24

このコマンドを実行するとカレントにkraken2_DBディレクトリが作成され、その中にゲノムやtaxonomyがダウンロードされてデータベースがビルドされる。ハードウェアと回線速度にもよるが、1日近くかかると思っておいた方がよい。

 

 

実行方法

各サンプルのシークエンシングリードを順番に指定する。フルパスでの指定には対応していないので相対パスで指定する。拡張子はfqやfastqなどに対応している。圧縮している場合、解凍してから使用する(本家kraken2はgzipped fastqに対応している)。サンプルが多いときはワイルドカードと組み合わせて指定する。

mkdir out_dir
kraken2 --paired --threads 20 -db kraken_db --confidence 0.5 --output out_dir/ cleandata/sample1_clean_R1.fq cleandata/sample1_clean_R2.fq ... cleandata/samplen_clean_R1.fq cleandata/samplen_clean_R2.fq

#ワイルドカード指定
kraken2 --paired --threads 20 -db kraken_db --confidence 0.5 --output out_dir/ fastq_dir/sample*.fastq

本家の一部のオプションは対応していない。そしてエラーメッセージは常に”ペアのリードを指定してください”と出るが、対応していないオプションを指定した時もこのメッセージが出るので注意。別のエラーである可能性がある。

 

目的の分類群が存在するかしないか調査のため、ゼロカウントも全て出力。

mkdir kraken2_taxonomic_profiling
kraken2 --use-names --report-zero-counts --paired --threads 20 -db <path>/<to>/kraken2-full-database/ --output kraken2_taxonomic_profiling/ --report report.txt fastq_dir/SRR*_{1,2}.fastq

 

kraken2_taxonomic_profiling/

カレントパスにはreport.txtが出力される。

 

コメント

サンプル名がSRAなどの番号だと、コマンドへのファイル名提供方法によっては奇数サンプルのファイルしか出力されなくなる。ペアエンドはあらかじめマージし、--pairedなしで指定したほうがよいかもしれない。

 

  • folkはv2.0.9が最新。本家kraken2のより新しいv2.1.3で作成したDBを読み込むと、ロード途中でフリーズした。(間違い)
  • より新しいversionのkrakenで作ったデータベースでも、このkraken2 folkインストール時に./install_kraken2.shコマンドで指定して作っていれば認識する。複数ファイルを認識しないなら、./install_kraken2.shで指定したデータベース内にできるkraken2が使われていないことをまず疑う。
  • 本folkより新しいversionのkraken2で作ったtranslateデータベースも認識する。
  • v2.0.9でデータベースをダウンロードしてビルドするとエラーが起きるようになった。原因はダウンロードURLが変更されているため。rsync_from_ncbi.plがダウンロードスクリプトなので、findでrsync_from_ncbi.plのパスを確認し(find ~/mambaforge/ -name rsync_from_ncbi.pl)し、v2.1.3のrsync_from_ncbi.plと置き換えてからビルドすれば回避できる。

引用

GitHub - daydream-boost/kraken2: The second version of the Kraken taxonomic sequence classification system

 

The second version of the Kraken taxonomic sequence classification system

https://github.com/DerrickWood/kraken2

 

参考

https://github.com/DerrickWood/kraken2/issues/87

 

*1

最新のEPYC計算機構成でRAM diskを試す事も思いつくが、いずれにせよRAM diskは共用の計算機環境では使えない