Long readのアセンブラとsmall plasmidのアセンブル (2)
はじめに
細菌ゲノムのLong readを用いたアセンブラでアセンブル結果（特にsmall plasmid）に注意が要るかもしれないケース、その2です。
Salmonella属のゲノムのアセンブルでの例です。
以前の記事（Klebsiella pneumoniaeのゲノムのアセンブルでのケース）
Long readのアセンブラとsmall plasmidのアセンブル
動作確認 (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
-emboss 6.6.0
-blast 2.16.0
-bakta 1.11
-BRIG 0.95
-jupyter notebook（bashカーネル インストール済み）
手順
1. データの用意
1-1. 使用するデータ
Microbiology Resource Announcements に掲載の、こちらの文献から例にとります。
Complete Genome Sequence of a Colistin-Susceptible Salmonella enterica Serovar Minnesota Strain Harboring mcr-9 on an IncHI2/IncHI2A Plasmid, Isolated from Chicken Meat
こちらの文献中のSalmonella enterica serotype Minnesota SA18578株のデータを使用します。
文献によるとゲノムサイズがそれぞれおよそ352Kbp、129Kbp、3.2Kbp、3.4Kbpの4つのプラスミドを保有する株で、3.2Kbpの小さいプラスミドに薬剤耐性遺伝子（qnrB19）を持っています。
文献にIlluminaのリードとOxford NanoporeのリードのSRA Runのアクセッション番号が記載されているので、sra-toolsを使用してリードのfastqファイルをダウンロードします。
1-2. 作業ディレクトリの作成
今回の例では任意の場所に「SA18578」という名前でディレクトリを作成して移動しておきます。
mkdir path_to/SA18578
cd path_to/SA18578
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
実行
#Oxford Nanoporeのリードをダウンロード
fasterq-dump SRR16069012 -p
#Illuminaのペアエンドリードをダウンロード
fasterq-dump SRR16071937 -p --split-files
conda deactivate #仮想環境から出る
- -p ：show progress
- --split-files：write reads into different files
pigzを使用してダウンロードしたファイルを圧縮しておきます。
#pigzのインストール（MacでHomebrew導入後）
brew install pigz
#ダウンロードしたファイルの圧縮
pigz SRR16069012.fastq
pigz SRR16071937_1.fastq
pigz SRR16071937_2.fastq
変数を設定して代入しておきます。
genomesize=5000000
long=SRR16069012.fastq.gz
short1=SRR16071937_1.fastq.gz
short2=SRR16071937_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. ロングリードのアダプタートリミング
Porechop_ABIを使用してOxford Nanoporeのロングリードのアダプタートリミングを実行します。
GitHub
#仮想環境の作成とインストール
mamba create -n porechop_abi
conda activate porechop_abi
mamba install bioconda::porechop_abi
オプションの確認
porechop_abi --help
実行
porechop_abi --threads $threads -i $long -o porechop.fastq.gz
conda deactivate #仮想環境から出る
2-3. ロングリードのトリミング
Filtlongを使用してロングリードのトリミングを実行します。
GitHub
#仮想環境の作成とインストール
mamba create -n filtlong
conda activate filtlong
mamba install bioconda::filtlong
オプションの確認
filtlong --help
実行
filtlong --min_length 1000 --min_mean_q 9 porechop.fastq.gz > filtlong_reads_1000.fastq
conda deactivate #仮想環境から出る
- --min_length ：minimum length threshold
- --min_mean_q：minimum mean quality threshold
Filtlong実行後
2-4. ショートリードのトリミング
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-5. ロングリードのサブサンプリング
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 --nano-raw filtlong_reads_1000_rasusa100.fastq.gz \
--out-dir flye_out_filtlong\
--genome-size $genomesize --threads $threads \
conda deactivate #仮想環境から出る
- --nano-raw：ONT reads with oldler chemistries, pre R9 Guppy5 (10-20% error)
アセンブル結果
| Flye | ||
|---|---|---|
| total length | 5297389 | |
| contigs | 3 | |
| contig size (circular) | 4816853 (True) | contig_1 |
| 351250 (True) | contig_2 | |
| 129286 (True) | contig_3 |
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 | 5341878 | |
| contigs | 5 | |
| contig size (circular) | 4820600 (True) | utg000002c |
| 351643 (True) | utg000001c | |
| 129401 (True) | utg000003c | |
| 35340 (True) | utg000004c | |
| 4894 | utg000005l |
アセンブリグラフ
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 | 5312971 | |
| contigs | 7 | |
| contig size (circular) | 4826786 (True) | 1 |
| 352079 (True) | 2 | |
| 90683 | 3 | |
| 35441 (True) | 4 | |
| 3372 (True) | 5 | |
| 3191 (True) | 6 | |
| 1419 | 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
#出力先のフォルダ名を変更
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
アセンブル結果
染色体ゲノム（4818110bp）は前半のFlyeによるアセンブル部分の結果、それ以外の6つのコンティグは後半のUnicyclerによるアセンブル部分の結果です。
| Plassembler | ||
|---|---|---|
| total length | 5342856 | |
| contigs | 7 | |
| contig size (circular) | 4818110 (True) | chromosome |
| 352059 (True) | 1 | |
| 129453 (True) | 2 | |
| 35441 (True) | 3 | |
| 3372 (True) | 4 | |
| 3191 (True) | 5 | |
| 1230 | 6 |
アセンブリグラフ（プラスミド）
アセンブル結果のまとめ：
| Flye | Miniasm/Minipolish | Unicycler | Plassembler | |
|---|---|---|---|---|
| total length | 5297389 | 5341878 | 5312971 | 5342856 |
| contigs | 3 | 5 | 7 | 7 |
| contig size (circular) | 4816853 (True) | 4820600 (True) | 4826786 (True) | 4818110 (True) |
| 351250 (True) | 351643 (True) | 352079 (True) | 352059 (True) | |
| 129286 (True) | 129401 (True) | 90683 | 129453 (True) | |
| 35340 (True) | 35441 (True) | 35441 (True) | ||
| 4894 | 3372 (True) | 3372 (True) | ||
| 3191 (True) | 3191 (True) | |||
| 1419 | 1230 |
今回の環境・設定での実行時間の比較：
| Flye | Miniasm/Minipolish | Unicycler | Plassembler | |
|---|---|---|---|---|
| real | 46m43.315s | 11m33.770s | 92m52.260s | 38m15.897s |
| user | 55m3.685s | 78m36.933s | 553m5.823s | 143m4.773s |
| sys | 99m4.738s | 0m38.738s | 7m10.104s | 2m47.964s |
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) プラスミドの配列
7. 3191bpのsmall plasmidの検証
文献の3.2Kbpのプラスミドの配列 pSA18578_3 (CP080516.1)のfastaファイルをダウンロードして、今回アセンブルされた3191bpのプラスミドを調べてみます。
7-1. multi-fastaファイルの分割
EMBOSSのseqretsplitコマンドを使用して、Unicyclerの出力結果のassembly.fastaファイルを個々のfastaファイルに分割します。
#embossの仮想環境の作成とインストール
mamba create -n emboss
conda activate emboss
mamba install bioconda::emboss
#multi-fastaファイルを分割
seqretsplit path_to/unicycler_out_filtlong/assembly.fasta out
conda deactivate #仮想環境から出る
7コンティグのmulti-fastaファイルから、コンティグごとのfastaファイルに分割されました。
3191bpの「6.fasta」を使用します。
7-2. Local BLAST
BLASTをローカル環境で実行します。
BLASTが入っていない場合はインストールしてPATHを通しておきます。
conda(mamba)を使用してインストールする例↓
mamba install bioconda::blast #blastのインストール
which blastn #blastnにPATHが通っているか確認
which makeblastdb #makeblastdbにPATHが通っているか確認
#データベースの作成
makeblastdb -in path_to/pSA18578_3.fasta -dbtype nucl
#blastnの実行
blastn -db path_to/pSA18578_3.fasta \
-query path_to/6.fasta \
-out blast_result.txt
#-outfmt オプションをつけてタブ区切り形式で出力する例
blastn -db path_to/pSA18578_3.fasta \
-query path_to/6.fasta \
-out blast_result_tab.txt -outfmt 7
makeblastdb
- -dbtype：<String, 'nucl', 'prot'> Molecule type of target db
blastn
- -db：<String> BLAST database name
- query：<File_In> Input file name
- -outfmt：<String> alignment view options:
0 = Pairwise,
1 = Query-anchored showing identities,
2 = Query-anchored no identities,
3 = Flat query-anchored showing identities,
4 = Flat query-anchored no identities,
5 = BLAST XML,
6 = Tabular,
7 = Tabular with comment lines,
8 = Seqalign (Text ASN.1),
9 = Seqalign (Binary ASN.1),
10 = Comma-separated values,
11 = BLAST archive (ASN.1),
12 = Seqalign (JSON),
13 = Multiple-file BLAST JSON,
14 = Multiple-file BLAST XML2,
15 = Single-file BLAST JSON,
16 = Single-file BLAST XML2,
17 = Sequence Alignment/Map (SAM),
18 = Organism Report
blast_result.txt
blast_result_tab.txt
文献のプラスミドは3228bpで、今回アセンブルされたプラスミドと37bpの差がありましたが、37塩基ぶんの配列が2箇所にヒットしていたので、それによる違いのようです。
7-3. BRIG
Web版のBLASTでmegablastを実行しました。
BLAST （megablast検索結果の上位ヒット部分）
上位ヒットの株のfastaファイルをダウンロードしてBRIGで描図してみました。
BRIG
おわりに
-
ロングリードによるアセンブルでは、Flye、Miniasm/Minipolish では3.2Kbp、3.4Kbpの、2つのサイズの小さいプラスミドが検出されませんでした。
-
ロングリードによるアセンブルで得られた配列からは、ResFinderによる検索ではqnrB19、PlasmidFinderではCol typeのプラスミドが検出されませんでした。
-
プラスミドのアセンブル結果を軽視・無視できないようなケースでは特にロングリードのアセンブラのアセンブル結果には慎重になる必要があるかもしれません。
-
今回、アセンブル後のPolishingの工程には触れていません。
-
今回使用した各プログラムのオプション・パラメータの値も一例であり、これが最適解とは限らず、プログラムの選択・パラメータの調整などで、より良い結果が得られる余地はあると思います。
- 129Kbpのプラスミドのアセンブル結果の検証に関しては、今回はsmall plasmidのアセンブルを焦点としているため、本稿では扱いませんでした。
参考
-
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