🧬
GenomeFace: 43,000個の微生物ゲノムでトレーニングされたディープラーニングベースのメタゲノムビナー
GenomeFace: a deep learning-based metagenome binner trained on 43,000 microbial genomes
GenomeFace: 43,000個の微生物ゲノムでトレーニングされたディープラーニングベースのメタゲノムビナー
はじめに
ビニングはメタゲノムアセンブルゲノム(MAG)を作成する際に、同じゲノム由来と予測される配列をグルーピングする操作です。
http://onlinelibrary.wiley.com/doi/full/10.1002/mbo3.1298
微生物群集の機能、相互作用、進化のダイナミクスを理解するためには、ゲノムの高精度な再構築が欠かせません。
GenomeFaceは、DNAの構成と環境内の生物種の多様性という2つの要素に基づくニューラルネットワークを活用し、ディープラーニングによるビニングアプローチを提供します。
この手法では、ニューラルネットワークによって生成された埋め込みをクラスタリングし、マーカー遺伝子（marker genes）を活用して最適なシーケンスクラスターを選定します。本アプローチは、4つのシミュレートデータセットで評価され、従来の計算負荷が高いツールと比較して、より多くのほぼ完全なゲノムを高速に再構築できることが示されました。
さらに、3つのケースでテストを行い、スケーラビリティと精度を実証しました。他のビニングツールと比較すると、ほぼ完全なゲノムを47〜183%多く生成し、3,000以上の新たな候補種のゲノムを特定しました。これは、既知の細菌の系統樹が約4%拡大する規模に相当します。
また、2025年2月にBiorxivで公開されたディープラーニングベースのビニングツールのベンチマーク論文では、COMEBinとGenomeFaceが精度面で一貫して高いパフォーマンスを示しました。加えて、MetaBAT2と比較して、GenomeFaceは処理速度の面でも優れた成績を記録しました。こうした背景から、現在開発中のMetaBAT3では、GenomeFaceの技術を統合する計画も進められているようです。
以下に、GenomeFaceの実行に必要な準備と実行例を示します。
GenomeFaceのインストール
GenomeFace はパッケージ キャッシュに約 4GB、環境用にさらに 10GB程度必要とすることに注意が必要です。また、GenomeFaceはCUDAを使用するためNVIDIAのGPUが必要になることにも注意してください。
mambaで仮想環境を作成してインストールします。
mamba create -p $PWD/genomeface \
genomeface cuda-version=11.7 \
-c nvidia -c rapidsai -c conda-forge -c bioconda -c https://portal.nersc.gov/project/m1982/genomeface_repo
下準備
下記のステップで作業を進めます。
- 環境構築
- Fastqの取得
- KneaddataによるデータQC
- megahitによるショートリードアセンブリ
- jgi_summarize_bam_contig_depthsによるcoverageの計算
- pyrodigalとpyhmmrを用いたSingle copy merker遺伝子の予測
上記5ステップが終了したらGenomeFaceでビニングを実行します。
1. 環境構築
mambaで必要なツールをインストールした仮想環境を作成します。
mamba create -p $PWD/prep_genomeface -c bioconda -c conda-forge \
pyhmmer pyrodigal kneaddata megahit metabat2 strobealign samtools numpy biopython pandas
kneaddata用にReleasesページからTrimmomaticの実行ファイルをダウンロード&解凍しておきます。
wget https://github.com/usadellab/Trimmomatic/files/5854859/Trimmomatic-0.39.zip
unzip Trimmomatic-0.39.zip
rm Trimmomatic-0.39.zip
# 確認
java -jar Trimmomatic-0.39/trimmomatic-0.39.jar -version
2. Fastqの取得
SRA Explorerでアクセッション番号を検索し、チェックボックスにチェック ⇢ Add X to collectionをクリックします。
右上のX saved datasetsよりBash script for downloading FastQ filesに記載されているcurlのコマンドフレーズを実行してfastqファイルを取得します。
以下はその例です。
mkdir 00RawFastq
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/002/SRR18585102/SRR18585102_1.fastq.gz -o 00RawFastq/SRR18585102_R1.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/002/SRR18585102/SRR18585102_2.fastq.gz -o 00RawFastq/SRR18585102_R2.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/013/SRR18585213/SRR18585213_1.fastq.gz -o 00RawFastq/SRR18585213_R1.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/013/SRR18585213/SRR18585213_2.fastq.gz -o 00RawFastq/SRR18585213_R2.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/020/SRR18585220/SRR18585220_1.fastq.gz -o 00RawFastq/SRR18585220_R1.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/020/SRR18585220/SRR18585220_2.fastq.gz -o 00RawFastq/SRR18585220_R2.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/091/SRR18585191/SRR18585191_1.fastq.gz -o 00RawFastq/SRR18585191_R1.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/091/SRR18585191/SRR18585191_2.fastq.gz -o 00RawFastq/SRR18585191_R2.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/019/SRR18585219/SRR18585219_1.fastq.gz -o 00RawFastq/SRR18585219_R1.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/019/SRR18585219/SRR18585219_2.fastq.gz -o 00RawFastq/SRR18585219_R2.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/015/SRR18585215/SRR18585215_1.fastq.gz -o 00RawFastq/SRR18585215_R1.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/015/SRR18585215/SRR18585215_2.fastq.gz -o 00RawFastq/SRR18585215_R2.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/018/SRR18585218/SRR18585218_1.fastq.gz -o 00RawFastq/SRR18585218_R1.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/018/SRR18585218/SRR18585218_2.fastq.gz -o 00RawFastq/SRR18585218_R2.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/017/SRR18585217/SRR18585217_1.fastq.gz -o 00RawFastq/SRR18585217_R1.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/017/SRR18585217/SRR18585217_2.fastq.gz -o 00RawFastq/SRR18585217_R2.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/016/SRR18585216/SRR18585216_1.fastq.gz -o 00RawFastq/SRR18585216_R1.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/016/SRR18585216/SRR18585216_2.fastq.gz -o 00RawFastq/SRR18585216_R2.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/014/SRR18585214/SRR18585214_1.fastq.gz -o 00RawFastq/SRR18585214_R1.fastq.gz
curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR185/014/SRR18585214/SRR18585214_2.fastq.gz -o 00RawFastq/SRR18585214_R2.fastq.gz
3. KneaddataによるデータQC
KneaddataのリファレンスDBがT2Tに対応していたのでそちらを使用します。以下のコマンドでリファレンスDBをダウンロードします。
mkdir 01Ref
wget https://huttenhower.sph.harvard.edu/kneadData_databases/Homo_sapiens_hg39_T2T_Bowtie2_v0.1.tar.gz \
-o 01Ref/Homo_sapiens_hg39_T2T_Bowtie2_v0.1.tar.gz
# 解凍
tar xvf 01Ref/Homo_sapiens_hg39_T2T_Bowtie2_v0.1.tar.gz -C 01Ref/
rm 01Ref/Homo_sapiens_hg39_T2T_Bowtie2_v0.1.tar.gz
以下のコマンドでkneaddataを実行して、ホストゲノムを除去します。
mkdir 02kneaddata
array=(SRR18585102 SRR18585213 SRR18585220 SRR18585191 SRR18585219 SRR18585215 SRR18585218 SRR18585217 SRR18585216 SRR18585214)
for acc in ${array[@]};do
kneaddata \
--trimmomatic Trimmomatic-0.39/ \
-t 32 \
-i1 "00RawFastq/${acc}_1.fastq.gz" \
-i2 "00RawFastq/${acc}_2.fastq.gz" \
--reference-db "01Ref/bowtie2-index" \
--output-prefix "${acc}" \
--log "${acc}"_kneaddata.log \
--trimmomatic-options "ILLUMINACLIP:$PWD/Trimmomatic-0.39/adapters/NexteraPE-PE.fa:2:30:10 SLIDINGWINDOW:4:20 MINLEN:50" \
-o "02kneaddata/${acc}"
pigz -p 32 "02kneaddata/${acc}/*.fastq"
done
4. megahitによるショートリードアセンブリ
ビニング対象のconfig.fastaをmegahitで作成します。KneadDataでQCとホスト除去済のデータがインプットファイルになります。
megahitによるショートリードアセンブリは以下のコマンドを実行します。
mkdir 03Assemble
array=(SRR18585102 SRR18585213 SRR18585220 SRR18585191 SRR18585219 SRR18585215 SRR18585218 SRR18585217 SRR18585216 SRR18585214 )
for acc in ${array[@]}; do
mamba run -p $PWD/prep_genomeface megahit \
-1 02kneaddata/${acc}/${acc}_paired_1.fastq.gz \
-2 02kneaddata/${acc}/${acc}_paired_2.fastq.gz \
-t 32 \
-o 03Assemble/${acc}
done
5. jgi_summarize_bam_contig_depthsによるcoverageの計算
strobealignでcontigにmappingしてsamtoolsでcoordinate sort後、metabat2に内包されているjgi_summarize_bam_contig_depthsでcoverageを計算します。
mkdir 03Assemble/strobealign_results
array=(SRR18585102 SRR18585213 SRR18585220 SRR18585191 SRR18585219 SRR18585215 SRR18585218 SRR18585217 SRR18585216 SRR18585214 )
for acc in ${array[@]}; do
mamba run -p $PWD/prep_genomeface strobealign \
-t 32 \
--eqx \
03Assemble/${acc}/final.contigs.fa \
02kneaddata/${acc}/${acc}_paired_1.fastq.gz \
02kneaddata/${acc}/${acc}_paired_2.fastq.gz | \
$PWD/prep_genomeface/bin/samtools view -bS -@ 8 | \
$PWD/prep_genomeface/bin/samtools sort -O BAM -@ 8 -o 03Assemble/strobealign_results/${acc}.bam
samtools index 03Assemble/strobealign_results/${acc}.bam
done
jgi_summarize_bam_contig_depthsによるcoverageの計算は以下のコマンドを実行します。
mkdir 03Assemble/jgi_summarize_bam_contig_depths_results
array=(SRR18585102 SRR18585213 SRR18585220 SRR18585191 SRR18585219 SRR18585215 SRR18585218 SRR18585217 SRR18585216 SRR18585214 )
for acc in ${array[@]}; do
mamba run -p $PWD/prep_genomeface jgi_summarize_bam_contig_depths \
--outputDepth 03Assemble/jgi_summarize_bam_contig_depths_results/${acc}_assembly_depth.tsv \
03Assemble/strobealign_results/${acc}.bam
done
6. pyrodigalとpyhmmrを用いたSingle copy marker遺伝子の予測
ベンチマーク論文で使用されているスクリプトを利用します。
git cloen https://github.com/soedinglab/ビニング_benchmarking.git
エラーが発生したので以下の部分を修正しました。
def predict_genes(entries):
""" predict genes """
annotations = []
proteins = []
for entry in entries:
genes = gene_finder.find_genes(str(entry.seq).encode())
for idx, gene in enumerate(genes):
protein = pyhmmer.easel.TextSequence(name=(str(idx) + "_" + str(entry.id)).encode(),
sequence=gene.translate()).digitize(alphabet)
proteins.append(protein)
hits = pyhmmer.hmmsearch(hmms, proteins, cpus=1)
for hit in hits:
for domain_hit in hit:
if domain_hit.included and not domain_hit.duplicate:
p_nm = domain_hit.name.decode()
record_id = p_nm[p_nm.index('_')+1:]
- if domain_hit.score >= cutoff[hit.query_name.decode()]:
+ if domain_hit.score >= cutoff[hit.query.name.decode()]:
annotations.append(\
- Result(p_nm, record_id, hit.query_name.decode(), domain_hit.score))
+ Result(p_nm, record_id, hit.query.name.decode(), domain_hit.score))
以下のコマンドでSingle copy marker遺伝子の予測を実行します。
mkdir 03Assemble/pyrodigal_results
array=(SRR18585102 SRR18585213 SRR18585220 SRR18585191 SRR18585219 SRR18585215 SRR18585218 SRR18585217 SRR18585216 SRR18585214 )
for acc in ${array[@]}; do
mkdir 03Assemble/pyrodigal_results/${acc}
mamba run -p $PWD/prep_genomeface python ビニング_benchmarking/util/pyrodigal_prediction.py \
--seq 03Assemble/${acc}/final.contigs.fa \
--outdir 03Assemble/pyrodigal_results/${acc}
done
GenomeFaceの実行
USAGE GenomeFace
usage: genomeface [-h] -i INPUT_FILE -a ABUNDANCE_FILE -g MARKERS_FILE -o
OUTPUT_FOLDER [-m MIN_CONTIG_LENGTH] [-p MIN_PURITY] [-s]
[-b]
GenomeFace Prerelease
- A next-generation tool for metagenomic ビニング, using deep learning and
multi-GPU accelerated clustering. Ideal for large-scale, real-world data.
The Exabiome Project (Lawrence Berkeley National Laboratory)
- Contact rlettich@lbl.gov for issues or unexpected poor performance.
optional arguments:
-h, --help show this help message and exit
-i INPUT_FILE Input FASTA file containing metagenome assembly
(optionally gzipped).
-a ABUNDANCE_FILE MetaBAT 2 style TSV file containing abundance data.
Typically produced by the `jgi_summarize_depths` or
`coverm` programs.
-g MARKERS_FILE Input TSV file describing marker genes found on each
contig. Can be produced by the included `markersgf`
program.
-o OUTPUT_FOLDER Output folder for writing bin FASTA files
-m MIN_CONTIG_LENGTH Minimum contig length to be considered for ビニング
(default: 1500).
-p MIN_PURITY Minimum marker gene estimated % purity for selecting
clusters for output. Balances Precsion / Recall
(default: 85)
-s Specifies that the input FASTA is multiple single
sample assemblies, concatenated. Contig names should
be Prefixed by per assembly name With ending in 'C'.
e.g. asm_oneCsequence5
-b Use Bayesian Optimizaiton to Optimize Balance Between
Compositional and Abundance Distances (Default:
False).
CUDA Devices (Clustering Acceleration):
- Device 0: NVIDIA GeForce RTX 3070 Ti
Tensorflow Enabled Devices (Neural Network Acceleration):
- /physical_device:GPU:0
Notes:
- By default, coassembly is assumed. For concatenated single-sample assemblies, use the [-s] flag.
Examples:
genomeface -i coassembly.fa.gz -a abundance.tsv -g markers.tsv -o ./output
genomeface -i concatenated_single_sample_assemblies.fa.gz -a abundance.tsv -g markers.tsv -o ./output -s -m 1000
GenomeFaceによるビニングは以下のコマンドで実行します。
mkdir 04Genomeface
array=(SRR18585102 SRR18585213 SRR18585220 SRR18585191 SRR18585219 SRR18585215 SRR18585218 SRR18585217 SRR18585216 SRR18585214)
for acc in ${array[@]}; do
mamba run -p $PWD/genomeface genomeface \
-i 03Assemble/${acc}/final.contigs.fa \
-a 03Assemble/jgi_summarize_bam_contig_depths_results/${acc}_assembly_depth.tsv \
-g 03Assemble/pyrodigal_results/${acc}/marker_hits \
-o 04Genomeface/${acc} \
-m 1000
done
04Genomefaceフォルダ下に作成されたアクセッションナンバーがGenomeFaceの出力フォルダです。
数値のみのファイル名のファイルはコンティグ配列がまとめられたFasta形式のファイルです。
>k141_65097
TCTATAAGACCAGCATGTGTCCGCTTGCGAAGGTCATGCGCCGCGAGTTAAAGAAGCGCGGAGTGAAGAAATTAAAGGTGGTTTATTCCACAG...
>k141_14206
TTACTTGCAAAATATATATTTTTAGTTGAACTTAAAGATTGCAGTCCCGTACGCCGGCAGCGGATACGCAAAGGAGTATGGCTGGCCATCGCACTCGC...
bins.tsvはBin Fastaファイルとcontig名の対応表です。最小コンティグ配列長や Precsion / Recallの値で足切りされた結果のBin Fastaが出力されているので、対応表にあるすべてのコンティグがBin Fastaとして出力されるわけではありません。
0 k141_14150
1 k141_36791
2 k141_87729
3 k141_33961
4 k141_5662
5 k141_14151
6 k141_84900
7 k141_50942
8 k141_42450
9 k141_2833
10 k141_65094
11 k141_53773
12 k141_67921
5402 k141_39624
12406 k141_59437
15 k141_31133
16 k141_73581
17 k141_84901
17035 k141_11323
19 k141_2835
20 k141_79243
21 k141_36792
22 k141_31134
23 k141_25479
...
以上がGenomeFace実行と結果取得までの流れです。
Discussion