Bio Variant Calling Structural Variant Calling

Version Compatibility

Reference examples tested with: Manta 1.6+, DELLY 1.2+, GRIDSS 2.13+, smoove 0.2+, SvABA 1.1+, bcftools 1.19+, samtools 1.19+, SURVIVOR 1.0.7+, Truvari 4.0+, Sniffles2 2.2+, cuteSV 2.0+

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

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

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Note: DELLY 1.x and Truvari 4.x renamed flags/defaults from earlier lines (Truvari's alt-sequence-similarity param was --pctsim, now --pctseq); always confirm against --help before trusting old muscle memory.

Structural Variant Calling

"Call structural variants from my WGS data" -> Reconstruct large rearrangements (DEL, INS, INV, DUP, BND) that no single read reports directly, by fusing evidence from how alignments break, disagree, or deviate from expectation.

• CLI (short read): configManta.py (Manta), delly call, smoove call (LUMPY), gridss, svaba run
• CLI (long read): sniffles, cuteSV, pbsv, dipcall.kit/pav (assembly-based)

The governing principle: SV discovery is signal reconstruction, not base-calling

An SV caller is not a pileup genotyper. SVs are never observed directly by a read; they are inferred from four orthogonal signals, and each caller fuses only a subset. Its blind spots follow directly from which signals it omits (Mahmoud 2019 Genome Biol 20:246 is the canonical signal-to-blind-spot map).

1. Discordant read pairs (RP) - a pair whose insert size or orientation departs from the library distribution. A DEL inflates the apparent insert; a tandem DUP gives everted (+/+/-/-) pairs; an INV flips one strand; an interchromosomal event splits mates across chromosomes. RP never sees the junction base, so it localizes a breakpoint only to fragment-size uncertainty (~+/-300-500 bp).
2. Split reads / soft-clips (SR) - one read aligning in two pieces to non-contiguous positions, recorded via the BAM SA (supplementary alignment) tag. The clip point IS the breakpoint at base resolution (~+/-0-10 bp). Requires a read that straddles the junction with unique anchor on both sides - impossible inside a repeat longer than the read or for an INS larger than (read length - anchor).
3. Read depth (RD) - copy-number shifts local coverage (het DEL ~0.5x, DUP ~1.5x, hom DEL ~0). The ONLY signal for large CNVs whose breakpoints fall in unmappable repeats; blind to balanced events (INV, balanced BND) that leave dosage unchanged. Resolution is bin-size-limited (100 bp-1 kb), the worst of any signal.
4. Local assembly (AS) - reconstruct the breakpoint-spanning contig de novo from clipped/discordant reads, then realign it. Recovers the exact junction sequence, microhomology, and inserted bases; resolves events no single read spans. Most powerful and most expensive; separates GRIDSS/Manta/SvABA from LUMPY.

The real tradeoff is sensitivity vs breakpoint resolution (not vs specificity)

• RP-heavy calling maximizes SENSITIVITY (a pair spanning a junction is easy, works at low depth) but yields IMPRECISE breakpoints (+/-hundreds bp, CIPOS=-289,301).
• SR / assembly calling maximizes RESOLUTION (base-precise, CIPOS=0,0) but LOSES sensitivity wherever a read cannot cleanly straddle the junction (repeats, low depth, large INS).

Below ~30x the SR signal thins (fewer reads straddle any junction) and callers silently drift into the low-resolution RP-only regime. A caller demanding SR confirmation reports beautiful breakpoints and misses the repeat-mediated events that matter clinically; a caller accepting RP-only calls finds more but hands back wide confidence intervals.

Caller selection

Consensus recipe: run DELLY + Manta + SvABA + GRIDSS2 (the latter GRIPSS-filtered first, since its raw VCF is a breakpoint graph, not a callset) and require >=2/4 agreement - singletons are enriched for false positives, so intersection trades a little recall for a large precision gain. Caveat: this caps INS recall (LUMPY-class blind spots drag the intersection down); for INS-heavy work prefer one assembly caller + a graph genotyper, or long reads. Methods evolve - verify current ensemble practice against tool docs before committing.

SV detection limits by platform

Insertions are a physics limit, not a tuning failure. Placing and sizing an INS requires reads spanning both junctions plus the novel bases; when the INS exceeds read length (Illumina 150 bp), no read spans it and the inserted sequence is unrecoverable without assembly - and short-read assembly fails again when the insertion is a repeat (mobile element, VNTR). GIAB HG002 Tier 1 contains MORE INS than DEL (7,281 INS vs 5,464 DEL >=50 bp; Zook 2020 Nat Biotechnol 38:1347), yet INS are the hardest short-read class. Ebert 2021 (Science 372:eabf7117) found 68% of 107,590 assembly-discovered SVs were missed by short reads. Do NOT benchmark short-read INS against a long-read truth set and blame the caller - ~30-50% is the ceiling by construction.

The VCF representation minefield

The VCF spec represents SVs three ways and callers disagree on all three - this is where careful people lose days.

• SVTYPE - the class (DEL/INS/DUP/INV/BND).
• END (INFO) - the other breakpoint on the same chromosome. POS is the last unaffected anchor base, so a DEL at POS=1000, END=2000 removes 1001-2000. Off-by-one here corrupts every length calc. (VCF 4.4 began deprecating INFO/END for symbolic alleles; bcftools/GATK still emit and consume the classic field.)
• SVLEN - SIGNED length: DEL is NEGATIVE, INS/DUP positive by historical convention (Manta, DELLY follow it; some tools emit absolute values). NEVER filter on raw SVLEN without ABS() or a naive SVLEN >= 50 drops every deletion.
• CIPOS / CIEND - confidence intervals encoding the RP-vs-SR resolution story: an SR-resolved breakpoint is CIPOS=0,0, an RP-only one is CIPOS=-289,301. This is the ground truth for "how much do I trust this breakpoint" and exactly what a population merger must respect.
• IMPRECISE (flag) - set when the breakpoint is RP/depth-derived, not junction-resolved; its absence implies PRECISE. An all-IMPRECISE callset cannot be merged tightly.

Breakend (BND) notation - the part everyone gets wrong. A rearrangement whose two sides are not a simple same-chromosome span is written as PAIRED BND records linked by MATEID, with an ALT string whose bracket direction encodes the join orientation:

2   321681  bnd_V  T  ]13:123456]T   MATEID=bnd_U
13  123456  bnd_U  A  A[2:321681[    MATEID=bnd_V

One reciprocal translocation is 2 BND records; chromothripsis is a graph of dozens. The SAME biological inversion can appear as <INV> in Manta, as 2+ BND records in GRIDSS, and as a DEL+DUP artifact pair in a naive RD caller - a merger that does not understand this triple-counts or drops the event.

Manta

Manta builds a genome-wide breakend association graph from RP+SR, then does targeted local assembly of each candidate breakend and realigns the contig for base-resolution breakpoints. candidateSmallIndels.vcf.gz is the recommended indel-candidate input for Strelka2.

configManta.py --bam sample.bam --referenceFasta reference.fa --runDir manta_run
manta_run/runWorkflow.py -j 8
# results/variants/: diploidSV.vcf.gz (germline), candidateSV.vcf.gz (unscored superset),
#   candidateSmallIndels.vcf.gz (Strelka2 input)

WES/panel and RNA need explicit modes - default depth filtering assumes WGS-uniform coverage and silently drops true SVs at high-depth targeted loci:

# --exome disables the high-depth filter; --rna handles split alignments across splice junctions
configManta.py --bam sample.bam --referenceFasta reference.fa --exome \
    --callRegions regions.bed.gz --runDir manta_exome

Tumor-normal somatic mode adds somaticSV.vcf.gz:

configManta.py --tumorBam tumor.bam --normalBam normal.bam \
    --referenceFasta reference.fa --runDir manta_somatic
manta_somatic/runWorkflow.py -j 8

DELLY

BCF output by default; the scalable cohort pattern is site-list-then-regenotype (Section below), not one big multi-BAM call.

delly call -g reference.fa -o sv_calls.bcf sample.bam
bcftools view sv_calls.bcf > sv_calls.vcf

# Per-type calling (INS recovers only small, SR-anchored insertions)
delly call -t DEL -g ref.fa -o deletions.bcf sample.bam   # also DUP, INV, BND, INS

# Somatic: call tumor+normal, then classify with a sample sheet
delly call -g reference.fa -o svs.bcf tumor.bam normal.bam
printf 'tumor\ttumor\nnormal\tcontrol\n' > samples.tsv
delly filter -f somatic -o somatic_svs.bcf -s samples.tsv svs.bcf

smoove (LUMPY, with genotyping + depth annotation)

Nobody runs raw lumpyexpress anymore. smoove wraps LUMPY + svtyper + duphold with sane defaults and a cohort workflow. duphold annotates each DEL/DUP with three depth fold-changes: over the whole event (DHFC), over GC/mappability-matched bins (DHBFC), and over the immediate flanks (DHFFC) - a cheap, high-value filter for RP-only calls lacking depth support. duphold's own recommended false-positive filters are DHFFC < 0.7 for deletions (flanking fold-change, its most robust metric) and DHBFC > 1.3 for duplications.

smoove call --name sample --fasta reference.fa --outdir smoove_out -p 8 sample.bam
# smoove_out/sample-smoove.genotyped.vcf.gz

# Filter RP-only false positives that lack depth corroboration (duphold's recommended fields)
bcftools view -i '(SVTYPE!="DEL" && SVTYPE!="DUP") || (SVTYPE="DEL" && FMT/DHFFC<0.7) || (SVTYPE="DUP" && FMT/DHBFC>1.3)' \
    smoove_out/sample-smoove.genotyped.vcf.gz > smoove.dhfc.vcf

GRIDSS: the raw VCF is a breakpoint graph, not a callset

GRIDSS2 is the only short-read caller doing genome-wide break-end assembly (positional de Bruijn graph) before calling, and the only one reporting single break-ends (SGL) where just one side maps - viral integrations, centromeric/telomeric junctions. But its raw VCF is deliberately noisy low-level breakpoints; it is NOT usable directly. The intended chain is GRIDSS -> GRIPSS (filtering/linkage) -> PURPLE (purity/ploidy/copy-number) -> LINX (interprets breakpoints into chromothripsis, breakage-fusion-bridge, fusions).

# --assembly must be a writable path; reference needs .fai and .dict
gridss --reference reference.fa --output gridss_raw.vcf \
    --assembly gridss_assembly.bam --threads 8 sample.bam

# Somatic: GRIDSS2 tumor+normal, then GRIPSS filtering
gridss --reference reference.fa --output gridss_raw.vcf --assembly asm.bam \
    --labels normal,tumor --threads 8 normal.bam tumor.bam
# GRIPSS (launched as `java -jar gripss.jar` in practice) also requires PON inputs
# (-pon_sgl_file, -pon_sv_file, -known_hotspot_file); confirm the full flag set with its docs.
gripss -sample tumor -reference normal -ref_genome reference.fa \
    -ref_genome_version 38 -pon_sgl_file sgl.pon -pon_sv_file sv.pon \
    -vcf gridss_raw.vcf -output_dir gripss_out/

Genotyping is not discovery (do NOT union discovery VCFs into a population matrix)

Discovery answers "is there an SV here, and what/where?"; genotyping answers "what is each sample's GT (0/0, 0/1, 1/1) here?". These are different statistical problems. A sample recorded 0/0 in a discovery-VCF union may simply not have had that event DISCOVERED in it - a false "missing", not a true reference genotype. Building an allele-frequency-quality cohort matrix requires FORCE-GENOTYPING every sample at every discovered site.

Correct cohort workflow: discover per-sample -> merge to a non-redundant SITE list -> re-genotype ALL samples at ALL sites -> merge genotypes.

# DELLY's native squared-off pattern: per-sample call -> merge SITES -> regenotype at sites
delly call -g ref.fa -o s1.bcf s1.bam        # ... one per sample
delly merge -o sites.bcf s1.bcf s2.bcf s3.bcf
delly call -g ref.fa -v sites.bcf -o s1.geno.bcf s1.bam   # regenotype each at the union sites
bcftools merge -m id -Ob -o cohort.bcf s1.geno.bcf s2.geno.bcf s3.geno.bcf

Population merging: the merge parameters ARE the result

"The same SV" across samples/callers is a fuzzy concept because breakpoints disagree by CIPOS margins. Position-only merging inflates allele frequency by up to 2.2x versus sequence-aware merging, because it collapses distinct nearby alleles into one (English 2022 Genome Biol 23:271). For any AF-dependent analysis (constraint, association, catalogs), always report the merger and its parameters.

# SURVIVOR: max_dist=1000 min_callers=2 type_agree=1 strand_agree=1 est_dist=0 min_size=50
ls manta.vcf delly.vcf gridss.vcf smoove.vcf > vcf_list.txt
SURVIVOR merge vcf_list.txt 1000 2 1 1 0 50 merged.vcf   # single-sample caller consensus only