# Long Read Sequencing

> Basecall ONT POD5 with Dorado, align with Minimap2, assemble with Flye/Hifiasm, call SVs with Sniffles2. Use when basecalling nanopore reads, doing long-read assembly, SV calling, or ONT methylation/isoform analysis.

- Skill: `pavel-kravchenko/long-read-sequencing` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/long-read-sequencing`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/long-read-sequencing/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

| Task | Tool | Key Command |
|------|------|-------------|
| ONT basecalling | Dorado | `dorado basecaller hac pod5/` |
| Read QC | NanoStat | `NanoStat --fastq reads.fq.gz` |
| Read filtering | NanoFilt | `NanoFilt -q 10 -l 1000` |
| Alignment | Minimap2 | `minimap2 -ax map-ont ref.fa reads.fq.gz` |
| HiFi alignment | Minimap2 | `minimap2 -ax map-hifi ref.fa reads.fq.gz` |
| cDNA alignment | Minimap2 | `minimap2 -ax splice reads.fq.gz` |
| De novo assembly (ONT) | Flye | `flye --nano-hq reads.fq.gz --genome-size 5m` |
| De novo assembly (HiFi) | Hifiasm | `hifiasm -o out.asm -t 16 reads.fq.gz` |
| Assembly polishing | Medaka | `medaka_consensus -i reads.fq.gz -d assembly.fa` |
| Assembly QC | QUAST | `quast.py assembly.fa -r ref.fa` |
| SV calling | Sniffles2 | `sniffles --input aligned.bam --vcf svs.vcf` |
| Methylation | modkit | `modkit pileup aligned.bam meth.bed --cpg` |
| Isoform analysis | bambu (R) | `bambu(reads='aln.bam', annotations=gtf)` |

## 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.
```bash
# 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.
```bash
# 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.
```bash
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.
```bash
# 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.
```r
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

| Property | ONT R9/R10 | PacBio HiFi |
|----------|-----------|-------------|
| Read length | Typically 5–50 kb | Typically 15–25 kb |
| Raw accuracy | 97–99% (R10) | >99.9% |
| Throughput | High (P2 Solo: 80 Gb) | Moderate (Sequel II: 160 Gb) |
| Methylation | Direct (native DNA) | 5mC with Kinetics |
| Cost per Gb | Low | Higher |

## 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
```python
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
```python
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
```python
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