Long Read Sequencing

Long-Read Sequencing

When to Use

Use this skill when:

• Processing Oxford Nanopore (ONT) or PacBio HiFi reads from POD5/FAST5/FASTQ
• Performing de novo genome or metagenome assembly (Flye, Hifiasm)
• Detecting structural variants (deletions, insertions, inversions, duplications) with Sniffles2
• Analyzing CpG methylation from ONT signal (modkit) or full-length transcript isoforms (bambu)
• Comparing ONT vs PacBio HiFi read quality, N50, or platform trade-offs

Version Compatibility

• Dorado ≥0.7 (basecalling, POD5 input), pod5-tools ≥0.3
• NanoPack (NanoStat/NanoFilt/NanoPlot) ≥1.4
• Minimap2 ≥2.28, samtools ≥1.20
• Flye ≥2.9, Hifiasm ≥0.19, Medaka ≥1.11
• Sniffles2 ≥2.4, modkit ≥0.3
• bambu (Bioconductor) ≥3.4, R ≥4.3
• Python ≥3.10 (pandas, biopython) for parsing outputs

Prerequisites

• pip install pandas biopython for output parsing shown below
• Conda/bioconda environment with dorado, minimap2, samtools, flye, hifiasm, sniffles, modkit, nanofilt, nanostat on PATH
• Familiarity with SAM/BAM (bio-alignment-files-sam-bam-basics) and VCF (bio-variant-calling-vcf-basics) formats

Quick Reference

Key Patterns

Goal: Take raw ONT POD5 signal to a sorted, indexed BAM ready for variant/SV calling. Approach: basecall with Dorado, drop low-quality/short reads with NanoFilt, align with Minimap2's map-ont preset, sort and index with samtools.

# Basecalling with Dorado (high accuracy model)
dorado basecaller hac pod5_data/ | samtools fastq > reads.fastq.gz

# QC
NanoStat --fastq reads.fastq.gz --outdir nanostat/
NanoFilt -q 10 -l 1000 reads.fastq.gz > reads_filtered.fastq.gz

# Alignment
minimap2 -ax map-ont hg38.fa reads_filtered.fastq.gz | \
    samtools sort -o aligned.bam && samtools index aligned.bam

Goal: Assemble a genome de novo from long reads and polish it to reduce residual error. Approach: use Flye for ONT reads (--nano-hq) or Hifiasm for PacBio HiFi (higher base accuracy needs no polishing pass), then run Medaka on ONT drafts to correct remaining indel/substitution errors.

# Flye for ONT genome assembly
flye --nano-hq reads_filtered.fastq.gz --genome-size 3g \
    --out-dir flye_out/ --threads 16

# Hifiasm for PacBio HiFi
hifiasm -o sample.asm -t 16 hifi_reads.fastq.gz
awk '/^S/{print ">"$2"\n"$3}' sample.asm.bp.p_ctg.gfa > assembly.fasta

# Polish ONT assembly with Medaka
medaka_consensus -i reads_filtered.fastq.gz -d flye_out/assembly.fasta \
    -o medaka/ -t 8 -m r1041_e82_400bps_hac_v4.2.0

Goal: Call structural variants (deletions, insertions, duplications, inversions, translocations) from long-read alignments. Approach: Sniffles2 clusters split/soft-clipped alignment signatures across reads; require --minsupport reads before emitting a call, then filter the VCF for PASS + support threshold.

sniffles --input aligned.bam --vcf svs.vcf \
    --reference hg38.fa --threads 8 --minsupport 5

# Keep only high-confidence calls
bcftools view -i 'FILTER="PASS" && INFO/SUPPORT>=5' svs.vcf > svs_filtered.vcf

Goal: Detect CpG methylation directly from ONT signal without bisulfite conversion. Approach: basecall with a modification-aware model (adds MM/ML tags), align (tags survive alignment), then pileup per-CpG methylation frequency with modkit.

# Basecall with methylation model
dorado basecaller hac,5mCG_5hmCG pod5/ > calls_mod.bam

