Bio Applied Ont Processing

ONT Long-Read Data Processing

When to Use

• Basecalling raw Oxford Nanopore POD5/FAST5 signal into sequence + quality (Dorado)
• Choosing a basecalling model tier (fast/hac/sup) or duplex mode for an accuracy/throughput tradeoff
• QC'ing ONT reads with NanoStat/NanoPlot and filtering low-quality/short reads with NanoFilt
• Aligning ONT reads with Minimap2's map-ont preset and interpreting samtools flagstat
• Extracting 5mC/5hmC methylation calls from modification-aware basecalling BAMs with modkit

Version Compatibility

Dorado ≥0.7 (CUDA/Metal), NanoPack2 (NanoStat/NanoPlot/NanoFilt) ≥1.40, minimap2 ≥2.26, samtools ≥1.19, modkit ≥0.2, Python ≥3.10 with numpy/pandas/matplotlib for downstream QC analysis.

Prerequisites

• pip install nanostat nanoplot nanofilt pandas numpy matplotlib
• Dorado is a standalone GPU binary, not a pip package: https://github.com/nanoporetech/dorado
• minimap2 and samtools (conda/bioconda) for alignment; modkit (ONT) for methylation pileup
• Prior concepts: bio-sequence-io-fastq-quality, bio-alignment-files-sam-bam-basics

ONT Technology Overview

ONT sequences DNA by threading a strand through a protein nanopore; a constant voltage drives ionic current, and each k-mer occupying the constriction disrupts current in a characteristic way. The resulting picoampere time-series ("squiggle") encodes sequence. R9.4.1 chemistry uses a 5-mer sensing region; R10.4.1's dual-reader pore uses a 9-mer window, improving homopolymer resolution and raw accuracy.

POD5 (Apache Arrow-based, columnar) is the current signal format, replacing legacy FAST5 (HDF5-based). Dorado reads POD5 natively; convert old data with pod5 convert fast5.

Basecalling and QC

Goal: Turn raw POD5 signal into filtered, QC'd FASTQ/BAM ready for alignment. Approach: Run Dorado with a model tier matched to the use case (fast for screening, hac for routine genomics, sup/duplex for clinical-grade accuracy), then gate reads with NanoStat/NanoPlot/NanoFilt before alignment.

# Basecall with the high-accuracy model (standard choice); output is unaligned BAM
# which preserves move tables and per-read metadata as BAM tags.
dorado basecaller hac pod5_data/ > calls.bam

# Modification-aware basecalling: append the mod code to the model name.
# Probabilities land in the MM/ML BAM tags per the SAM spec.
dorado basecaller hac,5mCG_5hmCG pod5_data/ > calls_modcall.bam

# Duplex mode pairs template+complement strands from one molecule (~Q30 on duplex reads)
dorado duplex hac pod5_data/ > calls_duplex.bam

# Convert to FASTQ for tools that don't read BAM directly
samtools fastq calls.bam | gzip > calls.fastq.gz

# NanoStat: text summary (N50, mean/median Q, %>Q10/Q15/Q20)
NanoStat --fastq calls.fastq.gz --outdir nanostat_out/ --threads 4

# NanoPlot: interactive HTML (length histogram, Q distribution, length-vs-Q dot plot)
NanoPlot --fastq calls.fastq.gz --outdir nanoplot_out/ --plots dot --N50

# NanoFilt: streaming filter. For WGS assembly: -q 8 -l 3000. For variant calling: -q 10 -l 1000.
NanoFilt -q 10 -l 1000 calls.fastq.gz | gzip > calls_filtered.fastq.gz

Quality reference: Q10 = 90% per-base accuracy (R9.4.1 minimum), Q15 = 96.8% (R10.4.1 median), Q20 = 99% (R10.4.1 sup/duplex).

Goal: Reproduce NanoStat-style summary stats and diagnostic plots in Python when you need programmatic access (e.g. batching across runs) instead of the CLI report. Approach: Compute N50 by sorting lengths descending and finding where the cumulative sum crosses half the total; plot length/quality distributions to spot bimodal populations.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


def summarize_ont_reads(lengths: np.ndarray, mean_q: np.ndarray) -> dict:
    """Compute NanoStat-style summary metrics for a set of ONT reads.

    Parameters
    ----------
    lengths : per-read length in bp
    mean_q  : per-read mean Phred quality

    Returns
    -------
    dict of summary statistics, including N50 read length.
    """
    sorted_len = np.sort(lengths)[::-1]
    cumsum = np.cumsum(sorted_len)
    n50 = sorted_len[np.searchsorted(cumsum, cumsum[-1] / 2)]
    return {
        "total_reads": len(lengths),
        "total_bases": int(lengths.sum()),
        "mean_length": lengths.mean(),
        "median_length": np.median(lengths),
        "n50": int(n50),
        "mean_q": mean_q.mean(),
        "pct_q10": (mean_q >= 10).mean() * 100,
        "pct_q15": (mean_q >= 15).mean() * 100,
        "pct_q20": (mean_q >= 20).mean() * 100,
    }


