ロングリードシーケンスは、複雑な微生物群集の特徴を明らかにする上で大きな可能性を秘めているが、ロングリード専用に設計された分類学的プロファイリングツールはまだ不足している。ここでは、ロングリードのユニークな特性を生かしたマーカーベースの分類学的プロファイラであるMelonを紹介する。Melonは、計算時間を短縮するために2段階の分類スキームを採用し、あいまいなリードを処理するために期待値最大化ベースのポスト補正モジュールを備えている。Melonは、模擬サンプルとシミュレーションサンプルの両方において、他の既存ツールと比較して優れた性能を達成した。廃水メタゲノミックサンプルを用いて、Melonが全ゲノムコピーの信頼性の高い推定値、および種レベルの分類学的プロファイルを提供できることを示し、Melonの応用可能性を実証した。
インストール
ubuntu20.04で環境を作ってテストした。
https://github.com/xinehc/melon
#conda
mamba create -n melon -c conda-forge -c bioconda melon
conda activate melon
> melon -h
usage: melon -d DIR -o DIR [-t INT] [-k DIR] [--skip-profile]
[--skip-clean] [-m INT] [-e FLOAT] [-i FLOAT] [-s FLOAT] [-n INT] [-p
FLOAT] [-a INT] [-c FLOAT] FILE [FILE ...]
Melon: metagenomic long-read-based taxonomic identification and quantification
positional arguments:
FILE Input fasta <*.fa|*.fasta> or fastq
<*.fq|*.fastq> file, gzip optional <*.gz>.
required arguments:
-d DIR, --db DIR Unzipped database folder, should contains
<prot.fa>, <nucl.*.fa> and <metadata.tsv>.
-o DIR, --output DIR Output folder.
optional arguments:
-t INT, --threads INT
Number of threads. [128]
-k DIR, --db-kraken DIR
Unzipped kraken2 database for pre-filtering of
non-prokaryotic reads. Skip if not given.
--skip-profile Skip profiling, output only estimated total
genome copies.
--skip-clean Skip cleaning, keep all temporary <*.tmp> files.
additional arguments:
-m INT Max. number of target sequences to report
(--max-target-seqs/-k in diamond). [25]
-e FLOAT Max. expected value to report alignments
(--evalue/-e in diamond). [1e-15]
-i FLOAT Min. identity in percentage to report
alignments (--id in diamond). [0]
-s FLOAT Min. subject cover to report alignments
(--subject-cover in diamond). [75]
-n INT Max. number of secondary alignments to report
(-N in minimap2). [2147483647]
-p FLOAT Min. secondary-to-primary score ratio to
report secondary alignments (-p in minimap2). [0.9]
additional arguments for EM:
-a INT Terminal condition - max. iteration. [1000]
-c FLOAT Terminal condition - epsilon (precision). [1e-10]
データベース
wget -q --show-progress https://figshare.com/ndownloader/files/42694969/database.tar.gz
tar -zxvf database.tar.gz
cpu_count=$(python -c 'import os; print(os.cpu_count())')
diamond makedb --in database/prot.fa --db database/prot --quiet
#minimap2のindex作成
ls database/nucl.*.fa | sort | xargs -P $cpu_count -I {} bash -c '
filename=${1%.fa*}; \
filename=${filename##*/}; \
minimap2 -x map-ont -d database/$filename.mmi ${1} 2> /dev/null' - {}
＃不要なファイルを除去
rm -rf database/*.fa
テストラン
Loman LabのR10.3モックサンプルからランダムに選択され、品質管理された10,000の原核生物リード（（Porechopとnanoqで処理され、minimap2で真菌およびその他のリードが除去された）のラン。
wget -q --show-progress https://figshare.com/ndownloader/files/42847672/example.fa.gz
melon example.fa.gz -d database -o .
出力はJSONとTSV形式。
example.tsv
出力ファイル *.tsv には、個々の生物種の推定ゲノムコピー、対応する相対量、およびギャップ圧縮された ANI (マーカー遺伝子を含むリードと参照ゲノムクラスター間の平均塩基同一性) 値が含まれている。
実際に使う際は、品質管理され、汚染除去されたロングリードの使用が推奨されている。nanoq -q 10 -l 1000（最小品質スコア10、最小リード長1,000 bp）などで、低品質の生リードを除去する。サンプルにヒトDNAや既知の真核生物/ウイルスが多く含まれていることが判明している場合は、マッピングによる除去を考える。（コメント；汚染除去については、種の存在量を推定するため、最初から汚染と解釈されるリードは除いたほうが良いという事だと思います）
レポジトリより
- コンタミネーションの起源が不明な場合、または原核生物の平均ゲノムサイズを推定したい場合は、シンプルなプレフィルターモジュールを有効にする。
- プレフィルタリングモジュールを有効にするには、少なくともヒトと菌類を含むKrakenのデータベース（PlusPF、PlusPFP、またはそれらの上限付きバージョン）をダウンロードする（link）。
感想
知る限り、ショートリードだと塩基のkrakenやCentrifuge、KMCP、タンパク質レベルのkaijuなど豊富なプロファイリングツールが存在しますが、エラーが多く矛盾した系統シグナルを出しやすいロングリードの分類学的分類に関しては、実用的なツールがほとんど発表されていませんでした。このmelonはマーカーベースの分類プリファイラーで、最初にマーカーを含むリードを抽出し、その後アライメントベースのアプローチでコンパクト化したRefSeqゲノムデータベースにマッピングすることにより種レベルの分類学的存在量プロファイルを推定するようです。その意味ではmOTUs3やMetaPhlAn4に近いですが、ロングリードをネイティブに扱うための工夫があります。詳しくは論文を読んでください。使ってみて、とても使いやすく実用性も高いツールだなと感じました。
引用
Melon: metagenomic long-read-based taxonomic identification and quantification using marker genes
Chen Xi, Xiaole Yin, Xianghui Shi, Weifu Yan, Yu Yang, Lei Liu, Tong Zhang
bioRxiv, Posted December 18, 2023
関連
