How to Detect Structural Variants in Your Own Whole Genome

By the end of this guide you will have a VCF containing several thousand structural variant calls for one genome. Each call is annotated with its calling evidence, its allele frequency in population databases, and the genes it intersects. You need an aligned BAM or CRAM from 30x or deeper short-read whole-genome sequencing against GRCh38. The alignment must come from a caller-aware aligner, meaning bwa-mem2 or dragen-os rather than a soft-clip-averse aligner. You also need roughly 100 GB of scratch disk, 32 GB of RAM, and 8 or more cores. Wall time on a 16-core machine is four to six hours for the calling steps. Everything here runs from the command line. The tools are the following.

• samtools
• manta
• delly
• smoove (a wrapper around lumpy and svtyper)
• cnvpytor
• truvari
• AnnotSV

If your data came from long reads, skip to step 7, where the approach changes substantially and gets better.

Two facts frame the whole exercise. First, structural variants (SVs, conventionally any rearrangement of 50 bp or more: deletions, duplications, insertions, inversions, and translocations) account for more differing bases between two human genomes than single-nucleotide variants do, even though they are far less numerous. A typical genome carries on the order of 20,000 SV calls against the reference, and the large ones dominate the base-pair budget. Second, short-read SV calling is genuinely hard, and every caller you run will disagree with every other caller. The point of a multi-caller pipeline is not consensus for its own sake but a calibrated sense of which calls you should believe.

1. Check that your alignment can support SV calling

Short-read SV detection draws on four signal types.¹

• discordant read pairs whose insert size or orientation contradicts the reference
• split reads that align partly to one locus and partly to another
• read depth that rises or falls across an interval
• local assembly of unmapped or soft-clipped reads

Which signals are available to you is decided at alignment time, so verify this before calling anything.

samtools view -H sample.cram | grep '^@PG'
samtools stats -@8 sample.cram > sample.stats
grep '^SN' sample.stats | egrep 'insert size|average length|reads mapped|reads duplicated'

You want mean insert size in the 300–500 bp range with a standard deviation under about 15 percent of the mean. Average read length should be 150 bp, with mapped fraction above 0.98 and duplicate rate below 0.15. A wide insert-size distribution destroys discordant-pair sensitivity because callers model the distribution and flag outliers. If the standard deviation is 40 percent of the mean, the tails swallow real signal. Also confirm that supplementary alignments are present, since split-read evidence lives in them:

samtools view -c -f 2048 -@8 sample.cram   # should be millions, not zero

If that count is zero, your BAM was produced by a pipeline that discarded supplementary alignments, and you should realign from FASTQ. Keep soft clips. Do not apply any filter that removes reads with mapping quality below 20 at this stage, because callers use low-MAPQ reads in repeat-adjacent breakpoints deliberately.

2. Run Manta first

Manta is our default and the single caller we would keep if forced to pick one. It builds breakend graphs from discordant pairs, then does local de novo assembly at each candidate locus. That gives it base-pair-resolution breakpoints and, unusually for a short-read caller, real insertion sequences up to a few hundred base pairs. In the most systematic benchmark of short-read SV callers to date, Manta was among the top performers across deletions, duplications, and inversions, with the important caveat that no single algorithm dominated across all SV types and size ranges.²

configManta.py \
  --bam sample.cram \
  --referenceFasta GRCh38.fa \
  --callRegions calling_regions.bed.gz \
  --runDir manta_run

manta_run/runWorkflow.py -m local -j 16 -g 32

The --callRegions argument matters more than people expect. Pass a BED of the primary chromosomes with centromeres and the major segmental duplication blocks excluded, or Manta will spend hours assembling garbage in alpha satellite and hand you thousands of false duplications. Use the Manta-recommended approach of including chr1–22, X, Y, and M, and subtracting the ENCODE blacklist plus the GRCh38 centromere/heterochromatin track from UCSC.

Your outputs are in manta_run/results/variants/. Use diploidSV.vcf.gz. Do not use candidateSV.vcf.gz. The candidate file has no quality filtering and is meant as caller input, so it is not an answer. Expect 6,000–10,000 PASS records for one 30x genome, heavily skewed toward deletions under 1 kb.

3. Add two callers with different failure modes

Callers fail in correlated ways when they use the same evidence, so pick complements rather than more of the same. We run Delly (discordant pairs plus split reads, with a paired-end scoring model that handles inversions and translocations better than most) and Smoove (Lumpy’s probabilistic combination of evidence types, genotyped by svtyper).

delly call -g GRCh38.fa -x human.hg38.excl.tsv -o delly.bcf sample.cram
bcftools view delly.bcf -Oz -o delly.vcf.gz && bcftools index delly.vcf.gz

smoove call --name sample --fasta GRCh38.fa \
  --exclude exclude.bed --genotype --processes 16 \
  --outdir smoove_out sample.cram

Then add a read-depth caller, because paired-end and split-read methods systematically lose large events. Once a deletion or duplication exceeds a few tens of kilobases, the informative read pairs are a vanishing fraction of the library, while the depth signal gets cleaner with size. This complementarity between depth-based and pair-based detection has been the organizing principle of SV discovery since the first population-scale efforts.¹ CNVpytor (the maintained successor to CNVnator) is our choice for its speed and its use of both depth and B-allele frequency:

