Bio Workflows Longread Sv Pipeline

Version Compatibility

Reference examples tested with: minimap2 2.28+, Sniffles 2.2+, cuteSV 2.1+, pbsv 2.9+, dipcall 0.3+, bcftools 1.19+, samtools 1.19+, truvari 4.0+

Before using code patterns, verify installed versions match. If versions differ:

• CLI: <tool> --version then <tool> --help to confirm flags

Use minimap2 >= 2.28: the lr:hq accurate-read preset was added in 2.27, and 2.28 fixes the 2.27 --MD regression. Supply a reference-matched tandem-repeat BED to the caller - it is the single biggest false-positive lever in repeats. Truvari renamed the alt-sequence-similarity param from --pctsim to --pctseq at v4; confirm against truvari bench --help.

If code throws an error, introspect the installed tool and adapt the example to the actual API rather than retrying.

Long-Read SV Pipeline

"Detect structural variants from my long-read sequencing data" -> Chain basecalling, platform-matched minimap2 alignment, an SV caller selected by platform and goal, an optional assembly-based branch, cohort merging, and a parameterized Truvari benchmark - with the SV mechanism delegated to the component skills.

This is a workflow (orchestration) skill: it makes the stage-to-stage decisions and quality gates that connect the component skills. It does NOT re-teach how a caller sees an SV or the VCF representation minefield - that lives in variant-calling/structural-variant-calling and long-read-sequencing/structural-variants.

Why long reads for SV: the physics that justifies this pipeline

Run this pipeline instead of a short-read SV workflow for one mechanistic reason: a single long read (PacBio HiFi ~15-25 kb, ONT tens of kb to >Mb ultralong) physically spans the SV and both flanks in one molecule, turning SV detection from an inference-over-fragments problem into near-direct observation. Two consequences decide whether long reads earn their per-sample cost:

• Insertions become tractable. Placing and sizing an insertion needs reads that carry the novel bases; when an INS exceeds Illumina read length (150 bp) no short read spans it, so short-read INS recall is stuck at ~30-50% while long reads reach ~90%+. Ebert 2021 (Science 372:eabf7117) found 68% of 107,590 assembly-discovered SVs were missed by short reads. If insertions matter, this pipeline is the answer, not a tuning knob.
• Repeat-mediated junctions resolve. A 20 kb read anchored in unique sequence on both flanks spans a breakpoint buried in a 5 kb repeat that no 150 bp read can straddle, bringing segmental-duplication NAHR, mobile-element insertions, and VNTR/STR expansions into reach.

The governing pipeline principle: an SV call is a representation artifact of choices made upstream. The aligner preset, the tandem-repeat BED handed to the caller, and the Truvari matching parameters move precision/recall as much as the caller does. Chaining decisions - not the caller name - are what this skill is about.

Pipeline map

POD5/FAST5 (ONT only)
    |  [Step 0] Dorado basecall  -> long-read-sequencing/basecalling
    v            (model + methylation are IRREVERSIBLE choices; sup model for SV)
FASTQ (ONT / PacBio HiFi)
    |  [QC]     NanoPlot / NanoComp -> long-read-sequencing/long-read-qc
    v            gate: read N50 >10 kb, sane quality, chimera screen
[Step 1] minimap2 alignment       -> long-read-sequencing/long-read-alignment
    |            preset by platform/chemistry; -Y keeps breakpoint seq on split reads
    v            gate: mapping rate >90%, coverage >=15x
[Step 2] SV calling               -> long-read-sequencing/structural-variants
    |            Sniffles2 / cuteSV / pbsv, caller by platform + goal
    |            + tandem-repeat BED (biggest FP lever)      variant-calling/structural-variant-calling
    v
[Step 3, optional] assembly-based SV (dipcall / PAV against a phased diploid assembly)
    |            highest-quality callset; the way truth sets are built
    v
[Step 4] cohort merge             -> two-step Sniffles2 .snf (per-sample -> combine)
    |
    v
