Zenn
🦠

Autocycler:細菌ゲノムのロングリードコンセンサスアセンブリ

に公開
bioinfo
はじ環
ロングリード
バクテリアゲノムアセンブリ
tech

Autocycler: long-read consensus assembly for bacterial genomes

Autocycler:細菌ゲノムのロングリードコンセンサスアセンブリ

https://www.biorxiv.org/content/10.1101/2025.05.12.653612v1

はじめに

Autocyclerは、複数のロングリードアセンブリを自動的に統合し、高精度な細菌ゲノムのコンセンサスアセンブリを生成するコマンドラインツールです。

Trycyclerを基盤とし、ユーザーの手動操作を必要とせずに、アセンブリからコンパクトなDe Bruijnグラフを構築し、コンティグのクラスタリング、フィルタリング、オーバーラップのトリミングを行います。さらに、各遺伝子座で最も一般的なバリアントを選択してコンセンサス配列を決定します。

5つの細菌分離株のOxford Nanoporeリードを用いた評価では、Autocyclerは既存のアセンブラーや他の自動パイプライン、コンセンサスツールよりも優れた結果を示し、エラー率の低さと構造的な正確性の高さを実現しています。

単体アセンブラーの中では、CanuとFlyeが最も低い置換およびIndelエラー率を示しました。その他のツールでは、構造の不一致、プラスミドの欠落、環状コンティグの末端の重複などのアセンブリエラーがしばしば確認されました。

パイプラインでは、HybracterがPlassemblerを統合しているためか、プラスミドの欠落や重複といった構造エラーを回避していました。一方、プラスミドが大きい種（Enterobacter、Klebsiella）では、Unicyclerの影響か、プラスミドアセンブリ中に導入されたと思われるエラーが確認されました。

DragonflyeはFlye単体よりもエラー率が高かったものの、これはRaconによるポリッシングと、現在のNanoporeリードに対して適切ではない--nano-rawオプションがデフォルトで使用されているためと考えられます。--racon 0 --opts '-i1' --nanohqを使用するとFlye単体と同等の結果となりました。

MAECIはAutocycler以外で唯一のコンセンサスアセンブリツールで、Canu、Flye、wtdbg2の3つのアセンブラーを統合しています。ただし、今回の実験ではFlyeやCanuよりも高いアセンブリエラー率を示しました。

コンセンサスアプローチは、一貫して単体ツールで発生する小規模なエラーや構造的な不確実性を低減します。前身のTrycyclerは強固なフレームワークを提供していましたが、手動調整を必要とする点が制限的でした。この点において、Autocyclerはコンセンサスアセンブリの利点を完全に自動化されたワークフローにもたらし、正確な細菌ゲノムアセンブリを大規模に実現します。

Autocyclerのインストール

https://github.com/rrwick/Autocycler/wiki/Software-requirements-and-installation

Autocyclerは以下のステップで処理が進みます。

  1. Autocycler subsample
  2. Generating input assemblies
  3. Autocycler compress
  4. Autocycler cluster
  5. Autocycler trim
  6. Autocycler resolve
  7. Autocycler combine

ただ、Autocycler自体にはアセンブリツールが含まれていません。そのため、ステップ2ではAutocyclerが提供してくれているHelperスクリプトを使用して、複数のアセンブラでアセンブリゲノム配列を作成する必要があります。

condaを使用することでまとめてインストールすることができます。

bash
# Miniforge3の導入
wget "https://github.com/conda-forge/miniforge/releases/download/25.3.0-3/Miniforge3-25.3.0-3-Linux-x86_64.sh"
bash Miniforge3-25.3.0-3-Linux-x86_64.sh -b -p $PWD/miniforge3

# activate
. ${PWD}/miniforge3/etc/profile.d/conda.sh
. ${PWD}/miniforge3/etc/profile.d/mamba.sh

conda activate

# Create env for autocycler
git clone https://github.com/rrwick/Autocycler.git
cd Autocycler/scripts
mamba env create --file environment.yml --name autocycler
# activate
mamba activate autocycler

# copy assembly helper scripts into conda env
cp *.py *.sh "$CONDA_PREFIX"/bin
# download plassembler_db
plassembler download -d "$CONDA_PREFIX"/plassembler_db