cnvpytor -root sample.pytor -rd sample.cram
cnvpytor -root sample.pytor -his 1000 10000 100000
cnvpytor -root sample.pytor -partition 1000 10000 100000
cnvpytor -root sample.pytor -call 10000 > cnvpytor_10k.calls

Run the 10 kb bin size for confident large CNV calls and the 1 kb bin size if you want to reach smaller events, accepting far more noise. Depth callers give you copy number but poor breakpoints, typically accurate only to the bin size, so treat their coordinates as approximate and their copy-number estimates as the real output.

4. Merge, and treat the merge as a measurement

Merging is where naive pipelines go wrong. Two callers reporting “the same” deletion will place its breakpoints tens to hundreds of base pairs apart, and a merge that demands exact coordinates will report near-zero agreement while a merge that demands nothing will collapse distinct events. Use truvari collapse or SURVIVOR with explicit tolerances, and state them:

bcftools concat -a -Oz -o all_callers.vcf.gz \
  manta.diploidSV.vcf.gz delly.vcf.gz smoove_out/sample-smoove.genotyped.vcf.gz
bcftools sort all_callers.vcf.gz -Oz -o all_sorted.vcf.gz && tabix -p vcf all_sorted.vcf.gz

truvari collapse -i all_sorted.vcf.gz \
  -o merged.vcf -c collapsed.vcf \
  -r 500 -p 0.7 -P 0.7 --keep maxqual

Here -r 500 allows breakpoints to differ by 500 bp, -p 0.7 requires 70 percent reciprocal size similarity, and -P 0.7 requires 70 percent reciprocal overlap. Those are the parameters we use for short-read merges. If you tighten -r to 50 you will see agreement rates drop by a third for reasons that have nothing to do with whether the variants are real.

Now annotate each merged record with how many callers supported it and build tiers. Calls supported by all three pair/split callers with concordant genotypes are your high-confidence set, and in our experience they validate at well above 90 percent for deletions in unique sequence. Singleton calls from one caller are where nearly all the false positives live, and singleton duplications and inversions are worse than singleton deletions. Inversions in particular are a known weak point: they are usually flanked by inverted repeats, which means the breakpoints sit in exactly the sequence context short reads cannot resolve.³

5. Filter against population catalogs before you look at anything

Most of your SVs are common polymorphisms, and the fastest way to make the output tractable is to say so. The major structural variant alleles of the human genome have been characterized directly by long-read assembly of diverse genomes, which showed that a large fraction of SVs in any individual are recurrent alleles shared across populations rather than private events.⁴ Filtering to rare calls typically reduces a 7,000-call set to a few hundred.

bcftools annotate -a gnomad.v4.sv.sites.vcf.gz \
  -c INFO/AF:=INFO/AF -h af.hdr merged.vcf.gz -Oz -o merged.af.vcf.gz

SVAFotate BedProbs -v merged.af.vcf.gz -b SVAFotate_core_SV_popAFs.GRCh38.bed.gz \
  -o merged.svafotate.vcf -f 0.7

SVAFotate is more forgiving than bcftools annotate here because it does overlap-based matching with an adjustable reciprocal fraction, which is what you need given that your breakpoints will not match gnomAD’s exactly. Then annotate function:

AnnotSV -SVinputFile merged.svafotate.vcf \
  -genomeBuild GRCh38 -annotationMode both \
  -outputFile sample_annotsv.tsv

AnnotSV returns per-gene rows, a full-versus-split annotation mode so you can see both whole-SV and per-gene-overlap views, and an ACMG-style class from 1 to 5. Read that class as a triage aid rather than a verdict. Interpreting whether an SV affects health is clinical work. A call that lands in a known disease gene should go to a medical geneticist with the raw evidence attached. It should not go into a self-directed conclusion.

6. Look at the reads for every call you care about

No filter substitutes for looking. For each candidate you intend to take seriously, open the locus in IGV with the BAM loaded. Set “view as pairs” on and color alignments by insert size and pair orientation. A real heterozygous deletion shows a clean depth drop to roughly half, a cluster of read pairs with insert sizes larger than the library mode by the deletion length, and split reads whose clip points agree to the base. A tandem duplication shows depth increase plus pairs in everted (tail-to-head) orientation. An inversion shows two clusters of same-strand pairs, one at each breakpoint.

Artifacts have signatures too. If the “deletion” sits between two near-identical repeat copies and every supporting read has MAPQ 0, the event is probably a mismapping. If depth in the flanking regions is itself erratic, you are looking at a coverage artifact. These are common in GC-extreme regions and in the segmental duplications that cover roughly 5 percent of the genome and that generate a disproportionate share of both real and spurious CNV calls.⁵ For batch review, samplot renders IGV-like images non-interactively:

samplot plot -n sample -b sample.cram \
  -o chr7_del.png -c chr7 -s 100230000 -e 100236000 -t DEL