[Step 5] benchmark                -> Truvari vs GIAB HG002 Tier 1 + CMRG
                 an F1 is meaningless without refdist/pctsize/pctseq

Platform decision: ONT vs PacBio HiFi

The platform sets the preset, the caller options, and what bonus channels are available. Decide before basecalling.

R10.4 chemistry moved ONT simplex to ~Q20, which is why lr:hq (not the noisy-read map-ont) is the right preset for modern ONT - it rewrites the scoring/chaining model for accurate reads and runs faster at equal accuracy. Older R9 data still needs map-ont.

SV caller selection (by platform and goal)

Do not re-derive the caller mechanism here; pick by goal and hand tuning to the component skill.

Methods evolve; verify current best practice against each tool's docs before committing. Deeper caller tuning (cuteSV per-platform params, the tandem-repeat BED, aligner effects) lives in long-read-sequencing/structural-variants.

Step 0: Basecalling (ONT only)

Goal: Convert raw POD5/FAST5 signal into reads suitable for SV calling, capturing methylation if it will ever be needed.

Approach: Basecall with Dorado using a chemistry-matched sup (super-accuracy) model; request modified bases at basecall time because methylation cannot be recovered later. PacBio HiFi arrives as reads already, so this step is skipped.

# sup model maximizes accuracy for SV; 5mCG_5hmCG requested now (irreversible if omitted).
# See long-read-sequencing/basecalling for model selection and duplex.
dorado basecaller sup pod5_dir/ --modified-bases 5mCG_5hmCG > reads.bam
samtools fastq -T MM,ML reads.bam | gzip > reads.fastq.gz   # -T carries methylation tags through

Step 1: Alignment

Goal: Produce a sorted, indexed BAM whose split (supplementary) alignments retain the breakpoint sequence SV callers reconstruct from.

Approach: Align with the platform-matched minimap2 preset; keep -Y so supplementary alignments are soft-clipped (not hard-clipped), which preserves the junction bases on split reads. SV calling rides on supplementary, not secondary, alignments.

# ONT R10/Q20: lr:hq (accurate reads). Older R9: map-ont. HiFi: map-hifi. PacBio CLR: map-pb.
minimap2 -ax lr:hq -t 16 --MD -Y reference.fa reads.fastq.gz | \
    samtools sort -@ 4 -o aligned.bam
samtools index aligned.bam

QC checkpoint (gate before spending compute on calling):

samtools flagstat aligned.bam                                  # mapping rate should be >90%
samtools depth -a aligned.bam | awk '{s+=$3} END{print "mean cov:", s/NR}'
# Gate: >=15x for confident SV calling; below ~10x callers drift toward false negatives.

Step 2: SV calling

Goal: Call SVs (>=50 bp DEL/INS/DUP/INV/BND) from the aligned reads with a caller matched to the platform.

Approach: Run Sniffles2 (the default) with a reference-matched tandem-repeat BED - it clusters the repeat-driven false positives that otherwise dominate the callset. --minsvlen 50 enforces the GIAB >=50 bp SV convention (Sniffles2 defaults to 35).

# Sniffles2: the tandem-repeat BED is the single biggest false-positive lever in repeats.
sniffles --input aligned.bam --reference reference.fa \
    --tandem-repeats human_GRCh38_TR.bed \
    --vcf svs.vcf.gz --threads 8 --minsvlen 50 --output-rnames

cuteSV as an alternative - its defaults are NOT platform-appropriate, and --genotype is off by default:

# ONT param set shown. HiFi: 1000/0.9/1000/0.5. CLR: 100/0.3/200/0.5. See structural-variants.
mkdir -p work_dir   # cuteSV requires the work dir to pre-exist; it does not create it
cuteSV aligned.bam reference.fa svs.vcf work_dir/ --threads 8 --genotype \
    --max_cluster_bias_INS 100 --diff_ratio_merging_INS 0.3 \
    --max_cluster_bias_DEL 100 --diff_ratio_merging_DEL 0.3

Step 3 (optional): Assembly-based SV

Goal: Produce the highest-quality SV callset by comparing a phased diploid assembly to the reference, rather than inferring from read alignments.

