メタゲノムコンティグのクラスタリング（binning）を行う MetaDecoder

　メタゲノムコンティグを潜在的なゲノムにクラスタリングすることは、微生物集団の機能的役割を調査する上で重要なステップである。既存のアルゴリズムは、シミュレーションまたは実際のシーケンスデータセットでかなりの成功を収めている。しかし、複雑なメタゲノムからのコンティグを正確に分類することは依然として課題である。本著者らは、k-merの頻度とカバレッジに基づいてメタゲノムコンティグを分類できる新規クラスタリングアルゴリズム「MetaDecoder」を導入した。MetaDecoderは2層モデルとして構築され、第1層はGPUベースの修正ディリクレ過程ガウス混合モデル（DPGMM）である。これは各DPGMMクラスターの重みを制御し、小さなクラスター内のコンティグを動的に解消して残存クラスターに再割り当てすることで、過度な分割を回避する。第2層は、半教師ありk-mer頻度確率モデルと、単一コピーマーカー遺伝子に基づくカバレッジをモデル化する修正ガウス混合モデルで構成される。シミュレーションデータと実データを用いたベンチマークにより、MetaDecoderがメタゲノムコンティグの効果的クラスタリング手法として有望であることが実証された。

メタゲノムコンティグの効果的なクラスタリングと微生物データからの微生物群集再構築を目的としたGPUベースのMetaDecoderを開発した。シミュレーションデータと実データの両方へのMetaDecoder適用により、より完全なクラスターを生成し、汚染を低減できることが実証された。MetaDecoderを用いて、新規の高品質ゲノムを同定し、既存の細菌ゲノムカタログを拡充した。

インストール

condaで環境を作って導入した（ubuntu22）。GPUにも対応している(レポジトリ参照)｡

Github

mamba create -n MetaDecoder python=3 -y
conda activate MetaDecoder
#依存
pip3 install numpy scipy scikit-learn threadpoolctl

#本体
pip3 install -U https://github.com/liu-congcong/MetaDecoder/releases/download/v1.2.1/metadecoder-1.2.1-py3-none-any.whl

> metadecoder -h

usage: metadecoder [-h] [-v] {coverage,seed,cluster} ...

The MetaDecoder algorithm.

options:

-h, --help show this help message and exit

-v, --version show program's version number and exit

program:

{coverage,seed,cluster}

Run metadecoder coverage|seed|cluster --help for help.

coverage Calculate sequence coverages from BAM files.

Usage: metadecoder coverage -s sample1.bam ... sampleN.bam -o output.coverage.

seed Map the marker genes (HMMs) to sequences.

Usage: metadecoder seed -f input.fasta -o output.seed.

cluster Running the MetaDecoder to cluster sequences.

Usage: metadecoder cluster -f input.fasta -c input.coverage -s input.seed -o output.metadecoder

Liucongcong congcong_liu@icloud.com.

> metadecoder coverage -h

usage: metadecoder coverage [-h] -b str [str ...] -o str [--bin_size int] [--mapq int] [--aligned float] [--threads int]

options:

-h, --help show this help message and exit

-b, --bam str [str ...]

The bam formatted files with headers "@SQ".

-o, --output str The output coverage file.

--bin_size int The size (bp) of each bin.

The value should be a positive integer.

Default: 2000.

--mapq int The MAPQ threshold.

The value should be a non-negative integer.

Default: 0.

--aligned float The minimum aligned ratio calculted from CIGAR.

The value should be from 0 to 1.

Default: 0.95.

--threads int The number of threads.

The value should be a positive integer.

Default: 10.

> metadecoder seed -h

usage: metadecoder seed [-h] -f str -o str [--threads int] [--coverage float] [--accuracy float]

options:

-h, --help show this help message and exit

-f, --fasta str The fasta formatted assembly file.

-o, --output str The output seed file.

--threads int The number of threads.

The value should be a positive integer.

Default: 20.

--coverage float The min coverage of domain.

The value should be from 0 to 1.

Default: 0.5.

--accuracy float The min accuracy.

The value should be from 0 to 1.

Default: 0.6.

> metadecoder cluster -h

usage: metadecoder cluster [-h] -f str -c str [str ...] -s str -o str [--output_unclustered_sequences] [--clustering_probability float] [--weight float] [--kmer int] [--min_cluster_size int] [--min_dpgmm_size int] [--max_dpgmm_distance float] [--min_sequence_length int] [--outlier float]

