Long readのアセンブラとsmall plasmidのアセンブル
はじめに
細菌ゲノムのLong readを用いたアセンブラでアセンブル結果（特にsmall plasmid）に注意が要るかもしれないケースです。
Klebsiella pneumoniaeのゲノムのアセンブルでの例です。
動作確認 (Hardware)
-MacBook Pro
-チップ Apple M1 (Rosetta2 インストール済み)
-macOS Sonoma 14.7.5
-メモリ 16GB
動作確認 (Software)
-sra-tools 3.2.1
-pigz 2.8
-NanoPlot 1.44.1
-filtlong 0.2.1
-fastp 0.25.0
-rasusa 2.1.0
-flye 2.9.6-b1802
-miniasm 0.3
-minimap2 2.29
-minipolish 0.1.3
-any2fasta 0.4.2
-unicycler 0.5.1
-plassembler 1.6.2
-BandageNG 2025.5.1
-abricate 1.0.1
-samtools 1.21
-dnaapler 1.2.0
-bakta 1.11
-pyGenomeViz 1.6.1
-jupyter notebook（bashカーネル インストール済み）
手順
1. データの用意
1-1. 使用するデータ
Microbial Genomicsに掲載の、こちらの文献を例にとります。
IS26-mediated plasmid reshuffling results in convergence of toxin–antitoxin systems but loss of resistance genes in XDR Klebsiella pneumoniae from a chronic infection
こちらの文献中のMRSN 546052株のデータを使用します。
文献によるとゲノムサイズがそれぞれおよそ283.4Kbp、103.7Kbp、70.8Kbp、6.1Kbp、2.1Kbpの5つのプラスミドを保有する株で、6.1Kbpの小さいプラスミドに薬剤耐性遺伝子（blaOXA-232）を持っています。
NCBIのサイトにアクセスして、文献に記載されているBioProjectの番号（PRJNA830444）で検索します。
検索結果のページから「Genomes」の欄内の「SRA」へのリンクをクリックします。
Search resultsに「WGS of MRSN546052: Klebsiella pneumoniae」のPacbioとIlluminaのリードのデータへのリンクがあります。
リンク先のページでSRA Runのアクセッション番号がわかるので、sra-toolsを使用してリードのfastqファイルをダウンロードします。
1-2. 作業ディレクトリの作成
今回の例では任意の場所に「MRSN546052」という名前でディレクトリを作成して移動しておきます。
mkdir path_to/MRSN546052
cd path_to/MRSN546052
1-3. sra-toolsのインストール
conda（mamba）を使って仮想環境を作成しました。
その仮想環境内に入った後 sra-tools をインストールします。
（※自分の環境ではRosettaを使用しています、ARM64での動作は未検証）
mamba create -n sra-tools
conda activate sra-tools
mamba install bioconda::sra-tools
1-4. データをダウンロード
sra-toolsの仮想環境に入った状態でfasterq-dumpを使用してリードのデータをダウンロードします。
オプションの確認
fasterq-dump --help
実行
#Pacbioのリードをダウンロード
fasterq-dump SRR18887187 -p
#Illuminaのペアエンドリードをダウンロード
fasterq-dump SRR18887193 -p --split-files
conda deactivate #仮想環境から出る
- -p ：show progress
- --split-files：write reads into different files
（optional）pigzを使用してダウンロードしたファイルを圧縮しておきます。
#pigzのインストール（MacでHomebrew導入後）
brew install pigz
#ダウンロードしたファイルの圧縮
pigz SRR18887187.fastq
pigz SRR18887193_1.fastq
pigz SRR18887193_2.fastq
圧縮前後のサイズ比較
圧縮前のサイズ
圧縮後のサイズ
今回は圧縮後のデータを使用していきます。
変数を設定して代入しておきます。
genomesize=5500000
long=SRR18887187.fastq.gz
short1=SRR18887193_1.fastq.gz
short2=SRR18887193_2.fastq.gz
threads=8
2. 前処理（QC）
2-1. ロングリードのデータの確認
NanoPlotを使用してロングリードのデータを確認します。
GitHub
#仮想環境の作成とインストール
mamba create -n nanoplot
conda activate nanoplot
mamba install bioconda::nanoplot
オプションの確認
NanoPlot --help
実行
NanoPlot --fastq $long --loglength -t $threads -o NanoPlot_out
conda deactivate #仮想環境から出る
- --loglength：Additionally show logarithmic scaling of lengths in plots
NanoPlot結果
statsのテキストファイル
LengthvsQualityScatterPlot_loglength_dot のグラフ
Non_weightedLogTransformed_HistogramReadlength のグラフ
2-2. ロングリードのトリミング
Filtlongを使用してロングリードのトリミングを実行します。
GitHub
#仮想環境の作成とインストール
mamba create -n filtlong
conda activate filtlong
mamba install bioconda::filtlong
オプションの確認
filtlong --help
実行
filtlong --min_length 1000 --min_mean_q 9 $long > filtlong_reads_1000.fastq
conda deactivate #仮想環境から出る
- --min_length ：minimum length threshold
- --min_mean_q：minimum mean quality threshold
Filtlong実行後
2-3. ショートリードのトリミング
fastpを使用してショートリードのトリミングを実行します。
GitHub
#仮想環境の作成とインストール
mamba create -n fastp
conda activate fastp
mamba install bioconda::fastp
オプションの確認
fastp --help
実行
fastp -i $short1 -I $short2 \
-o short_paired_1.fastq.gz -O short_paired_2.fastq.gz \
-r -q 20 -t 5 -T 5 -w $threads -h report.html
conda deactivate #仮想環境から出る
- -r：move a sliding window from front to tail, if meet one window with mean quality < threshold, drop the bases in the window and the right part, and then stop
- -q：the quality value that a base is qualified. Default 15 means phred quality >=Q15 is qualified
- -t：trimming how many bases in tail for read1, default is 0
- -T：trimming how many bases in tail for read2. If it's not specified, it will follow read1's settings
- -w：worker thread number, default is 3
- -h：the html format report file name
report.html（一部）
2-4. ロングリードのサブサンプリング
rasusaを使用して平均カバレッジが100倍になるようにリードをランダムにサブサンプリングします。
GitHub
#仮想環境の作成とインストール
mamba create -n rasusa
conda activate rasusa
mamba install bioconda::rasusa
オプションの確認
rasusa reads --help
実行
rasusa reads -o filtlong_reads_1000_rasusa100.fastq.gz \
-g $genomesize -c 100 -s 42 -O g filtlong_reads_1000.fastq
conda deactivate #仮想環境から出る
- -g：Genome size to calculate coverage with respect to. e.g., 4.3kb, 7Tb, 9000, 4.1MB
- -c：The desired depth of coverage to subsample the reads to
- -s：Random seed to use
- -O：u: uncompressed; b: Bzip2; g: Gzip; l: Lzma; x: Xz (Lzma); z: Zstd
gzip圧縮されたファイル（filtlong_reads_1000_rasusa100.fastq.gz）が出力されます。
3. ロングリードによるアセンブル
ロングリードのアセンブラとして、
Benchmarking of long-read assemblers for prokaryote whole genome sequencing
こちらの文献を参考にして、FlyeとMiniasm/Minipolish の2つのプログラムを使用してみます。
3-1. Flyeによるロングリードのアセンブル
Flyeを使用してロングリードのアセンブルを実行します。
GitHub
#仮想環境の作成とインストール
mamba create -n flye
conda activate flye
mamba install bioconda::flye
オプションの確認
flye --help
実行
timeコマンドをつけて実行時間も測定しながら実行します。
time flye --pacbio-raw filtlong_reads_1000_rasusa100.fastq.gz \
--out-dir flye_out_filtlong \
--genome-size $genomesize --threads $threads \
conda deactivate #仮想環境から出る
- --pacbio-raw：PacBio regular CLR reads (<20% error)
アセンブル結果
| Flye | ||
|---|---|---|
| total length | 5801303 | |
| contigs | 4 | |
| contig size (circular) | 5342114 (True) | contig_1 |
| 284676 (True) | contig_3 | |
| 103686 (True) | contig_2 | |
| 70827 (True) | contig_4 |
BandageNGでアセンブリグラフを確認しました。
アセンブリグラフ
3-2. Miniasm/Minipolishによるロングリードのアセンブル
Miniasm/Minipolishを使用してロングリードのアセンブルを実行します。
GitHub (miniasm)
GitHub (Minipolish)
#仮想環境の作成とインストール
mamba create -n miniasm_minipolish
conda activate miniasm_minipolish
mamba install bioconda::minipolish
mamba install bioconda::miniasm
mamba install bioconda::any2fasta
Minipolishのリポジトリにあるスクリプト（miniasm_and_minipolish.sh）を利用するため、任意のディレクトリにMinipolishのリポジトリをcloneします。
git clone https://github.com/rrwick/Minipolish.git
実行
timeコマンドをつけて実行時間も測定しながら実行します。
time bash path_to/Minipolish/miniasm_and_minipolish.sh \
filtlong_reads_1000_rasusa100.fastq.gz $threads > miniasm_filtlong.gfa \
&& any2fasta miniasm_filtlong.gfa > miniasm_filtlong.fasta
conda deactivate #仮想環境から出る
アセンブル結果
| Miniasm/Minipolish | ||
|---|---|---|
| total length | 5802720 | |
| contigs | 4 | |
| contig size (circular) | 5343495 (True) | utg000002c |
| 284721 (True) | utg000003c | |
| 103684 (True) | utg000001c | |
| 70820 (True) | utg000004c |
アセンブリグラフ
4. Unicyclerによるハイブリッドアセンブル
Unicyclerを使用してshort-read-firstのハイブリッドアセンブルを実行します。
GitHub
#仮想環境の作成とインストール
mamba create -n unicycler
conda activate unicycler
mamba install bioconda::unicycler
オプションの確認
unicycler --help_all
実行
timeコマンドをつけて実行時間も測定しながら実行します。
time unicycler -1 short_paired_1.fastq.gz -2 short_paired_2.fastq.gz \
-l filtlong_reads_1000_rasusa100.fastq.gz -o unicycler_out_filtlong \
--keep 3 -t $threads --kmer_count 5
conda deactivate #仮想環境から出る
- --keep：Level of file retention
3 = keep all temp files and save all graphs (for debugging) - --kmer_count：Number of k-mer steps to use in SPAdes assembly (default: 8)
アセンブル結果
| Unicycler | ||
|---|---|---|
| total length | 5806790 | |
| contigs | 7 | |
| contig size (circular) | 5341111 (True) | 1 |
| 282920 (True) | 2 | |
| 103694 (True) | 3 | |
| 70829 (True) | 4 | |
| 5961 | 5 | |
| 2095 (True) | 6 | |
| 180 | 7 |
アセンブリグラフ
5. Plassemblerによるハイブリッドアセンブル
Plassemblerを使用してアセンブルを実行します。
Plassemblerは、最初にFlye（またはRaven）によるアセンブルを実行した後、閾値（デフォルトでは1Mb）より大きいサイズの配列を染色体ゲノムと判定して、染色体の配列にマッピングしたリードを除いて、残ったリードでUnicyclerを実行してプラスミドのアセンブルの精度をあげる、という仕組みのようです。
#仮想環境の作成とインストール
mamba create -n plassembler
conda activate plassembler
mamba install bioconda::plassembler
オプションの確認
plassembler run --help
任意のディレクトリにデータベースをダウンロードします。
plassembler download -d path_to/database
実行
timeコマンドをつけて実行時間も測定しながら実行します。
事前にFlyeのアセンブル結果があれば、そのディレクトリやファイルを指定することでFlyeによるアセンブルの工程をスキップすることもできるようです。
今回はFlyeによるアセンブルから実行します。
time plassembler run -d path_to/database \
-l filtlong_reads_1000_rasusa100.fastq.gz \
-1 short_paired_1.fastq.gz -2 short_paired_2.fastq.gz -t $threads \
--keep_fastqs --keep_chromosome --pacbio_model pacbio-raw
#出力先のフォルダ名を変更
mv plassembler.output plassembler.output_filtlong
conda deactivate #仮想環境から出る
- --keep_fastqs：Whether you want to keep FASTQ files containing putative plasmid reads and long reads that map to multiple contigs (plasmid and chromosome)
- --keep_chromosome：If you want to keep the chromosome assembly
- --pacbio_model：Pacbio model for Flye. Must be one of pacbio-raw, pacbio-corr or pacbio-hifi. Use pacbio-raw for PacBio regular CLR reads (<20 percent error), pacbio-corr for PacBio reads that were corrected with other methods (<3 percent error) or pacbio-hifi for PacBio HiFi reads (<1 percent error)
アセンブル結果
染色体ゲノム（5342079bp）は前半のFlyeによるアセンブル部分の結果、それ以外の6つのコンティグは後半のUnicyclerによるアセンブル部分の結果です。
| Plassembler | ||
|---|---|---|
| total length | 5808010 | |
| contigs | 7 | |
| contig size (circular) | 5342079 (True) | chromosome |
| 283172 (True) | 1 | |
| 103694 (True) | 2 | |
| 70829 (True) | 3 | |
| 5961 | 4 | |
| 2095 (True) | 5 | |
| 180 | 6 |
アセンブリグラフ（プラスミド）
アセンブル結果のまとめ：
| Flye | Miniasm/Minipolish | Unicycler | Plassembler | |
|---|---|---|---|---|
| total length | 5801303 | 5802720 | 5806790 | 5808010 |
| contigs | 4 | 4 | 7 | 7 |
| contig size (circular) | 5342114 (True) | 5343495 (True) | 5341111 (True) | 5342079 (True) |
| 284676 (True) | 284721 (True) | 282920 (True) | 283172 (True) | |
| 103686 (True) | 103684 (True) | 103694 (True) | 103694 (True) | |
| 70827 (True) | 70820 (True) | 70829 (True) | 70829 (True) | |
| 5961 | 5961 | |||
| 2095 (True) | 2095 (True) | |||
| 180 | 180 |
今回の環境・設定での実行時間の比較：
| Flye | Miniasm/Minipolish | Unicycler | Plassembler | |
|---|---|---|---|---|
| real | 71m40.887s | 8m55.033s | 117m46.019s | 49m1.376s |
| user | 62m49.411s | 59m24.014s | 677m9.123s | 86m47.247s |
| sys | 189m41.692s | 0m54.279s | 15m46.506s | 63m8.274s |
6. データベース検索
アセンブルされた配列を、abricateを使用してResFinder・PlasmidFinderのデータベースと照合します。
GitHub
#仮想環境の作成とインストール
mamba create -n abricate
conda activate abricate
mamba install bioconda::abricate
#ResFinder・PlasmidFinderのデータベースを最新のものにアップデート
abricate-get_db --db resfinder --force
abricate-get_db --db plasmidfinder --force
オプションの確認
abricate --help
実行
#ResFinderの実行例
abricate -db resfinder --minid 90 --mincov 60 path_to/assembly_fasta
#PlasmidFinderの実行例
abricate -db plasmidfinder --minid 95 --mincov 60 path_to/assembly_fasta
#ターミナル上で列を指定・整形して表示させる例
#abricate -db resfinder --minid 90 --mincov 60 path_to/assembly_fasta | \
# cut -f 2,6,7,9,10,11 | column -s $'\t' -t
conda deactivate #仮想環境から出る
- --minid：Minimum DNA %identity
- --mincov：Minimum DNA %coverage
minidおよびmincovの値は、web版のResfinderおよびPlasmidFinderの記事作成時点でのデフォルト値に準じました。
cutコマンドの-f オプションで抽出する列を指定して、columnコマンドでターミナル上で結果を見やすくするために表示の整形をすることもできます（下の出力結果）。
出力結果
ResFinder
Flye (ResFinder) 染色体とプラスミドの配列
Miniasm (ResFinder) 染色体とプラスミドの配列
Unicycler (ResFinder) 染色体とプラスミドの配列
Plassembler (ResFinder) 染色体の配列
Plassembler (ResFinder) プラスミドの配列
PlasmidFinder
Flye (PlasmidFinder) 染色体とプラスミドの配列
Miniasm (PlasmidFinder) 染色体とプラスミドの配列
Unicycler (PlasmidFinder) 染色体とプラスミドの配列
Plassembler (PlasmidFinder) プラスミドの配列
(文献のSupplementary Data)
7. 6141bpのsmall plasmidの検証
UnicyclerとPlassemblerでは5961bpと180bpの2コンティグで環状のグラフ（計 6141bp）を形成していました（複数コンティグのため「incomplete」判定）。
BandageNGで、この6141bpの環状のグラフをドラッグして選択後「Output」のメニューの「Save selected node sequences to FASTA」を選んで、multi-fastaファイルに配列を書き出します。
デフォルトでは「selected_sequences.fasta」の名前で保存されます。
この配列にリードをマッピングさせて、マッピングしたリードをsam/bamファイルから抽出してfastqファイルに変換後、それらのリードを使ってUnicyclerで再度アセンブルを実行します。
plassemblerの仮想環境に入っているminimap2およびsamtoolsを利用しました。
conda activate plassembler
#PacBioのロングリードのマッピング
minimap2 -t $threads -ax map-pb path_to/selected_sequences.fasta \
filtlong_reads_1000.fastq > long.sam
#マッピングされたリードの抽出とbamファイルに変換
samtools view -b -F 4 -@ $threads long.sam > long_mapped.bam
#fastqファイルに変換
samtools fastq -@ $threads long_mapped.bam > long_mapped.fastq
#ショートリードのマッピング
minimap2 -t $threads -ax sr path_to/selected_sequences.fasta \
short_paired_1.fastq.gz short_paired_2.fastq.gz > short.sam
#マッピングされたリードの抽出とbamファイルに変換
samtools view -b -f 2 -@ $threads short.sam > short_mapped.bam
#ペアエンドのfastqファイルに変換
samtools fastq -1 short_mapped_R1.fastq -2 short_mapped_R2.fastq \
short_mapped.bam
conda deactivate #仮想環境から出る
conda activate unicycler
unicycler -1 short_mapped_R1.fastq -2 short_mapped_R2.fastq \
-l long_mapped.fastq -o unicycler_out_mapped \
--keep 3 -t $threads --kmer_count 5
conda deactivate
unicycler.log（一部）とアセンブリグラフ
6141bpのcompleteな配列にアセンブルされました。
得られた配列をabricateを使用してResFinder・PlasmidFinderのデータベースと照合しました。
ResFinder・PlasmidFinder結果
CP006802（今回アセンブルしたプラスミドと同一株）（文献）との比較結果をpyGenomeVizで可視化してみました。
おわりに
-
ロングリードによるアセンブルでは、Flye、Miniasm/Minipolish ともに文献の283.4Kbp、103.7Kbp、70.8Kbpに相当する3つのプラスミドはアセンブルされましたが、6.1Kbpおよび2.1Kbpの、2つのサイズの小さいプラスミドが検出されませんでした。
-
ロングリードによるアセンブルで得られた配列からは、ResFinderによる検索ではblaOXA-232、PlasmidFinderではColKP-3が検出されませんでした。
-
プラスミドのアセンブル結果を軽視・無視できないようなケースでは特にロングリードのアセンブラのアセンブル結果には慎重になる必要があるかもしれません。
-
今回、アセンブル後のPolishingの工程には触れていません。
-
今回使用した各プログラムのオプション・パラメータの値も一例であり、これが最適解とは限らず、プログラムの選択・パラメータの調整などで、より良い結果が得られる余地はあると思います。
参考
-
Benchmarking of long-read assemblers for prokaryote whole genome sequencing
-
Assembling the perfect bacterial genome using Oxford Nanopore and Illumina sequencing
-
Missing & misassembled plasmids · Issue #706 · mikolmogorov/Flye
-
(preprint) Autocycler: long-read consensus assembly for bacterial genomes
Discussion