Reviewing 50 loci by eye takes about an hour and will change your opinion of your own pipeline more than any benchmark.

7. If you can, do this with long reads instead

Short-read SV calling has a ceiling set by read length, and the ceiling is low in exactly the regions where SVs cluster. Long reads from PacBio HiFi (15–20 kb, ~99.9 percent accurate) or Oxford Nanopore span most breakpoints and their flanking repeats, which turns breakpoint resolution from an inference problem into a read-alignment problem. Long-read studies routinely recover two to three times more SVs per genome than short-read pipelines, with the excess concentrated in insertions and in tandem-repeat regions where short reads are effectively blind.⁶ The characterization of common human SV alleles that makes population filtering possible was itself built on long-read assemblies for this reason.⁴

The pipeline is shorter:

minimap2 -ax map-hifi -t 16 --MD -Y GRCh38.fa hifi.fastq.gz \
  | samtools sort -@8 -o hifi.bam && samtools index hifi.bam

sniffles --input hifi.bam --vcf sniffles.vcf.gz \
  --reference GRCh38.fa --threads 16 --minsupport auto

cuteSV hifi.bam GRCh38.fa cutesv.vcf ./tmp \
  --max_cluster_bias_INS 1000 --diff_ratio_merging_INS 0.9 \
  --max_cluster_bias_DEL 1000 --diff_ratio_merging_DEL 0.5 -t 16

Run both Sniffles2 and cuteSV and merge with truvari collapse using tighter tolerances (-r 200) than you would for short reads, since long-read breakpoints agree far better. Sequence-resolved insertions come out of long reads directly, which is the single biggest gain: insertions are the SV class short reads miss most severely, because an inserted sequence longer than the read has no reference to align to.³

Common problems

Thousands of duplication calls concentrated near centromeres and acrocentric arms mean your exclude regions are wrong. Rebuild the BED from the UCSC centromere and gap tracks plus the ENCODE blacklist, and rerun. This is the most common cause of an SV callset that looks alarmingly large.

Very few calls, especially very few deletions under 500 bp, usually points to an insert-size problem or to an alignment that dropped supplementary records. Recheck samtools stats and the supplementary count from step 1 before blaming the caller.

Callers disagreeing on genotype for a variant they all detect is normal and not a bug. Genotyping an SV from short reads means estimating allele balance from a handful of spanning reads, and the estimate is noisy at 30x. If genotype matters to you, regenotype the merged sites uniformly with a dedicated tool (svtyper for pair-based sites, Paragraph or GraphTyper for graph-based regenotyping) rather than trusting per-caller genotypes.

Coordinates that disagree with a published variant by a few hundred base pairs are expected rather than wrong. Breakpoint precision varies by evidence type: split-read and assembly-based calls can be exact, while depth-based calls are only as precise as the bin size.⁷ Compare by overlap, never by exact position.

A large “deletion” that is really a reference error, or a call that appears in every sample you ever run, is a reference-assembly artifact. GRCh38 contains regions that differ from most human haplotypes, and the reference allele is sometimes the minor one. Population allele frequency near 1.0 in gnomAD-SV is the tell.

Finally, resist the urge to interpret. A structural variant that deletes exons of a gene you recognize is a measurement. Its consequence depends on the allele, the tissue, the transcript, and a great deal of clinical context. Take the call, the read-level evidence, and the annotation to a clinical geneticist and let them tell you what it means.

Footnotes

1. Can Alkan, Bradley P. Coe, Evan E. Eichler. Genome structural variation discovery and genotyping. Nature Reviews Genetics, 2011. https://doi.org/10.1038/nrg2958 ↩ ↩²

2. Shunichi Kosugi, Yukihide Momozawa, Xiaoxi Liu, et al. Comprehensive evaluation of structural variation detection algorithms for whole genome sequencing. Genome Biology, 2019. https://doi.org/10.1186/s13059-019-1720-5 ↩

3. Steve S. Ho, Alexander E. Urban, Ryan E. Mills. Structural variation in the sequencing era. Nature Reviews Genetics, 2019. https://doi.org/10.1038/s41576-019-0180-9 ↩ ↩²

4. Peter A. Audano, Arvis Sulovari, Tina A. Graves-Lindsay, et al. Characterizing the Major Structural Variant Alleles of the Human Genome. Cell, 2019. https://doi.org/10.1016/j.cell.2018.12.019 ↩ ↩²

5. Jeffrey M. Kidd, Tina Graves, Tera L. Newman, et al. A Human Genome Structural Variation Sequencing Resource Reveals Insights into Mutational Mechanisms. Cell, 2010. https://doi.org/10.1016/j.cell.2010.10.027 ↩

6. Wouter De Coster, Christine Van Broeckhoven. Newest Methods for Detecting Structural Variations. Trends in Biotechnology, 2019. https://doi.org/10.1016/j.tibtech.2019.02.003 ↩

7. Lorenzo Tattini, Romina D’Aurizio, Alberto Magi. Detection of Genomic Structural Variants from Next-Generation Sequencing Data. Frontiers in Bioengineering and Biotechnology, 2015. https://doi.org/10.3389/fbioe.2015.00092 ↩