This tutorial will guide you through the typical steps of analyzing read data from an isolate genome using ONT sequencing data and Illumina.
Overview
Time estimation: 6H
Version: main
Last update: 2026-06-17
Questions:Objectives:
How do I basecall ONT data?
What is the difference between FAST5 and POD5 files?
How can we inspect and convert raw Nanopore data formats?
How do I assemble ONT data
How do I map ONT data?
How do I polish my ONT assembly using Illumina data?
How do run a hybrid assembly?
How do I compare assemblies?
How do I annotate my genome?
Understand the physical and electrical principles of Nanopore sequencing.
Explain the transition from early HMM-based basecallers to modern Transformer-based models.
Interpret Phred quality scores and differentiate between Fast, HAC, and SUP models.
Compare the structural differences between HDF5-based FAST5 and Arrow-based POD5.
Use command-line tools to inspect legacy raw data.
Convert legacy FAST5 files into the modern POD5 standard.
Prerequisites
- Please do the linux introduction before this tutorial.
- Basic understanding of the ONT and Illumina sequencing technology.
- We assume you have successfully connected to an instance in the de.NBI cloud with the software pre-installed. Otherwise you will need to install the required tools on your own and make sure you have sufficient resources available.
- Throughout the course we assume you are working on data downloaded to a volume under
/vol/longread/, we create a link~/workdir/to that folder, if you are working somewhere else, adjust the~/workdirlink to that location and all commands should work as outlined in the course.- We also assume that you have a machine with 28 cores available, if not - adjust the commands that specify a certain number of threads / cores accordingly.
First, create a link to /vol/longread (or the folder in which you want to work during the course) and switch to that directory:
ln -s /vol/longread/ ~/workdir
cd ~/workdir
You might need to change the permissions of /vol/longread, for example (in the cloud setup we use for the on-site course) with:
sudo chown ubuntu:ubuntu /vol/longread/
(Adjust accordingly to your setup)
IMPORTANT: Some software is installed within a python virtual environment, you need to activate it with:
source ~/longread/bin/activate
If some tool cannot be executed during this tutorial - make sure the environment is active! Indicated byt (longread) in your commandline.
Next, we download our tutorial dataset and extract it:
cd ~/workdir
wget https://s3.bi.denbi.de/cmg/mgcourses/longread2026/coursedata.tar.gz
tar -xzvf coursedata.tar.gz
Have a quick look at the content of the coursedata folder:
ls -l ~/workdir/coursedata/
It contains the following components:
Unlike traditional sequencing-by-synthesis methods, Oxford Nanopore Technologies (ONT) sequences native DNA or RNA molecules directly. This allows for the observation of full, unbroken fragments of genomic material, providing unprecedented insights into structural variations, repetitive regions, and base modifications.
The journey from a biological molecule to a digital sequence file relies on a highly sophisticated combination of nanoscale biophysics and deep learning algorithms.
.fast5, nowadays .pod5).What's in a Squiggle? Multiple nucleotides (usually 4 to 6) occupy the narrowest part of the pore constriction at any given millisecond. This means the electrical signal does not correspond to a single base at a time, but rather a moving k-mer context. Translating these overlapping signals requires powerful machine learning models.
Basecalling is the computational process of translating this noisy, continuous electrical current into a discrete sequence of characters (A, C, G, T/U).
The software responsible for basecalling has undergone a massive paradigm shift over the last decade, moving from simple statistical frameworks to cutting-edge deep learning architectures optimized for modern graphics processors (GPUs).
Early basecalling algorithms relied on Hidden Markov Models (HMMs) and simple, shallow neural networks.
.fast5 format (based on the HDF5 standard).Guppy revolutionized long-read bioinformatics by transitioning to GPU acceleration and deep learning frameworks.
Dorado is ONT’s current state-of-the-art production basecaller. Written in optimized C++ for peak hardware performance, it replaces Guppy entirely.
.pod5 file format, significantly speeding up file input/output (I/O).To indicate the reliability of each called base, basecallers output a quality score using the standard Phred scale. The Phred quality score (\(Q\)) is logarithmically linked to the base-calling error probability (\(P\)) by the formula:
\[Q = -10 \log_{10}(P)\]What does Q20 mean? If a base is called with a quality score of 20 (Q20), what is the probability that this base call is an error?
Solution Rearranging the formula: \(P = 10^{-\frac{Q}{10}} = 10^{-\frac{20}{10}} = 10^{-2} = 0.01\) A score of Q20 indicates a 1% error probability, meaning the base call has a 99.0% accuracy.
As algorithms shifted from HMMs to advanced Transformers, and chemistry moved from the R9.4.1 pore to the refined geometry of the R10.4.1 pore, accuracy underwent a dramatic transformation. This shift is often described as entering the “Q20 Era.”
| Era / Pore & Software Combination | Typical Raw Q-Score | Base Call Accuracy | Significance |
|---|---|---|---|
| Early R9 + Albacore | Q7 – Q10 | ~80.0% – 90.0% | High error rates; required heavy consensus bioinformatics to fix systematic mistakes. |
| Late R9.4.1 + Guppy (HAC) | Q13 – Q15 | ~95.0% – 96.8% | The standard for several years; highly capable for structural variants but struggled with homopolymers. |
| Modern R10.4.1 + Dorado (HAC) | Q20+ | ≥ 99.0% | The Q20 Milestone: Raw reads match traditional sequencing accuracy, enabling reliable single-nucleotide variant (SNV) calling. |
| Modern R10.4.1 + Dorado (Duplex) | Q30+ | ≥ 99.9% | High-fidelity long reads, ideal for ultra-accurate rare variant detection and de novo assembly. |
Choosing the Right Dorado Model Dorado offers three different model types depending on your computational budget and laboratory goals:
- Fast: Lightweight models optimized for real-time quality control checks during a run.
- High Accuracy (HAC): The balanced standard baseline for most biological applications.
- Super Accuracy (SUP): Compute-intensive models that require heavy GPU power but yield the highest possible raw Q-scores and the most accurate homopolymer resolution.
Before raw electrical signals (“squiggles”) can be basecalled into sequence characters, they must be stored on disk. Oxford Nanopore Technologies (ONT) has utilized two primary formats for this purpose: FAST5 and POD5.
In this practical exercise, you will peek inside a legacy .fast5 file using standard system tools, convert it to the optimized .pod5 format, and view its tabular metadata on the command line.
Ensure you have hdf5-tools (supplying h5ls) and the official pod5 Python package ready in your environment:
# Install the hdf5-tools command line tools if needed
sudo apt install hdf5-tools
# Install the pod5 command line utility if needed
pip install pod5
Because .fast5 files are binary HDF5 containers, standard text-processing commands like cat, head, or less will output unreadable machine code. Instead, we use specialized HDF5 system utilities like h5ls to explore their internal directories.
Run the following command to recursively view the internal structure of a sample file:
cd ~/workdir/
h5ls -r coursedata/raw/rawdata.fast5
Understanding HDF5 Structure Running this reveals a deeply nested tree layout. Each individual sequencing read gets its own path containing specific arrays like
/Raw/Signal(the continuous stream of picoampere measurements) and/channel_id(the metadata tracking the specific pore block). Extracting files out of this nested design in parallel strains compute resources.
Because modern basecallers like Dorado are heavily optimized for Apache Arrow layouts, converting older data assets to .pod5 is a crucial initial preprocessing step.
Execute the format conversion using the command line toolkit:
pod5 convert fast5 coursedata/raw/rawdata.fast5 --output coursedata/raw/rawdata.pod5
Converting multiple files If you are dealing with an entire run directory containing hundreds of legacy files, you can point the tool to a folder instead:
pod5 convert fast5 ./legacy_fast5_dir/ --output ./modern_pod5_dir/
Speeding up processing via Multi-threading The
pod5 converttool splits workloads efficiently. On multi-core systems, you can append the--threadsor-tflag (e.g.,--threads 4) to execute concurrent data writing routines.
Once converted, you can examine the file contents without extracting the heavy signal traces. The pod5 view tool acts as a rapid extraction program, rendering raw file attributes into a clean, human-readable tabular output.
Run the standard view command:
pod5 view coursedata/raw/rawdata.pod5
To slice the data or quickly check specific fields across your run, you can instruct the program to isolate target metrics and pipe them directly into a tab-separated (.tsv) file:
pod5 view coursedata/raw/rawdata.pod5 --include "read_id,channel,num_samples,end_reason" --output summary.tsv
Checking Data Footprints Execute a file size check on your local disk using
ls -lh coursedata/raw/rawdata.*. What differences do you notice between the two formats?Solution You will instantly see that the
.pod5variant is smaller than its.fast5source file. Columnar storage models compress highly repetitive data streams far more tightly while preserving identical raw electrical signal arrays.
Now that we understand the raw data formats and have our datasets ready, we will perform the actual basecalling using Dorado, Oxford Nanopore’s official state-of-the-art basecaller.
Optional: How to install Dorado manually If Dorado is not pre-installed on your workshop instance, you can download, extract, and add the executable to your environment’s
PATHwith these commands:# Download the Linux x64 binary wget https://cdn.oxfordnanoportal.com/software/analysis/dorado-0.9.6-linux-x64.tar.gz # Extract the archive tar -xzvf dorado-0.9.6-linux-x64.tar.gz # Append the bin folder to your PATH environment variable export PATH=$PATH:$(pwd)/dorado-0.9.6-linux-x64/bin
Which dorado version to use? We are using an older version (0.9.6) of dorado that fits to our data. It is the latest version that still supports 4kHz sample rate. Newer versions of dorado (>1.0.0) don’t support 4kHz models/data anymore. For newer data, watch out for the latest version of dorado and make sure to use the latest model since every generation of models improves the basecall accuracy.
Dorado relies on deep-learning neural network models specifically trained on different pore types, sequencing speeds, and accuracy requirements. You can view all officially available models directly via the command line:
dorado download --list
The models are grouped into distinct categories based on their design and use case:
[molecule]_[pore]_[speed]_[accuracy]@[version]. For example, dna_r10.4.1_e8.2_400bps_fast@v4.1.0 tells us it is for DNA, using R10.4.1 pores, running at a translocation speed of 400 bases per second, utilizing the Fast neural network architecture._5mCG_5hmCG or _6mA).herro-v1 or models containing _polish) optimized for post-processing steps or consensus correction. They are used to maximize final sequence accuracy by refining draft assemblies or correction arrays.In addition, dorado models are categorized into three main performance tiers:
Let’s download both the Fast simplex model (which we will run completely) and the SUP simplex model (to test its performance impact):
# Download the Fast model specified for our R10.4.1 chemistry
dorado download --model dna_r10.4.1_e8.2_400bps_fast@v4.1.0
# Download the corresponding Super Accuracy (SUP) model
dorado download --model dna_r10.4.1_e8.2_400bps_sup@v4.1.0
We will now perform basecalling on our raw sample data using the lightweight fast model. By default, Dorado outputs alignments in SAM/BAM format. To force a traditional text-based output, we will pass the --emit-fastq flag.
Run the basecaller on our raw POD5 file/directory:
dorado basecaller dna_r10.4.1_e8.2_400bps_fast@v4.1.0 ~/workdir/coursedata/raw/rawdata.pod5 --emit-fastq > ~/workdir/fast_reads.fastq
Once completed, let’s inspect the first few lines of your freshly generated FASTQ sequence file:
head -n 4 ~/workdir/fast_reads.fastq
Analyzing the FASTQ structure Can you identify the four core lines making up a standard FASTQ record? What do the characters in the fourth line represent?
Solution A standard FASTQ block contains:
@followed by the unique read identifier metadata.- The actual base sequence characters (A, C, G, T).
- A
+separator line.- The ASCII-encoded Phred quality scores (\(Q\)) mapping directly to the error probability of each called base.
Now, let’s observe the massive computational difference when running the high-complexity SUP model. Launch the following command:
dorado basecaller dna_r10.4.1_e8.2_400bps_sup@v4.1.0 ~/workdir/coursedata/raw/rawdata.pod5 --emit-fastq > ~/workdir/sup_reads.fastq
Observing Performance Profiles Notice how the processing speed drops dramatically compared to the Fast run. Without a high-end dedicated graphics processor (GPU), running the Super Accuracy model on large datasets can take a massive amount of time.
Once you have observed this behavior for a few moments, terminate the running process manually by pressing
Ctrl + Cin your terminal window.
For the downstream analysis sections of this course, we want to work with a complete dataset. We have already basecalled several raw data files (located in coursedata/ont/). To make processing streamlined (since many tools only take a single input file), we will concatenate all our prepared Nanopore FASTQ files into a single file using cat.
Execute the following command to merge your ONT reads:
cat ~/workdir/coursedata/ont/*.fastq.gz > ~/workdir/coursedata/ont.fastq.gz
Verifying file concatenation You can verify the success of your merge step by listing the files and evaluating their sizes:
ls -lh ~/workdir/coursedata/ont.fastq.gz ~/workdir/coursedata/ont/*.fastq.gz
While we will use the full pre-prepared datasets for our hands-on exercises, it is important to know how to manage adapter trimming and quality filtering directly during the basecalling stage using Dorado.
Unlike other tools where you must explicitly turn trimming on, Dorado enables trimming of adapters, primers, and barcodes by default for DNA sequencing. You can modify or disable this behavior using specific flags:
--no-trim (or --trim none): Skips all trimming entirely, keeping raw adapter/barcode sequences in your reads.--trim adapters: Restricts trimming to only sequencing adapters, leaving barcodes or primers intact.--min-qscore [value]: Discards any read whose mean Phred quality score falls below your specified threshold (default is 0, meaning no filtering).# Example reference command (Do not run during this workshop):
# This skips default barcode/primer trimming and filters out reads with a mean Q-score below 10:
# dorado basecaller dna_r10.4.1_e8.2_400bps_hac@v5.0.0 input_pod5/ --emit-fastq --trim adapters --min-qscore 10 > filtered_reads.fastq
To evaluate the success of our sequencing run, we must look at quality metrics such as read length distributions and Phred score frequencies. We will use FastQC and NanoPlot to analyze our long-read ONT data, and compare the profile to short-read Illumina data.
First, let’s run FastQC on our combined Nanopore reads:
Optional: How to install fastqc Run this (or follow instructions on the github linked in the Appendix)
sudo apt install fastqc
mkdir ~/workdir/fastqc
fastqc -t 8 ~/workdir/coursedata/ont.fastq.gz -o ~/workdir/fastqc/
To view the generated interactive HTML quality report, open it in a web browser using firefox:
firefox ~/workdir/fastqc/ont_fastqc.html
Since FastQC was originally designed for short reads, it can struggle with the highly variable read lengths of Nanopore data. Therefore, we use NanoPlot, a tool specifically tailored for long-read technologies.
Run NanoPlot to generate a comprehensive quality report:
Optional: How to install NanoPlot Run this (or follow instructions on the github linked in the Appendix)
pip install nanoplot
NanoPlot --fastq ~/workdir/coursedata/ont.fastq.gz -o ~/workdir/nanoplot_ont
NanoPlot creates a dedicated output folder filled with various diagnostic plots (PNG/PDF formats) and a master HTML report summary. Let’s open the main report:
firefox ~/workdir/nanoplot_ont/NanoPlot-report.html
Comparing FastQC and NanoPlot Look at both browser tabs. What key metrics does NanoPlot show that are missing or misrepresented in FastQC?
Solution NanoPlot provides explicit long-read metrics such as N50 read length, active pore counters over time, and read length vs. read quality scatter plots. FastQC often flags long reads as “failed” due to length variations because it expects uniform short-read lengths.
fastplong is an ultra-fast, all-in-one quality control and preprocessing tool explicitly designed for long-read sequencing datasets (Nanopore and PacBio). It serves as the long-read counterpart to the widely used short-read tool fastp.
In a single multi-threaded execution pass, fastplong can perform:
Optional: How to install fastplong Run this (or follow instructions on the github linked in the Appendix)
wget http://opengene.org/fastplong/fastplong chmod a+x ./fastplong sudo mv fastplong /usr/local/bin/
Let’s run fastplong to enforce a minimum length requirement of 1,000 bp (-l 1000) and generate a localized performance dashboard:
fastplong -i ~/workdir/coursedata/ont.fastq.gz \
-o ~/workdir/coursedata/ont_trimmed.fastq.gz \
-l 1000 \
-h ~/workdir/fastplong_report.html
Open the interactive HTML summary generated natively by fastplong:
firefox ~/workdir/fastplong_report.html
To visualize the direct impact of our preprocessing step, let’s run FastQC and NanoPlot a second time, pointing them to our freshly filtered output file:
# Re-run FastQC on trimmed data
fastqc -t 8 ~/workdir/coursedata/ont_trimmed.fastq.gz -o ~/workdir/fastqc/
firefox ~/workdir/fastqc/ont_trimmed_fastqc.html ~/workdir/fastqc/ont_fastqc.html
# Re-run NanoPlot on trimmed data
NanoPlot --fastq ~/workdir/coursedata/ont_trimmed.fastq.gz -o ~/workdir/nanoplot_ont_trimmed/
firefox ~/workdir/nanoplot_ont_trimmed/NanoPlot-report.html ~/workdir/nanoplot_ont/NanoPlot-report.html
Analyzing Post-Filtering Changes Compare the browser tabs for the raw NanoPlot report (
nanoplot_ont) vs. the filtered report (nanoplot_ont_trimmed). What happened to the total read count, the mean read length, and the N50 metric?Solution Because we enforced a minimum length filter of 1,000 bp, all short fragments were eliminated. This causes the total read count to decrease, while both the mean read length and the N50 value increase. The overall mean Phred quality profile also shifts upwards, as short, degraded molecules are successfully purged from the dataset.
Now, let’s run FastQC on our paired-end Illumina data to see the differences in throughput and quality distributions:
fastqc ~/workdir/coursedata/illumina/Barcode11_TSLF_S10_L001_R1_001.fastq.gz \
~/workdir/coursedata/illumina/Barcode11_TSLF_S10_L001_R2_001.fastq.gz \
-t 8 \
-o ~/workdir/fastqc/
Open both resulting reports to inspect the short-read profiles:
firefox ~/workdir/fastqc/Barcode11_TSLF_S10_L001_R1_001_fastqc.html ~/workdir/fastqc/Barcode11_TSLF_S10_L001_R2_001_fastqc.html
To identify variants or verify the structure of our sequenced strain, we need to align our reads to the provided reference genome (reference.fasta). Long reads and short reads require completely different alignment algorithms.
Create a folder to store our mappings first:
mkdir ~/workdir/mappings/
Optional: How to install minimap2 Run this (or follow instructions on the github linked in the Appendix)
curl -L https://github.com/lh3/minimap2/releases/download/v2.31/minimap2-2.31_x64-linux.tar.bz2 | tar -jxvf - export PATH=$PATH:$(pwd)/minimap2-2.31_x64-linux/ chmod a+x ./fastplong sudo mv fastplong /usr/local/bin/
Minimap2 is the standard for aligning long, error-prone reads. We will use the -ax map-ont preset optimized for Oxford Nanopore data.
minimap2 -t 28 -ax map-ont ~/workdir/coursedata/reference.fasta ~/workdir/coursedata/ont_trimmed.fastq.gz > ~/workdir/mappings/ont_mapped.sam
For short reads, we use BWA (Burrows-Wheeler Aligner). Before aligning, BWA requires building an index of the reference genome.
Optional: How to install bwa Run this (or follow instructions on the github linked in the Appendix)
sudo apt install bwa
# Index the reference genome
bwa index ~/workdir/coursedata/reference.fasta
# Align the paired-end short reads using BWA-MEM
bwa mem -t 28 ~/workdir/coursedata/reference.fasta \
~/workdir/coursedata/illumina/Barcode11_TSLF_S10_L001_R1_001.fastq.gz \
~/workdir/coursedata/illumina/Barcode11_TSLF_S10_L001_R2_001.fastq.gz \
> ~/workdir/mappings/illumina_mapped.sam
Our alignment outputs are currently in SAM (Sequence Alignment Map) format, which is human-readable text but highly inefficient for visualization tools. We need to convert them to binary BAM format, sort them by genomic coordinates, and build a coordinate index (.bai) for quick random-access loading in IGV (Integrative Genomics Viewer).
Optional: How to install samtools Run this (or follow instructions on the github linked in the Appendix)
sudo apt install samtools
# Convert SAM to BAM
samtools view -S -b ~/workdir/mappings/ont_mapped.sam > ~/workdir/mappings/ont_mapped.bam
# Sort the BAM file by coordinate position
samtools sort ~/workdir/mappings/ont_mapped.bam -o ~/workdir/mappings/ont_sorted.bam
# Index the sorted BAM file
samtools index ~/workdir/mappings/ont_sorted.bam
# Convert SAM to BAM
samtools view -S -b ~/workdir/mappings/illumina_mapped.sam > ~/workdir/mappings/illumina_mapped.bam
# Sort the BAM file by coordinate position
samtools sort ~/workdir/mappings/illumina_mapped.bam -o ~/workdir/mappings/illumina_sorted.bam
# Index the sorted BAM file
samtools index ~/workdir/mappings/illumina_sorted.bam
Now that both datasets are fully optimized and indexed, we will visually inspect and compare our mapping data to explore how long reads perform against short reads across trickier genomic features.
Optional: How to install IGV Run this: and start IGV via
igv.sh.wget https://data.broadinstitute.org/igv/projects/downloads/2.19/IGV_Linux_2.19.8_WithJava.zip unzip IGV_Linux_2.19.8_WithJava.zip #start IGV_Linux_2.19.8/igv.sh
IGV
~/workdir/coursedata/reference.fasta.~/workdir/ont_sorted.bam~/workdir/illumina_sorted.bamInvestigating Structural Differences Browse through the alignment tracks. What is the main difference in terms of sequencing errors? What specific problem is caused by the ONT errors?
Solution Illumina has less errors and they are usullay substitutions, ONT has a high amount of insertions. Which causes a problem for gene prediction and annotation due to frameshifts.
Genome assembly is the process of piecing together massive amounts of short or long DNA fragments (reads) to reconstruct the original underlying chromosome. Because we do not use a reference genome during de novo assembly, the algorithms rely strictly on sequence overlaps. Two main algorithmic paradigms dominate the field:
Flye is a specialized de novo assembler designed for long, error-prone reads. Instead of building a classic OLC overlap graph (which scales poorly with high read depths), Flye constructs an unpolished repeat graph. It collapses complex genomic repeats into single edges, and then utilizes the long span of individual read paths to accurately untangle and separate those repeat copies.
Optional: How to install Flye Run this (or follow instructions in github):
pip install setuptools git clone https://github.com/fenderglass/Flye cd Flye python setup.py install
SPAdes (St. Petersburg Genome Assembler) is the gold standard for bacterial short-read assemblies. It relies on multi-sized De Bruijn Graphs (combining multiple $k$-mer lengths) to simultaneously optimize specificity and sensitivity, providing robust performance across single-isolate cultures and single-cell sequencing.
Optional: How to install Flye Run this (or follow instructions in github):
wget https://github.com/ablab/spades/releases/download/v4.3.0/SPAdes-4.3.0-Linux.tar.gz tar -xzvf SPAdes-4.3.0-Linux.tar.gz export PATH=PATH:$(pwd)/SPAdes-4.3.0-Linux/bin/
QUAST (Quality Assessment Tool) is an evaluation utility that calculates structural metrics (like contig counts, N50 value, and total length) and identifies misassemblies by aligning your assembled contigs back against a trusted reference genome.
Optional: How to install Quast Run this (or follow instructions in github):
wget https://github.com/ablab/quast/releases/download/quast_5.3.0/quast-5.3.0.tar.gz tar -xzvf quast-5.3.0.tar.gz cd quast-5.3.0 ./setup.py install
Bandage (Bioinformatics Application for Navigating De Novo Assembly Graphs Easily) is a graphical interface utility that reads Graphical Assembly Graph (.gfa) files. It allows you to see how contigs connect to one another, helping you determine whether your bacterial genome successfully assembled into a single closed circular chromosome.
Optional: How to install Bandage Run this (or follow instructions in github):
wget https://github.com/rrwick/Bandage/releases/download/v0.8.1/Bandage_Ubuntu_dynamic_v0_8_1.zip unzip Bandage_Ubuntu_dynamic_v0_8_1.zip sudo mv Bandage /usr/local/bin/ #might be necessary: sudo apt install libqt5svg5
In this section, we will run separate long-read and short-read assembly pipelines, statistically benchmark their outputs against our reference, inspect their connectivity graphs, and align the draft contigs visually.
Because modern Dorado basecalled data achieves exceptional accuracy (entering the Q20 standard), we use Flye’s high-fidelity option (--nano-hq) to generate our draft genome:
flye --nano-hq ~/workdir/coursedata/ont.fastq.gz --out-dir ~/workdir/flye_output --threads 28
Next, we generate a corresponding short-read assembly utilizing our paired-end Illumina datasets:
spades.py -1 ~/workdir/coursedata/illumina/Barcode11_TSLF_S10_L001_R1_001.fastq.gz \
-2 ~/workdir/coursedata/illumina/Barcode11_TSLF_S10_L001_R2_001.fastq.gz \
-o ~/workdir/spades_output --threads 28
We can now run a direct comparative evaluation between both assembly results using our known genome sequence as a reference:
quast.py ~/workdir/flye_output/assembly.fasta \
~/workdir/spades_output/contigs.fasta \
-r ~/workdir/coursedata/reference.fasta \
-o ~/workdir/quast_output
Open the interactive QUAST HTML summary document in your browser to view the benchmark comparison:
firefox ~/workdir/quast_output/report.html
Analyzing Assembly Metrics Look at the metric comparisons in the QUAST report. Which assembly contains fewer total contigs? Which possesses a higher N50 score? What does this tell you about the power of long reads?
Solution Typically, the Flye long-read assembly will result in significantly fewer contigs (often a single continuous contig for a closed bacterial chromosome) and a drastically higher N50 score approaching the true size of the genome. The SPAdes short-read assembly is usually split across multiple fragments because short fragments cannot resolve genomic repeats.
Statistical metrics only tell half the story. We need to look at the assembly graphs to see the structure of our contigs.
Bandage
~/workdir/flye_output/assembly_graph.gfa and click Draw graph.~/workdir/spades_output/assembly_graph_with_scaffolds.gfa and click Draw graph.Interpreting Graph Topologies In the Flye window, you should see a single, beautiful, interconnected closed loop representing the intact circular bacterial chromosome. In contrast, the SPAdes graph will likely display a highly fragmented web of disjointed paths and isolated nodes, highlighting where the short-read assembly broke down at repeat boundaries.
Finally, we want to align our assembled fasta contigs back against the reference genome to visually spot missing structural parts or mismatches in IGV. We use minimap2 with the -ax asm5 preset, which is optimized for aligning highly accurate genome assemblies.
# Map the Flye assembly contigs
minimap2 -t 28 -ax asm5 ~/workdir/coursedata/reference.fasta ~/workdir/flye_output/assembly.fasta > ~/workdir/mappings/flye_vs_ref.sam
samtools view -S -b ~/workdir/mappings/flye_vs_ref.sam | samtools sort -o ~/workdir/mappings/flye_vs_ref_sorted.bam
samtools index ~/workdir/mappings/flye_vs_ref_sorted.bam
# Map the SPAdes assembly contigs
minimap2 -t 28 -ax asm5 ~/workdir/coursedata/reference.fasta ~/workdir/spades_output/scaffolds.fasta > ~/workdir/mappings/spades_vs_ref.sam
samtools view -S -b ~/workdir/mappings/spades_vs_ref.sam | samtools sort -o ~/workdir/mappings/spades_vs_ref_sorted.bam
samtools index ~/workdir/mappings/spades_vs_ref_sorted.bam
~/workdir/coursedata/reference.fasta) is actively loaded.~/workdir/mappings/flye_vs_ref_sorted.bam~/workdir/mappings/spades_vs_ref_sorted.bamEven though modern ONT R10.4.1 chemistry combined with Dorado pushes raw read accuracy into the Q20 (>99%) range, long-read assemblies can still retain minor systematic errors. These errors are most frequently found in homopolymer runs (e.g., long stretches of AAAA), manifesting as small insertions or deletions (indels). To fix these remaining micro-errors, we can perform a process called polishing using highly accurate Illumina short reads.
We will use Polypolish, a short-read polishing tool designed specifically for long-read assemblies.
How Polypolish Avoids False Corrections Traditional polishers take all short-read alignments and use a consensus to alter the assembly. However, in repetitive genomic regions, short reads frequently misalign to the wrong repeat copy, leading the polisher to introduce errors rather than fix them.
Polypolish solves this by examining the alternative alignments for each short read. If a read can map to multiple places in the assembly, Polypolish will only propose a correction if all possible target sites agree on the mismatch. If the mapping is ambiguous, it leaves the sequence untouched, preventing false corrections in repeat boundaries.
To ensure Polypolish operates effectively, we must execute a specific multi-step pipeline:
bwa mem with the -a option. This option forces the aligner to output all possible alignment locations for a read, not just the single best hit.polypolish filter to calculate the expected insert size of read pairs and filter out low-confidence alignments.polypolish polish to correct the assembly using the filtered pileups.Optional: How to install Polypolish Run this (or follow instructions in github):
wget https://github.com/rrwick/Polypolish/releases/download/v0.6.1/polypolish-linux-x86_64-musl-v0.6.1.tar.gz tar -xzvf polypolish-linux-x86_64-musl-v0.6.1.tar.gz sudo mv polypolish /usr/local/bin/
First, let’s build the BWA index of our long-read genome assembly and map both Illumina forward and reverse files completely independently using the required -a flag:
# Index the Flye draft genome
bwa index ~/workdir/flye_output/assembly.fasta
# Create directory for polypolish files
mkdir polypolish
# Map R1 and R2 forward/reverse reads completely independently with the -a flag
bwa mem -t 28 -a ~/workdir/flye_output/assembly.fasta ~/workdir/coursedata/illumina/Barcode11_TSLF_S10_L001_R1_001.fastq.gz > ~/workdir/polypolish/polypolish_r1.sam
bwa mem -t 28 -a ~/workdir/flye_output/assembly.fasta ~/workdir/coursedata/illumina/Barcode11_TSLF_S10_L001_R2_001.fastq.gz > ~/workdir/polypolish/polypolish_r2.sam
Next, we pass our independent raw SAM files into Polypolish’s filtering subcommand. This evaluates read pairing distances to clear away non-specific background mappings:
polypolish filter --in1 ~/workdir/polypolish/polypolish_r1.sam --in2 ~/workdir/polypolish/polypolish_r2.sam --out1 ~/workdir/polypolish/filtered_r1.sam --out2 ~/workdir/polypolish/filtered_r2.sam
Now, we provide the original unpolished Flye assembly along with both freshly filtered alignment tracks to create our refined fasta file:
polypolish polish ~/workdir/flye_output/assembly.fasta ~/workdir/polypolish/filtered_r1.sam ~/workdir/polypolish/filtered_r2.sam > ~/workdir/polypolish/flye_polished.fasta
Instead of assembling long reads first and polishing them later, a hybrid assembly combines both data types simultaneously into a single algorithmic workflow.
We will use the hybrid mode of SPAdes. The SPAdes hybrid approach works as follows:
This approach combines the single-nucleotide accuracy of short reads with the structural spanning power of long reads seamlessly.
Execute the hybrid SPAdes pipeline by supplying both your paired-end short reads and your combined long-read datasets:
spades.py -1 ~/workdir/coursedata/illumina/Barcode11_TSLF_S10_L001_R1_001.fastq.gz \
-2 ~/workdir/coursedata/illumina/Barcode11_TSLF_S10_L001_R2_001.fastq.gz \
--nanopore ~/workdir/coursedata/ont.fastq.gz \
-o ~/workdir/spades_hybrid_output --threads 28
We now have four distinct assembly variants tracking our target genome. Let’s run a final comparative evaluation with QUAST to see how polishing and hybrid strategies alter genome completeness and accuracy metrics.
The 4 assembly variants to evaluate are:
flye_output/assembly.fasta (ONT Long-Reads Only)spades_output/scaffolds.fasta (Illumina Short-Reads Only)flye_polished.fasta (ONT Long-Reads Polished with Short-Reads)spades_hybrid_output/contigs.fasta (Hybrid Co-Assembly)Run QUAST with all four assembly files against the true reference genome:
quast.py ~/workdir/flye_output/assembly.fasta \
~/workdir/spades_output/contigs.fasta \
~/workdir/polypolish/flye_polished.fasta \
~/workdir/spades_hybrid_output/contigs.fasta \
-l "flye,spades,polypolish,hybrid_spades" \
-r ~/workdir/coursedata/reference.fasta \
-t 28 \
-o ~/workdir/quast_final_output
Open the resulting dashboard summary report in your browser:
firefox ~/workdir/quast_final_output/report.html
Evaluating the Impact of Polishing and Hybridization Compare the column profiles of the unpolished Flye assembly vs. the polished Flye assembly. Look at metrics like “mismatches per 100 kbp” or “indels per 100 kbp”. What changes do you observe? How does the Hybrid assembly compare in contig count?
Solution Polishing with Polypolish typically causes a significant drop in the number of indels per 100 kbp compared to raw Flye contigs, which often restores disrupted open reading frames and increases the total number of fully recovered genes. The Hybrid SPAdes assembly often improves in contiguity, but depending on repeat complexity, it may still contain a few more contig fragments than Flye’s completely closed loop structure.
Once you have successfully assembled and polished a bacterial chromosome, it consists simply of a long, uncharacterized string of nucleotides (A, C, G, T). To make this data useful for biological research, you must perform genome annotation. This process involves identifying the structural features of the genome—such as protein-coding sequences (CDS), transfer RNAs (tRNAs), and ribosomal RNAs (rRNAs)—and assigning functional biological identities to them based on sequence similarity to known databases.
In this module, we will compare two popular tools used for this task:
Optional: How to install Prokka It’s quite complicated to install without conda/docker/singularity. Check out the github and use one of those methods.
Optional: How to install Bakta It’s quite complicated to install without conda/docker/singularity. Check out the github and use one of those methods.
We will run both annotators on our polished long-read assembly (flye_polished.fasta) and evaluate how their structural findings and functional naming conventions differ.
Execute Prokka by specifying an output directory and a custom file prefix:
prokka --cpus 28 --outdir ~/workdir/prokka_output --prefix prokka_ont ~/workdir/polypolish/flye_polished.fasta
Don;t worry about the Could not run command: tbl2asn message if it appears. We don’t need the asn file.
Unlike Prokka, Bakta relies on a separate, heavy database containing millions of curated proteins. For this workshop, this database has been pre-staged for you. Run Bakta using the following command:
conda activate bakta
bakta --threads 28 --db ~/bakta_db/db-light --output ~/workdir/bakta_output ~/workdir/polypolish/flye_polished.fasta
conda deactivate
source ~/longread/bin/activate
Both tools generate various standardized outputs, including GFF3, GenBank, and FASTA files. To quickly benchmark their structural predictions, we can review the text-based summary logs (.txt) produced by each pipeline.
Use cat to print out both overview profiles in your terminal:
# View the Prokka summary report
cat ~/workdir/prokka_output/prokka_ont.txt
# View the Bakta summary report
cat ~/workdir/bakta_output/bakta_ont.txt
Analyzing Annotation Discrepancies Look closely at the total counts of Coding Sequences (CDS), tRNAs, and rRNAs in both outputs. Are the numbers identical? If they differ, what could cause one tool to predict more genes than the other?
Solution Even though they use the same underlying software for core gene finding (Prodigal), the total counts often differ slightly. Bakta uses stricter structural filters and a much larger database, allowing it to accurately split overlapping reading frames, filter out false positive predictions, and identify specialized elements like pseudogenes or small non-coding RNAs that Prokka completely misses.
A major difference lies in how specifically proteins are named. Let’s use grep to check how many genes were left uncharacterized as “hypothetical protein” in both annotation suites:
# Count hypothetical proteins in Prokka's GFF output
grep -c "hypothetical protein" ~/workdir/prokka_output/prokka_ont.gff
# Count hypothetical proteins in Bakta's GFF output
grep -c "hypothetical protein" ~/workdir/bakta_output/bakta_ont.gff
Interpreting Naming Quality You will notice that Bakta significantly reduces the fraction of “hypothetical protein” labels compared to Prokka. Thanks to its modern reference integration with UniProt and RefSeq, Bakta can assign definitive, functional gene names to sequences where Prokka could only find vague, outdated family matches.
Once individual genomes are annotated, the next logical milestone is to explore how multiple strains or species relate to one another. For this downstream phase, we shift from localized command-line annotation to web-based comparative genomics using EDGAR (Efficient Database framework for comparative Genome Analyses).
EDGAR is a fully automated high-throughput platform tailored for the deep comparative analysis of prokaryotic genomes. Users upload their fully annotated genome files (such as the .gff or GenBank files generated by Bakta) into public or password-protected private projects.
The underlying pipeline performs intensive all-versus-all sequence alignments across all selected strains. By evaluating BLAST Score Ratios (BSR), EDGAR accurately determines orthology relational paths to delineate specific genomic subsets:
Furthermore, EDGAR calculates average nucleotide identity (ANI) metrics and renders publication-ready visualizations, including Venn diagrams, UpSet plots, and synteny maps mapping gene order conservation across syntenic chromosomal layouts.
Key Points
Basecalling translates raw ionic current disruptions (‘squiggles’) into nucleotide sequences.
The evolution from Albacore to Guppy and Dorado represents a major shift from CPUs to GPU-accelerated deep learning.
Modern R10.4.1 chemistry paired with Dorado pushes accuracies past Q20 (>99%), matching or exceeding short-read standards.
FAST5 uses a hierarchical structure that causes I/O bottlenecks during high-throughput basecalling.
POD5 uses a flat, columnar format (Apache Arrow) that offers faster multi-threaded access and smaller file sizes.
Author(s): Nils Kleinbölting