Bio Long Read Sequencing Structural Variants

Version Compatibility

Reference examples tested with: Sniffles 2.2+, cuteSV 2.1+, minimap2 2.28+, 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

Results depend on inputs that outlive the binary version - record them:

• The reference-matched tandem-repeat BED supplied to the caller (Sniffles --tandem-repeats) drives the FP rate in repeats more than any other setting. Record which TR BED was used.
• Benchmark numbers depend on the region set + TR handling + Truvari params; record all three.
• cuteSV parameters are platform-specific (ONT vs HiFi vs CLR); the defaults are not platform-appropriate.

If code throws an error, introspect the installed tool (sniffles --help, cuteSV --help) and adapt the example to the actual API rather than retrying.

Long-Read Structural Variants

"Find structural variants in my long reads" -> Map with the SV-ready preset (soft-clipped supplementaries), call with a TR-aware caller, and benchmark stating the region set and Truvari params.

• CLI: sniffles --input aln.bam --vcf svs.vcf --reference ref.fa --tandem-repeats TR.bed

Long reads are the killer app for SVs: a single read spans the breakpoint (within-read CIGAR or split alignment) and resolves repeats short reads cannot. By convention SV = >=50 bp; the 30-100 bp range is a VNTR-dominated gray zone where callers disagree most.

The Single Most Important Modern Insight -- An SV Call Is a Representation Artifact as Much as a Biological Fact

In tandem repeats and segmental duplications, the same biological event has many valid VCF encodings - a deletion can be written as the reciprocal insertion on the other allele, and a VNTR expansion's breakpoints slide freely across repeat units. Consequently:

1. The tandem-repeat BED, the aligner, and the Truvari parameters decide precision/recall as much as the caller does. A claim like "caller X has F1 0.95" is meaningless without also stating the region set, the TR BED supplied to the caller, and the Truvari params - change any one and the number moves more than the gap between callers.
2. Without a TR BED, one event fragments into several false-positive calls with inconsistent breakpoints. --tandem-repeats makes clustering repeat-aware (widening the merge window inside annotated TRs) - the single biggest FP-reduction lever, not a nicety.
3. truvari refine exists precisely to re-harmonize representations within TR regions; benchmarking TR-dense regions without it systematically understates recall.

Caller Taxonomy

Decision Tree by Scenario

Alignment for SV Calling

Map with minimap2 (the modern default; NGMLR is a higher-precision/slower legacy niche for Sniffles). Use the platform preset and keep soft-clipped supplementary alignments - split-read callers reconstruct breakpoints from the clipped sequence on those records.

minimap2 -ax map-ont --MD -Y ref.fa ont.fq.gz | samtools sort -o aln.bam && samtools index aln.bam
#   -Y keeps SEQ on supplementaries (the SV substrate); --MD for cuteSV; map-hifi/map-pb for PacBio

Sniffles2 - germline and the .snf population workflow

# Single sample (always supply --reference for INS sequence and --tandem-repeats for repeats)
sniffles --input aln.bam --vcf svs.vcf --reference ref.fa --tandem-repeats human_GRCh38_TR.bed

# Cohort: per-sample .snf signature index, then merge + joint-genotype
sniffles --input s1.bam --snf s1.snf --reference ref.fa --tandem-repeats TR.bed
sniffles --input s2.bam --snf s2.snf --reference ref.fa --tandem-repeats TR.bed
sniffles --input s1.snf s2.snf --vcf cohort.vcf --reference ref.fa

# Force-call / regenotype a known SV set in a new sample
sniffles --input new.bam --genotype-vcf known_svs.vcf --vcf genotyped.vcf

# Single-sample low-VAF / mosaic (NOT a tumor-normal caller)
sniffles --input tumor.bam --vcf mosaic.vcf --mosaic

The .snf is a binary signature index (NOT a VCF - never bcftools it); it retains sub-threshold signatures so the merge re-genotypes an SV even in a sample that did not independently pass support.

cuteSV - the per-platform parameter trap

cuteSV's defaults are not platform-appropriate; the README gives distinct sets by error rate. --genotype is OFF by default. Positional args: cuteSV <bam> <ref> <out.vcf> <work_dir>. Force-calling moved to the separate cuteFC tool.

mkdir cutesv_work
cuteSV aln.bam ref.fa cutesv.vcf cutesv_work --genotype \
  --max_cluster_bias_INS 100 --diff_ratio_merging_INS 0.3 \
  --max_cluster_bias_DEL 100 --diff_ratio_merging_DEL 0.3   # ONT set

Benchmarking with Truvari

truvari bench --base giab_tier1.vcf.gz --comp calls.vcf.gz \
  --includebed tier1_regions.bed --pctseq 0.7 --refdist 500 --passonly -o bench/
truvari refine bench/        # re-harmonize TR-region representations for a fair comparison

--pctseq (default 0.7) compares the actual inserted/deleted sequence, not just coordinates - set 0 for depth-based callers lacking alt sequence, keep 0.7 for long-read callers. Region set dominates the headline: Tier1 (resolvable INS/DEL >=50 bp) overstates whole-genome performance; CMRG reflects hard clinical loci. Tier1 v0.6 is INS/DEL only - do not report INV recall against it.

Per-Method Failure Modes

One VNTR fragments into many false positives

Trigger: calling in tandem repeats without a TR BED. Mechanism: the breakpoint slides across repeat units, scattering signatures. Symptom: several calls with inconsistent breakpoints where one event exists. Fix: supply --tandem-repeats to the caller; truvari refine when benchmarking.