# Simulate an R10.4.1 run: log-normal read lengths, ~normal quality around Q16.5
rng = np.random.default_rng(42)
n_reads = 5000
lengths = np.clip(rng.lognormal(mean=9.6, sigma=0.9, size=n_reads).astype(int), 200, 500_000)
mean_q = np.clip(rng.normal(16.5, 2.5, size=n_reads), 5, 30)

stats = summarize_ont_reads(lengths, mean_q)
for k, v in stats.items():
    print(f"{k:>15}: {v:,.2f}" if isinstance(v, float) else f"{k:>15}: {v:,}")

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].hist(lengths / 1000, bins=60, color="steelblue", edgecolor="white", linewidth=0.3)
axes[0].axvline(stats["n50"] / 1000, color="crimson", linestyle="--", label=f"N50={stats['n50']/1000:.1f} kb")
axes[0].set(xlabel="Read length (kb)", ylabel="Count", xscale="log", title="Length distribution")
axes[0].legend()
axes[1].hist(mean_q, bins=40, color="darkorange", edgecolor="white", linewidth=0.3)
axes[1].axvline(10, color="red", linestyle="--", label="Q10")
axes[1].axvline(15, color="green", linestyle="--", label="Q15")
axes[1].set(xlabel="Mean read quality", ylabel="Count", title="Quality distribution")
axes[1].legend()
plt.tight_layout()
plt.show()

Alignment with Minimap2

Goal: Map filtered long reads to a reference and assess mapping quality. Approach: Use minimap2's minimizer-based seeding with the map-ont preset (higher mismatch tolerance than short-read aligners, no splice scoring), sort/index with samtools, then read flagstat.

# map-ont: tuned for ONT genomic reads. Use map-hifi for PacBio HiFi, splice for cDNA/RNA.
minimap2 -ax map-ont -t 8 hg38.fa calls_filtered.fastq.gz \
  | samtools sort -o ont_aligned.bam -@ 8
samtools index ont_aligned.bam

# Key flagstat fields for long reads:
#   mapped (%)            -> primary mapped reads; expect 95-99% same-species genomic DNA
#   supplementary          -> chimeric/split reads, ~1-5% (SV signature; <0.1% for short reads)
#   paired-in-sequencing    -> always 0 (long reads are single-end)
samtools flagstat ont_aligned.bam

Methylation from Modified Basecalling

Goal: Get per-CpG 5mC/5hmC frequency from a modification-aware basecalling run. Approach: Basecall with a _5mCG_5hmCG model, align (MM/ML tags survive alignment), then pileup with modkit to bedMethyl.

dorado basecaller hac,5mCG_5hmCG pod5_data/ > calls_mod.bam
minimap2 -ax map-ont -t 8 --MD hg38.fa calls_mod.bam | samtools sort -o mod_aligned.bam
samtools index mod_aligned.bam

# bedMethyl columns: chrom start end name score strand thickStart thickEnd rgb coverage pct_modified
modkit pileup mod_aligned.bam methylation.bed --ref hg38.fa --cpg --combine-strands --threads 8

def load_bedmethyl(path: str, min_coverage: int = 10) -> "pd.DataFrame":
    """Load a modkit bedMethyl file and filter by minimum read coverage.

    Column 10 is coverage, column 11 is percent methylated (0-100).
    """
    import pandas as pd
    cols = ["chrom", "start", "end", "name", "score", "strand",
            "thickStart", "thickEnd", "rgb", "coverage", "pct_modified"]
    df = pd.read_csv(path, sep="\t", header=None, names=cols, usecols=range(11))
    return df[df["coverage"] >= min_coverage]

Pitfalls

• Dorado output is unaligned BAM by default — piping straight to FASTQ loses move tables and MM/ML modification tags needed for methylation calling; keep the BAM.
• Model tier tradeoff: sup is 5–10× slower than hac — don't default to sup for large runs unless accuracy (e.g. clinical SNP calling) demands it.
• FAST5 is legacy — convert to POD5 (pod5 convert fast5) before re-basecalling; Dorado does not read FAST5 directly.
• Adapters/barcodes are not trimmed by Dorado basecalling — demultiplex/trim separately (dorado demux) before downstream analysis.
• MM/ML tags are alignment-tool-dependent — confirm your aligner/sort step preserves BAM tags before running modkit; --MD is required by some pileup tools.
• Always sort before index (samtools sort then samtools index) — flagstat/pileup tools assume coordinate-sorted, indexed BAM.
• Supplementary alignment rate of 1–5% is normal for long reads (SV signature), not a QC failure — don't compare directly to short-read (<0.1%) expectations.

See Also

• long-read-sequencing — broader ONT/PacBio workflow (assembly, SV calling, isoform analysis)
• bio-applied-assembly-sv — Flye/Hifiasm assembly and Sniffles2 SV calling from these reads
• dna-methylation — bisulfite-based methylation analysis as an alternative to nanopore 5mC
• bio-alignment-files-sam-bam-basics — BAM/SAM tag structure and coordinate-sorting fundamentals