Autocyclerの公式Wikiで報告されているようなエラーは、今回のインストールでは確認されませんでした。

USAGE Autocycler
help
     /\        | |                      | |          
    /  \  _   _| |_ ___   ___ _   _  ___| | ___ _ __ 
   / /\ \| | | | __/ _ \ / __| | | |/ __| |/ _ \ '__|
  / ____ \ |_| | || (_) | (__| |_| | (__| |  __/ |   
 /_/    \_\__,_|\__\___/ \___|\__, |\___|_|\___|_|   
                               __/ |                 
                              |___/

a tool for generating consensus bacterial genome assemblies
Documenation: https://github.com/rrwick/Autocycler/wiki

Usage: autocycler <COMMAND>

Commands:
  clean       manual manipulation of the final consensus assebly graph
  cluster     cluster contigs in the unitig graph based on similarity
  combine     combine Autocycler GFAs into one assembly
  compress    compress input contigs into a unitig graph
  decompress  decompress contigs from a unitig graph
  dotplot     generate an all-vs-all dotplot from a unitig graph
  gfa2fasta   convert an Autocycler GFA file to FASTA format
  resolve     resolve repeats in the the unitig graph
  subsample   subsample a long-read set
  table       create TSV line from YAML files
  trim        trim contigs in a cluster

Options:
  -h, --help     Print help
  -V, --version  Print version

Fastqの取得

今回、論文ではNaoporeの単離ゲノムシーケンスデータが利用されていました。

Nanoporeのデータ:
PromethION 2 SoloでSQKRBK114.96を使ってシーケンスされたものです。ベースコールはDorado v0.9.5 10 をsup@v5.0.0で実施しています。

以下の種のゲノムシーケンスデータです。

Species Biosample ONT reads Illumina reads sequence size (bp)
Enterobacter hormaechei SAMN31246718 SRR32486094 SRR21866175 4958409, 343098, 80745, 2495
Klebsiella pneumoniae SAMN20033487 SRR32486200 SRR15097896 5369557, 330109, 285776, 3514, 1240
Listeria innocua SAMN46906078 SRR32486273 SRR32484996 2883001, 89544
Providencia rettgeri SAMN46905808 SRR32486258 SRR32458788 4318567, 147239
Shigella flexneri SAMN10100853 SRR32486158 SRR7885912 4583932, 222364, 6790, 5057, 3835, 3181, 2088, 1240

以下のfastqファイルのうち今回はSRR32486273_Listeria_ONT_1を使用します。

mkdir fastq_nanopore

aria2c -x 16 -ctrue ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR324/000/SRR32486200/SRR32486200_1.fastq.gz -o fastq_nanopore/SRR32486200_Klebsiella_ONT_1.fastq.gz
aria2c -x 16 -ctrue ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR324/058/SRR32486158/SRR32486158_1.fastq.gz -o fastq_nanopore/SRR32486158_Escherichia_Shigella_ONT_1.fastq.gz
aria2c -x 16 -ctrue ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR324/073/SRR32486273/SRR32486273_1.fastq.gz -o fastq_nanopore/SRR32486273_Listeria_ONT_1.fastq.gz
aria2c -x 16 -ctrue ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR324/094/SRR32486094/SRR32486094_1.fastq.gz -o fastq_nanopore/SRR32486094_Entobacter_ONT_1.fastq.gz
aria2c -x 16 -ctrue ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR324/058/SRR32486258/SRR32486258_1.fastq.gz -o fastq_nanopore/SRR32486258_Enterobacteriales_ONT_1.fastq.gz

Autocyclerの実行

https://github.com/rrwick/Autocycler/wiki/Fully-automated-assembly

1. Autocycler subsample

mkdir Results
reads="fastq_nanopore/SRR32486273_Listeria_ONT_1.fastq.gz"
threads=32
genome_size=$(genome_size_raven.sh "$reads" "$threads")  # can set this manually if you know the value
ravan output
[raven::] loaded 565098 sequences 29.629916s
[raven::Graph::Construct] minimized 0 - 565098 / 565098 28.270089s
[raven::Graph::Construct] mapped sequences 1124.414242s
[raven::Graph::Construct] annotated piles 0.659261s
[raven::Graph::Construct] removed contained sequences 1.675777s
[raven::Graph::Construct] removed chimeric sequences 0.022267s
[raven::Graph::Construct] minimized 0 - 685 / 685 1.949135s
[raven::Graph::Construct] mapped valid sequences 0.317903s
[raven::Graph::Construct] updated overlaps 0.003249s
[raven::Graph::Construct] removed false overlaps 0.019420s
[raven::Graph::Construct] stored 916 nodes 0.065157s
[raven::Graph::Construct] stored 5536 edges 0.000000s
[raven::Graph::Construct] 1157.803549s
[raven::Graph::Assemble] removed transitive edges 0.005913s
[raven::Graph::Assemble] removed tips and bubbles 0.002364s
[raven::Graph::Assemble] removed long edges 0.020054s
[raven::Graph::Assemble] 0.189339s
[racon::Polisher::Polish] minimized 0 - 5 / 5 0.164475s
[racon::Polisher::Polish] mapped sequences 9.487565s
[racon::Polisher::Polish] found 558916 overlaps 9.709817s
[racon::Polisher::Polish] reverse complemented sequences 0.005071s
[racon::Polisher::Polish] aligned 558916 / 558916 overlaps [================] 35.659629s
[racon::Polisher::Polish] prepared 5941 window placeholders 2.370365s
[racon::Polisher::Polish] called consensus for 5941 / 5941 windows [================] 81.906086s
[racon::Polisher::Polish] 129.715118s
[racon::Polisher::Polish] minimized 0 - 5 / 5 0.141824s
[racon::Polisher::Polish] mapped sequences 10.714563s
[racon::Polisher::Polish] found 560039 overlaps 10.915791s
[racon::Polisher::Polish] reverse complemented sequences 0.002922s
[racon::Polisher::Polish] aligned 560039 / 560039 overlaps [================] 30.580647s
[racon::Polisher::Polish] prepared 5965 window placeholders 2.379284s
[racon::Polisher::Polish] called consensus for 5965 / 5965 windows [================] 81.306597s
[racon::Polisher::Polish] 125.240930s
[raven::] 1442.865148s
autocycler subsample
# Step 1: subsample the long-read set into multiple files
autocycler subsample --reads "$reads" --out_dir Results/subsampled_reads --genome_size "$genome_size"
autocycler subsample 出力
出力
Starting autocycler subsample (2025-05-18 19:39:46)
    This command subsamples a long-read set into subsets that are maximally independent from each other.

Settings:
  --reads fastq_nanopore/SRR32486273_Listeria_ONT_1.fastq.gz
  --out_dir Results/subsampled_reads
  --genome_size 2980984
  --count 4
  --min_read_depth 25
  --seed 0

Input FASTQ:
  Read count: 565098
  Read bases: 3101717815
  Read N50 length: 8793 bp


Calculating subset size (2025-05-18 19:40:03)
    Autocycler will now calculate the number of reads to put in each subset.

Total read depth: 1040.5×
Mean read length: 5489 bp

Calculating subset sizes:
  subset_depth = 25 * log_2(4 * total_depth / 25) / 2
               = 92.2x
  reads per subset: 50096


Subsetting reads (2025-05-18 19:40:03)
    The reads are now shuffled and grouped into subset files.

subset 1:
  reads 1-50096
  Results/subsampled_reads/sample_01.fastq

subset 2:
  reads 141276-191371
  Results/subsampled_reads/sample_02.fastq

subset 3:
  reads 282550-332645
  Results/subsampled_reads/sample_03.fastq

subset 4:
  reads 423825-473920
  Results/subsampled_reads/sample_04.fastq


Finished! (2025-05-18 19:40:22)
    You can now assemble each of the subsampled read sets to produce a set of assemblies for input into Autocycler compress.

2. Generating input assemblies

アセンブル:

denovo assembly
# Step 2: assemble each subsampled file
mkdir Results/assemblies
for assembler in canu flye miniasm necat nextdenovo raven; do
    for i in 01 02 03 04; do
        "$assembler".sh Results/subsampled_reads/sample_"$i".fastq Results/assemblies/"$assembler"_"$i" "$threads" "$genome_size"
    done
done

# Optional step: remove the subsampled reads to save space
# rm subsampled_reads/*.fastq

3. Autocycler compress

# Step 3: compress the input assemblies into a unitig graph
autocycler compress -i Results/assemblies -a Results/autocycler_out -t $threads
出力
Compressed unitig graph: Results/autocycler_out/input_assemblies.gfa
Input assembly stats:    Results/autocycler_out/input_assemblies.yaml

4. Autocycler cluster

# Step 4: cluster the input contigs into putative genomic sequences
autocycler cluster -a Results/autocycler_out
出力
出力
Starting autocycler cluster (2025-05-18 22:52:54)
    This command takes a unitig graph (made by autocycler compress) and clusters the sequences based on their similarity. Ideally, each cluster will then contain sequences which can be
combined into a consensus.

Settings:
  --autocycler_dir Results/autocycler_out
  --cutoff 0.2
  --min_assemblies 6 (automatically set)
  --max_contigs 25


Pairwise distances (2025-05-18 22:52:54)
    Every pairwise distance between contigs is calculated based on the similarity of their paths through the graph.

53 sequences, 2809 total pairwise distances

Saving distance matrix:
  Results/autocycler_out/clustering/pairwise_distances.phylip


Clustering sequences (2025-05-18 22:52:54)
    Contigs are organise into a tree using UPGMA. Then clusters are defined from the tree using the distance cutoff.

Saving clustering tree:
  Results/autocycler_out/clustering/clustering.newick

Cluster 001:
  canu_01.fasta tig00000001 (2930317 bp)
  canu_02.fasta tig00000005 (2882997 bp)
  canu_03.fasta tig00000001 (2939366 bp)
  canu_04.fasta tig00000004 (2882999 bp)
  flye_01.fasta contig_2 (2882999 bp)
  flye_02.fasta contig_1 (2883001 bp)
  flye_03.fasta contig_1 (2882999 bp)
  flye_04.fasta contig_3 (2882999 bp)
  miniasm_01.fasta utg000001c (2882988 bp)
  miniasm_02.fasta utg000001c (2882974 bp)
  miniasm_03.fasta utg000001c (2882979 bp)
  miniasm_04.fasta utg000001c (2882978 bp)
  necat_01.fasta bctg00000000 (2883118 bp)
  necat_02.fasta bctg00000000 (2883004 bp)
  necat_03.fasta bctg00000000 (2883038 bp)
  necat_04.fasta bctg00000000 (2878876 bp)
  nextdenovo_01.fasta ctg000010 (2909158 bp)
  nextdenovo_02.fasta ctg000000 (2902302 bp)
  nextdenovo_03.fasta ctg000010 (2906003 bp)
  nextdenovo_04.fasta ctg000010 (2908102 bp)
  raven_01.fasta Utg742 (2881173 bp)
  raven_02.fasta Utg738 (2883326 bp)
  raven_03.fasta Utg792 (2882994 bp)
  raven_04.fasta Utg738 (2882436 bp)
  cluster distance: 0.001606
  passed QC

Cluster 003:
  canu_01.fasta tig00000003 (89544 bp)
  canu_02.fasta tig00000004 (89544 bp)
  canu_03.fasta tig00000002 (89544 bp)
  canu_04.fasta tig00000003 (89544 bp)
  flye_01.fasta contig_3 (89544 bp)
  flye_02.fasta contig_3 (89544 bp)
  flye_03.fasta contig_3 (89544 bp)
  flye_04.fasta contig_1 (89544 bp)
  miniasm_01.fasta utg000002c (89543 bp)
  miniasm_02.fasta utg000002c (89544 bp)
  miniasm_03.fasta utg000002c (89544 bp)
  miniasm_04.fasta utg000002c (89544 bp)
  necat_01.fasta bctg00000001 (134279 bp)
  necat_02.fasta bctg00000001 (89544 bp)
  necat_03.fasta bctg00000001 (89547 bp)
  necat_04.fasta bctg00000001 (89543 bp)
  nextdenovo_01.fasta ctg000000 (118089 bp)
  nextdenovo_04.fasta ctg000000 (110168 bp)
  raven_01.fasta Utg744 (89539 bp)
  raven_02.fasta Utg740 (89540 bp)
  raven_03.fasta Utg794 (89538 bp)
  raven_04.fasta Utg740 (89540 bp)
  cluster distance: 0.001783
  passed QC

Cluster 002:
  nextdenovo_03.fasta ctg000000 (116983 bp)
  failed QC: present in too few assemblies

Cluster 004:
  canu_04.fasta tig00000005 (54089 bp)
  failed QC: present in too few assemblies
  failed QC: contained within cluster 1

Cluster 005:
  canu_02.fasta tig00000006 (36586 bp)
  failed QC: present in too few assemblies
  failed QC: contained within cluster 1

Cluster 006:
  flye_01.fasta contig_1 (6020 bp)
  flye_02.fasta contig_2 (5416 bp)
  flye_03.fasta contig_2 (5194 bp)
  flye_04.fasta contig_2 (6276 bp)
  cluster distance: 0.153845
  failed QC: present in too few assemblies
  failed QC: contained within cluster 1


Finished! (2025-05-18 22:52:54)
    You can now run autocycler trim on each cluster. If you want to manually inspect the clustering, you can view the following files.

Pairwise distances:         Results/autocycler_out/clustering/pairwise_distances.phylip
Clustering tree (Newick):   Results/autocycler_out/clustering/clustering.newick
Clustering tree (metadata): Results/autocycler_out/clustering/clustering.tsv

5. Autocycler trim & 6. Autocycler resolve

# Steps 5 and 6: trim and resolve each QC-pass cluster
for c in Results/autocycler_out/clustering/qc_pass/cluster_*; do
    autocycler trim -c "$c"
    autocycler resolve -c "$c"
done
出力
Finished! (2025-05-18 22:54:44)
Final consensus graph: Results/autocycler_out/clustering/qc_pass/cluster_003/5_final.gfa

7. Autocycler combine

# Step 7: combine resolved clusters into a final assembly
autocycler combine -a Results/autocycler_out -i Results/autocycler_out/clustering/qc_pass/cluster_*/5_final.gfa
出力
Finished! (2025-05-18 22:56:00)
Combined graph: Results/autocycler_out/consensus_assembly.gfa
Combined fasta: Results/autocycler_out/consensus_assembly.fasta

Consensus assembly is fully resolved 😄

Autocyclerの評価方法

Autocyclerのcombineの出力

combineコマンドはすべてのクラスターを結合して最終的なアセンブリを生成します。各クラスターが単一のcontigにまとめられた場合、Consensus assembly is fully resolved 😄と出力されます。一方、1つ以上のクラスターが単一のシーケンスとしてまとめられなかった場合、One or more clusters failed to fully resolve 😟と出力されます。

ここで実行結果も判断することができます。

Bandageでアセンブリグラフを確認

Bandage/BandageNGでアセンブリグラフを視覚的に確認することで解析結果を確認することができます。

BandageNGの導入
git clone https://github.com/asl/BandageNG.git
cd BandageNG

mkdir build && cd build
cmake ..
make

# 起動
./BandageNG

BandageNGを起動したら、consensus_assembly.gfaを読み込ませて、Draw graphをクリックします。以下のようにゲノム配列とプラスミド配列と思われるものが環状の状態になっていれば問題ないです。

まとめ

本記事では、複数のロングリードアセンブリを統合して高品質なコンセンサス配列を自動生成するツール、Autocycler の概要と使用方法について解説しました。Trycyclerを基盤とし、手動操作を最小限に抑えた設計により、バクテリアゲノムのアセンブリを効率的に実現できる点が大きな特徴です。

また、複数のアセンブリツールを組み合わせることで、単一のツールでは得られなかった高精度な結果が得られる可能性も示されました。特に、Autocyclerによって自動的に選定されたアセンブリの組み合わせは、個別の評価を行う手間を省きつつ、精度の高いコンセンサス配列を導き出します。

導入にあたっては若干の依存関係の調整が必要となるものの、公式Wikiの手順に従えば、大きなトラブルなく環境構築が可能です。また、実際に提供されているデモデータを用いることで、動作確認も容易に行えます。

Autocyclerは、ロングリードを用いた微生物ゲノム解析における高品質アセンブリ生成の有力な選択肢となるツールです。今後の研究において、ゲノムアセンブリの自動化や品質向上を図る際の有効な手段として活用されることが期待されます。

Discussion