タイトルの通りの内容。kraken2の公式サイトにはGTDBのR220バージョンのデータベースが用意されていて､kraken2の菌叢解析にGTDB DBのデータベースを利用できるようになった｡しかし､20205年6月時点で最新のR226はまだ用意されていない =>
追記: 6月9日に公開されました。
(AWS Public Dataset Programの支援によって公開されている､Centrifuge , bowtie, hisatなどのpre-build indexも利用できる)
最新DBはブログを書いている間に公開されたが､kraken2のR226 DBを作成する手順を確認しておく。
１，R226のためのtaxdumpファイル群を作る必要がある。ここではgtdb_to_taxdumpを使わせてもらう。
mamba create -n gtdb-env python=3.10 -y
conda activate gtdb-env
pip install gtdb_to_taxdump
gtdb_to_taxdump.pyの実行。
https://data.gtdb.ecogenomic.org/releases/release226/226.0/
gtdb_to_taxdump.py \
https://data.gtdb.ecogenomic.org/releases/release226/226.0/ar53_taxonomy_r226.tsv.gz \
https://data.gtdb.ecogenomic.org/releases/release226/226.0/bac120_taxonomy_r226.tsv.gz \
> taxID_info.tsv
dumpファイル群ができる（写真はR226より昔のもの）。
２、dumpファイル群はfastaを配置予定のディレクトリのtaxonom/に配置する。
mkdir -p kraken2_R226DB/taxonomy
mv *dmp kraken2_R226DB/taxonomy/
３、GTDB fastaの準備。R226のgtdbtk_packageをダウンロードする。あるいはfastaファイルだけをダウンロードする。ここではgtdbtk_package全体をダウンロードして中のfastaを移動させる。それからfasta名をIDだけにリネームする。
https://ecogenomics.github.io/GTDBTk/installing/index.html
wget https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/auxillary_files/gtdbtk_package/full_package/gtdbtk_data.tar.gz
tar -xzvf gtdbtk_data.tar.gz
mkdir original_fasta
find gtdbtk_data/skani/database/GCA/ -type f -name "*.fna.gz" -exec mv {} original_fasta/
find gtdbtk_data/skani/database/GCF/ -type f -name "*.fna.gz" -exec mv {} original_fasta/
#rename fasta
cd original_fasta/
find . -name "*.gz" > list
cat list | while read line; do
newname=$(basename "$line" .genomic.fna.gz).fna.gz
mv "$line" "$(dirname "$line")/$newname"
done
4、kraken2でビルドするゲノムの各ヘッダーには、kraken2 TaxIDが付与されている必要がある。
例
original fasta
> head GCA_000008085.1.fna
header renamed fasta
> head GCA_000008085.1.fna
bashスクリプトを使ってfastaにIDを付与する (*1)｡カレントにstripped_names.txtを作成､
カレントに入力fastaディレクトリfasta_dirを配置して､スクリプトを実行する｡
cat kraken2_R226DB/taxonomy/names.dmp |sed 's/RS_//' |sed 's/GB_//' > stripped_names.txt
#実行権限をつけてラン
chmod +x add_taxid.sh
./add_taxid.sh
作成された"fasta_with_taxid"ディレクトリにtax ID付きのfastaファイルが書き出される｡
5､names.dmpから参照してtax IDを付与する。さきほど変換してtaxonomyに配置したnames.dmpからルックアップする(修正 proteinタグを除去)。
cat kraken2_R226DB/taxonomy/names.dmp |sed 's/RS_//' |sed 's/GB_//' > stripped_names.txt
mkdir kraken2_protein_db
＃1個ずつデータベースに追加していく､4で得たfasta_with_taxidディレクトリを指定
for f in fasta_with_taxid/*.fna; do
kraken2-build --add-to-library "$f" --db kraken2_protein_db --threads 20 --no-masking
done
#最後にビルドする(数時間はかかる)
kraken2-build --build --db kraken2_R226DB/ --threads 20 --no-masking
ビルドできたら下の.faaファイル(ヘッダーをリネームしたもの)などはすべて不要になる｡
6､ビルドできたら作ったDBを指定してkraken2を実行する｡kraken2マルチサンプル対応folkを使うなら､gzip解凍機能はないため､fastq.gzは事前に解凍してから使用する(高速にたくさんのfastqをプロファイリングするため生のfastqを要求する)
#1､ fq.gzの並列解凍(大容量SSD推奨)
parallel -j 10 gzip -d ::: *.fastq.gz &
#注 ペアの片方だけ存在するとkraken2 folkのpairedモードはエラーを出すので注意する
#2 kraken2 folkの実行(kraken2 folkはkraken2 DBを指定してビルドしたはず,
folkのバイナリはそのディレクトに配置されるのでそこにパスを通す)
conda activate kraken2folk
export PATH=/path/to/kraken2_R226DB/:$PATH
# run paired-endの場合
mkdir kraken2_result
kraken2 --paired --threads 20 -db kraken2_R226DB/ --confidence 0.2 --output kraken2_result/ *.fastq
注: kraken2マルチサンプル対応folkは､一部のfastqのペアの同期が壊れていると正しく動作しない｡大量のfastqを指定するときは､kraken2実行前にペアの同期が正常なのチェックするスクリプトを走らせておいたほうがいい｡
注意点
- kraken2のメモリ使用量はindexサイズとほぼ同じくらいとなる｡R226はindex サイズ644.0GBなので､メモリは650GB程度は要求する(実際は1TBメモリのマシンでピークメモリが670GB近く使った)｡メモリが足りない場合､メモリ使用量を低減するオプションを使わないとランできない｡
- GTDB R226データベースには､種を代表するゲノムが143,613個(要確認)含まれている｡これはGTDBに含まれる数十万以上のゲノムから種を代表するゲノムとして､1つずつ選ばれたものとなる｡そのうち､単離株はNCBIのrepresentative genome(約21000個)と同じ株(品質が低いもの以外はタイプ株)のゲノムが使われていると思われる｡残りのsp.株についてはUnculturedの仮想種のゲノムが使われている｡すべてのペア間はANI 95% species boundaryを満たしている｡
- GTDBの種代表ゲノムはDB更新に従って加速的に増している｡おそらく今後のGTDB更新でもまだまだ原核生物の種数は増えそうに思える(論文によって京を超える数から数十万と控えめまで幅広い予想がある)｡この14万3千個のゲノムのindexでもメモリを670GBほど使うことを考えると､GTDBの今後のバージョンアップでさらに種代表ゲノムが増えたとき､kraken2をGTDBと組み合わせて使うのは厳しくなるかもしれない(xeon と組わせる前提で､optaneサーバーメモリが2TB程度は比較的安く購入できるが､そもそもメモリとしては遅く､ハードウエアの誓約が多い)｡
- proteinモードのDBを使う場合､--confidenceのデフォルトに注意する｡Kraken 2 は protein DB では 6 フレーム翻訳を行い、15 aaなどの minimizer でリファレンスprotein DB と照合する。そのため､confidenceをデフォルトの0から0.2以上にすると､ 6 フレーム翻訳によってマッチ k-merの比率が1から大きく減ってしまう｡confidence=0.2以上で試してみると､既知ゲノムでも分配されないfastqが99％以上になった｡
- DBが巨大なため､高速なPCIe-SSDにDBを配置しても､DBの事前読み込みだけで10分近い時間がかかる(ハードウェアによって異なる)
関連
*1
add_taxid.shスクリプト
==================================================
NAMES_FILE="stripped_names2.txt"
IN_DIR="fasta_dir"
OUT_DIR="fasta_with_taxid"
mkdir -p "$OUT_DIR"
# 各 fasta ファイルを処理
for fasta in "$IN_DIR"/*.fna; do
acc=$(basename "$fna" .fna)
# taxid を抽出
taxid=$(awk -F '\t' -v acc="$acc" '$2 == acc { gsub(/[^0-9]/, "", $1); print $1 }' "$NAMES_FILE")
if [ -z "$taxid" ]; then
echo "taxid not found in $acc"
continue
fi
echo "taxid: $taxid assigned in $acc"
#ヘッダーに taxid を付加
sed "s/^>/>kraken:taxid|${taxid} /" "$faa" > "$OUT_DIR/${acc}.faa"
done