[--sampling_length1 int int] [--sampling_number1 int] [--sampling_length2 int int] [--sampling_number2 int] [--random_number int] [--disable_gpu] [--no_clusters]

options:

-h, --help show this help message and exit

-f, --fasta str The fasta formatted assembly file.

-c, --coverage str [str ...]

The coverage files generated by "metadecoder coverage".

-s, --seed str The seed file generated by "metadecoder seed".

-o, --output str The prefix of output files.

--output_unclustered_sequences

Generate "*.unclustered.fasta" file for all unclustered sequences.

Default: False.

--clustering_probability float

The minimum clustering probability.

The value should be from 0.0 to 1.0

Default: 0.0.

--weight float The scale factor of the weight in regularization term.

The coverage probabilistic model uses "--weight" / "number of sequencing samples" as the weight of priors.

The value should be a non-negative float.

Default: 1.

--kmer int The length of kmer.

The value should be a positive integer.

Default: 4.

--min_cluster_size int

The minimum size of a cluster to output.

The value should be a positive integer.

Default: 200000.

--min_dpgmm_size int The minimum size of a DPGMM cluster.

DPGMM stops shrinking the number of clusters, until the size of the smallest cluster is greater than or equal to this value.

The value should be a positive integer.

Default: 500000.

--max_dpgmm_distance float

Sequences in a DPGMM cluster will be removed, if the average Euclidean distance of all paired kmer frequencies is greater than this value.

We recommend that you do not change this value.

The value should be a positive float.

Default: 0.04.

--min_sequence_length int

Sequences with the length being greater than or equal to this value can be involved in the MetaDecoder algorithm.

The value should be a integer greater than 2000.

Default: 2500.

--outlier float The proportion of outliers that will be removed by the isolation forest.

The value should be from 0 to 0.5.

Default: 0.1.

--sampling_length1 int int

The range of the length of sub-seeds sampled from each seed used to train the classifier.

The value should be two positive integers from low to high.

Default: 1500 2500.

--sampling_number1 int

The number of sub-seeds sampled from each seed used to train the classifier.

The value should be a positive integer.

Default: 20.

--sampling_length2 int int

The range of the length of sub-seeds sampled from each extended seed used to train the classifier.

The value should be two positive integers from low to high.

Default: 3000 4000.

--sampling_number2 int

The number of sub-seeds sampled from each extended seed used to train the classifier.

The value should be a positive integer.

Default: 300.

--random_number int The random number generator seeded by the given integer

Default: 0.

--disable_gpu Disable the GPU if it does not have enough memory to accommodate the data.

Default: False.

--no_clusters Do not output fasta format cluster files.

Default: False.

実行方法

１、contig のカバレッジの計算｡MAGにマッピングして得たbamファイルを指定する(bam作成例 *1)｡

#シングルbam指定
metadecoder coverage -b SAMPLE1.bam --threads 10 -o METADECODER.COVERAGE

#複数bam指定
metadecoder coverage -b SAMPLE1.bam SAMPLE2.bam SAMPLE3.bam --threads 10 -o METADECODER.COVERAGE

-b The bam formatted files with headers "@SQ".
-o The output coverage file.
--bin_size The size (bp) of each bin. The value should be a positive integer. Default: 2000.
--threads The number of threads. The value should be a positive integer. Default: 10.

出力例

METADECODER.COVERAGE

2、シードとなるシングルコピーマーカー遺伝子をコンティグにマッピングする｡MAG配列を指定する(single-copy marker gene はクラスタリングの初期ガイドとして使われる)｡megahitからのMAGを使う場合､名前がユニークでないというエラーが出ることがあるのでリネームして使うのが無難*2｡

metadecoder seed --threads 20 -f ASSEMBLY.fasta -o METADECODER.SEED

-f The fasta formatted assembly file.
-o The output seed file.
--threads The number of threads. The value should be a positive integer. Default: 20.

出力例

3、クラスタリング｡1と2の出力を指定する｡

metadecoder cluster -f ASSEMBLY.fasta -c METADECODER.COVERAGE -s METADECODER.SEED -o METADECODER

#複数サンプルのカバレッジファイルを指定
metadecoder cluster -f ASSEMBLY.fasta -c *.METADECODER.COVERAGE -s METADECODER.SEED -o METADECODER