| name | long-read-sequencing |
| description | 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. |
| tool_type | bash |
| primary_tool | minimap2 |
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.
dorado basecaller hac pod5_data/ | samtools fastq > reads.fastq.gz
NanoStat --fastq reads.fastq.gz --outdir nanostat/
NanoFilt -q 10 -l 1000 reads.fastq.gz > reads_filtered.fastq.gz
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 --nano-hq reads_filtered.fastq.gz --genome-size 3g \
--out-dir flye_out/ --threads 16
hifiasm -o sample.asm -t 16 hifi_reads.fastq.gz
awk '/^S/{print ">"$2"\n"$3}' sample.asm.bp.p_ctg.gfa > assembly.fasta
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
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.
dorado basecaller hac,5mCG_5hmCG pod5/ > calls_mod.bam
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/')
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
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'] == ]
large_dels = deletions[deletions[] >= ]
()
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