# Align (MM/ML tags preserved) then pileup CpG methylation
minimap2 -ax map-ont --MD hg38.fa calls_mod.bam | samtools sort -o mod_aligned.bam
samtools index mod_aligned.bam
modkit pileup mod_aligned.bam methylation.bed --ref hg38.fa --cpg --combine-strands --threads 8

Goal: Quantify full-length transcript isoforms from long-read cDNA/direct-RNA alignments. Approach: align with Minimap2's splice preset, then run bambu in R to discover/quantify isoforms against a reference annotation.

library(bambu)
se <- bambu(reads='aligned.bam',
            annotations=gencode_gtf,
            genome=hg38_fa)
writeBambuOutput(se, path='bambu_output/')
# rowRanges(se) contains isoform coordinates
# assay(se, 'CPM') contains isoform-level expression

Technology Comparison

Pitfalls

• Basecall model selection — match model to flow cell (R9.4 vs R10.4) and kit chemistry
• Assembly genome size — always provide --genome-size to Flye for ploidy-aware assembly
• Coverage for assembly — aim for ≥50× for Flye; ≥30× for Hifiasm HiFi
• SV minimum support — default --minsupport 5 for Sniffles2; lower for low-coverage data
• Medaka model — use the correct model matching your basecaller version and flow cell

Code Templates

NanoStat QC Summary Parser

import subprocess
import re

def parse_nanostat(fastq_path):
    """Run NanoStat and return key metrics as dict."""
    out = subprocess.check_output(
        ['NanoStat', '--fastq', fastq_path, '-t', '4'],
        stderr=subprocess.DEVNULL).decode()
    metrics = {}
    for line in out.splitlines():
        m = re.match(r'(.+?):\s+([\d.,]+)', line.strip())
        if m:
            key = m.group(1).strip().lower().replace(' ', '_')
            metrics[key] = float(m.group(2).replace(',', ''))
    return metrics

stats = parse_nanostat('reads.fastq.gz')
print(f"N50: {stats.get('read_length_n50', 'N/A')} bp")
print(f"Mean Q: {stats.get('mean_read_quality', 'N/A')}")

Parse Sniffles2 VCF for SVs

import pandas as pd

def parse_sv_vcf(vcf_path):
    records = []
    with open(vcf_path) as f:
        for line in f:
            if line.startswith('#'):
                continue
            fields = line.strip().split('\t')
            chrom, pos, sv_id, ref, alt = fields[:5]
            info = dict(
                kv.split('=', 1) if '=' in kv else (kv, True)
                for kv in fields[7].split(';')
            )
            records.append({
                'chrom': chrom, 'pos': int(pos),
                'svtype': info.get('SVTYPE', ''),
                'svlen': abs(int(info.get('SVLEN', 0))),
                'support': int(info.get('SUPPORT', 0)),
                'af': float(info.get('AF', 0)),
            })
    return pd.DataFrame(records)

svs = parse_sv_vcf('sniffles_svs.vcf')
deletions = svs[svs['svtype'] == 'DEL']
large_dels = deletions[deletions['svlen'] >= 1000]
print(f"Large deletions (≥1 kb): {len(large_dels)}")

Assembly N50 Calculator

from Bio import SeqIO

def assembly_stats(fasta_path):
    lengths = sorted([len(r.seq) for r in SeqIO.parse(fasta_path, 'fasta')],
                     reverse=True)
    total = sum(lengths)
    cumsum = 0
    n50 = 0
    for l in lengths:
        cumsum += l
        if cumsum >= total * 0.5:
            n50 = l
            break
    return {
        'num_contigs': len(lengths),
        'total_length': total,
        'largest': lengths[0],
        'n50': n50,
    }

stats = assembly_stats('assembly.fasta')
print(f"N50 = {stats['n50']:,} bp | Total = {stats['total_length']:,} bp")

See Also

• bio-applied-genome-assembly — assembly algorithms, contiguity/quality metrics, annotation
• bio-applied-assembly-sv — structural variant calling and filtering in depth
• dna-methylation — bisulfite/short-read methylation calling and DMR analysis
• bio-applied-isoform-analysis — full-length transcript isoform quantification