Approach: Assemble the genome (hifiasm/verkko), then call variants from the two haplotype assemblies aligned to the reference. dipcall (the syndip method) and PAV (the HGSVC method) are the assembly-vs-reference callers; this is how the GIAB and HGSVC truth sets themselves are built. Use it when an assembly already exists or when callset quality outranks turnaround.

# dipcall needs two haplotype assemblies (hap1/hap2) plus minimap2/k8/htsbox on PATH.
run-dipcall prefix reference.fa hap1.fa hap2.fa > prefix.mak
make -j2 -f prefix.mak                                          # emits prefix.dip.vcf.gz + prefix.dip.bed

Step 4: Cohort merging (the two-step .snf design)

Goal: Build a joint-genotyped multi-sample SV matrix, not a union of per-sample discovery VCFs.

Approach: Sniffles2's population design processes each sample independently into a compact .snf, then combines the .snf files in a second pass - scaling linearly in N. A union of per-sample discovery VCFs is wrong: a sample recorded 0/0 may simply not have had that event discovered in it (a false missing), which corrupts allele frequencies. The two-step .snf combine force-genotypes every sample at every merged site.

# Pass 1: per-sample .snf (each sample processed once, independently).
for s in sample1 sample2 sample3; do
    sniffles --input ${s}.bam --reference reference.fa \
        --tandem-repeats human_GRCh38_TR.bed --snf ${s}.snf
done
# Pass 2: combine into a jointly genotyped cohort VCF (linear in N).
sniffles --input sample1.snf sample2.snf sample3.snf --vcf cohort.vcf.gz

For sequence-aware AF work across callsets, prefer Truvari collapse over position-only merging (position-only mergers inflate allele frequency by up to 2.2x; English 2022) - see variant-calling/structural-variant-calling for the merger decision table.

Step 5: Benchmarking - an SV F1 is meaningless without its parameters

Goal: Report a defensible, reproducible accuracy figure - not a number that looks good because of loose matching.

Approach: Truvari bench counts a call as a true positive only if it matches a truth variant under ALL of --refdist, --pctsize, and --pctseq simultaneously. Every one of these moves the score, so a bare F1 is uninterpretable. Report the full parameter set, run truvari refine for a harmonized re-comparison, and stratify by region.

# Report EVERY parameter. --pctseq 0 disables alt-sequence checking and quietly inflates INS scores.
truvari bench -b HG002_SV_Tier1.vcf.gz -c svs.vcf.gz -o bench_tier1/ --passonly \
    -f reference.fa --refdist 500 --pctsize 0.70 --pctseq 0.70 --sizemin 50   # -f persists reference to params.json
truvari refine bench_tier1/                                     # harmonized breakpoint re-comparison (needs the reference bench recorded)

# CMRG is NOT optional: Tier 1 EXCLUDES the medically relevant repetitive genes.
truvari bench -b HG002_CMRG_SV.vcf.gz -c svs.vcf.gz -o bench_cmrg/ --passonly \
    --refdist 500 --pctsize 0.70 --pctseq 0.70 --sizemin 50

Three escalating bars are routinely conflated, and a pipeline can pass the first while failing the ones that matter:

• Event detection - something of about the right type/size near the right place (loose refdist, no sequence check). Easy.
• Breakpoint accuracy - POS/END within a few bp (tight --refdist, --pctseq on). Matters at exon/splice boundaries.
• Genotype accuracy - the sample GT (het/hom) is correct (genotype-aware comparison). A caller can detect an event perfectly and still call het-as-hom, which is fatal for Mendelian analyses.

Region stratification is decisive: Tier 1 (Zook 2020 Nat Biotechnol 38:1347) is conservative isolated SVs, while CMRG (Wagner 2022 Nat Biotechnol 40:672) covers the repetitive medically relevant genes Tier 1 leaves out - where GRCh38 false duplications cause reference-specific misses that masking raised from 8% to 100% recall. A good Tier 1 F1 certifies nothing about the genes clinicians care about; run both.