Building a CNV Analysis Pipeline from Your Own Whole-Genome Data

By the end of this guide you will have a filtered, annotated set of copy number variants (CNVs) from your own whole-genome sequence: deletions and duplications ranging from a few hundred base pairs to whole chromosome arms, each with a copy number estimate, a quality score, a population frequency, and a list of overlapping genes. You need a 30x or deeper short-read WGS alignment (a CRAM or BAM against GRCh38, ideally the analysis set with alt contigs excluded from alignment), the matching reference FASTA and index, roughly 200 GB of scratch disk, 16 or more CPU cores, and 32 GB of RAM. Everything below runs on Linux with conda or Docker. A single sample takes two to six hours of wall clock time depending on how many callers you run.

A word on vocabulary before the commands. A single nucleotide polymorphism (SNP) is a one-base change at a fixed position. A CNV is a change in how many copies of a stretch of DNA you carry, from a 500 bp deletion inside one intron to a 2 Mb duplication spanning a dozen genes. The two are detected by completely different signals. SNPs come from the base calls in aligned reads. CNVs come from four kinds of evidence. Read depth is how many reads pile up in a window. Read pairs show how far apart mates land relative to the expected insert size, and split reads cross a breakpoint. B-allele frequency is the allele balance of nearby heterozygous SNPs. Short-read WGS can detect CNVs well, and better than exome or array data, because coverage is uniform across the genome and breakpoints often fall in mappable sequence. It detects them imperfectly, which is the reason this pipeline uses more than one caller.

1. Check that your alignment can support CNV calling

Depth-based CNV calling is a measurement of read counts per window, so it inherits every bias in library prep and alignment. Before calling anything, quantify those biases, because they set the floor on what you can resolve.

# Mean depth, coverage uniformity, insert size distribution
mosdepth --by 1000 --fast-mode -t 4 sample sample.cram \
  --fasta GRCh38_full_analysis_set_plus_decoy_hla.fa

samtools stats -@ 8 --reference GRCh38.fa sample.cram > sample.stats
grep -E "^SN" sample.stats | grep -E "insert size|reads mapped|reads duplicated"

Three numbers matter. Mean autosomal depth should be at least 30x. Below about 20x, small deletions between 1 and 10 kb start dropping out of read-depth callers, because the count difference between two copies and one copy in a 1 kb window is within Poisson noise. The duplicate rate should be under roughly 15 percent, since PCR duplicates inflate depth non-uniformly and PCR-free libraries are noticeably cleaner for this purpose. The insert size standard deviation should be tight relative to the mean (a mean near 400 bp with an SD of 80 to 120 bp is typical), because read-pair callers set their anomalous-insert threshold from that distribution and a fat tail turns into false deletion calls.

Also look at the per-1 kb coverage from mosdepth and compute the median absolute deviation across autosomes. If coverage swings wildly with GC content, you will need the GC correction step in section 3. Expect the caller to be conservative at telomeres, centromeres, and segmental duplications regardless.

2. Call structural variants from read pairs and split reads

Start with a caller that uses breakpoint evidence rather than depth, because those calls come with base-pair-resolution breakpoints and you can validate them by eye. Manta is the practical choice for germline short-read data: it is fast, it emits well-formed VCF, and its deletion and duplication calls carry paired-read and split-read counts you can filter on.

configManta.py \
  --bam sample.cram \
  --referenceFasta GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --runDir manta_out \
  --callRegions calling_regions.bed.gz

manta_out/runWorkflow.py -m local -j 16

The --callRegions file should be a bgzipped, tabix-indexed BED of the primary chromosomes with centromeres and the hg38 ENCODE blacklist subtracted. Without it Manta spends a large fraction of its runtime producing calls in regions you will discard anyway.

The useful output is manta_out/results/variants/diploidSV.vcf.gz. Split it into CNV-like classes and drop the low-evidence tail:

bcftools view -f PASS -i 'INFO/SVTYPE="DEL" || INFO/SVTYPE="DUP"' \
  manta_out/results/variants/diploidSV.vcf.gz \
  -Oz -o sample.manta.cnv.vcf.gz
bcftools index sample.manta.cnv.vcf.gz

Manta’s sensitivity drops sharply above roughly 100 kb for duplications, because a long duplication has few reads spanning its breakpoints relative to its span, and tandem duplications in repetitive sequence produce ambiguous mate orientations. That is exactly the size range where read depth works well, which is why the next step is not optional.

3. Call CNVs from read depth

Read-depth calling divides the genome into windows and normalizes the count in each window against GC content, mappability, and ideally a panel of other samples sequenced the same way. It then segments the normalized log2 ratio into regions of constant copy number. For a single germline genome without a matched panel, CNVpytor is the tool we reach for. It reads CRAM directly, does its own GC correction, and combines depth with B-allele frequency from your SNP calls. That resolves copy-neutral loss of heterozygosity and helps distinguish a real one-copy deletion from a mapping artifact.

# 1. Read depth into a .pytor file, 10 kb and 100 kb bins
cnvpytor -root sample.pytor -rd sample.cram -T GRCh38.fa
cnvpytor -root sample.pytor -his 1000 10000 100000
cnvpytor -root sample.pytor -partition 1000 10000 100000
cnvpytor -root sample.pytor -call 10000 > sample.cnvpytor.10k.tsv
cnvpytor -root sample.pytor -call 1000  > sample.cnvpytor.1k.tsv

# 2. Add B-allele frequency from an existing SNV VCF
cnvpytor -root sample.pytor -snp sample.snv.vcf.gz -sample SAMPLE_ID
cnvpytor -root sample.pytor -mask_snps
cnvpytor -root sample.pytor -baf 10000
cnvpytor -root sample.pytor -call combined 10000 > sample.cnvpytor.combined.tsv

The call table gives type, region, copy number estimate, and four statistics per event. Filter on them rather than trusting raw output. We require the e-value (column 5, the significance of the depth shift) below 1e-5. The fraction of the region that is in the mappability mask (column 9) must be below 0.5, and the ratio of within-region depth variance to genome-wide variance below about 0.2. In practice that removes most of the garbage at pericentromeric boundaries.

If your data came from a pipeline that already produced a depth-based CNV VCF, keep it. Vendor pipelines such as DRAGEN run a self-normalizing depth caller with a GC-corrected panel of normals and generally produce good large-event calls, and modern GPU-accelerated WGS pipelines have made running the full germline workflow including structural variants cheap enough that there is no reason to skip it on cost grounds.¹ Having a second depth-based call set is useful precisely because concordance between independent callers is the most informative quality signal available to you.

4. Merge the call sets and treat concordance as evidence

Different algorithm families find different variants, and the overlap between them is smaller than most people expect. In monozygotic twin quartets, discordance between genetically identical individuals is by construction a false positive. Work using them found that sequencing-based CNV callers varied widely in both sensitivity and precision, and that agreement between methods was limited.² A somatic benchmarking study on a hyper-diploid cancer genome reached a similar conclusion: performance depended heavily on event size, copy number state, and the specific tool, with no single caller dominating.³ The practical response is to run two or three callers from different families and rank calls by how many support them.

# Normalize both call sets to a common BEDPE-ish format, then merge
svtk vcf2bed sample.manta.cnv.vcf.gz manta.bed
awk 'BEGIN{OFS="\t"} $1!~/^#/ {print $2,$3,$4,$1,$6,$4-$3}' \
  sample.cnvpytor.combined.tsv | tr ':-' '\t\t' > cnvpytor.bed

# SURVIVOR for VCF-level merging with a 1 kb breakpoint tolerance
ls sample.manta.cnv.vcf.gz sample.cnvpytor.vcf > vcf_list.txt
SURVIVOR merge vcf_list.txt 1000 1 1 0 0 50 sample.merged.vcf

The arguments to SURVIVOR merge come in a fixed order. They are maximum breakpoint distance of 1000 bp, minimum number of supporting callers to output of 1, and require same type of 1. The remaining three are require same strand of 1, estimate distance from SV size of 0, and minimum SV size of 50. Setting the minimum supporting callers to 1 keeps singletons so you can filter later. Keep singletons in the file and record the supporting-caller count in the INFO field. A large duplication called only by depth is expected rather than suspicious, while a 2 kb deletion called only by depth deserves scrutiny.

For interval-level work with the resulting BED files, bedtools intersect -f 0.5 -r -a manta.bed -b cnvpytor.bed gives reciprocal 50 percent overlap, which is the standard concordance criterion in the CNV literature and the one to use when you compare your calls to any published set.

5. Annotate with population frequency and gene overlap

A raw CNV call means almost nothing until you know how common it is. A typical human genome carries a few thousand CNVs larger than 1 kb, and the overwhelming majority are common polymorphisms present in a substantial fraction of the population. Annotation is what reduces thousands of calls to the handful worth reading.

# AnnotSV: gene overlap, dosage sensitivity, population AF, ACMG-style scoring
AnnotSV -SVinputFile sample.merged.vcf \
  -genomeBuild GRCh38 \
  -annotationMode both \
  -outputFile sample.annotsv.tsv \
  -svtBEDcol 4 \
  -overlap 70

AnnotSV joins your calls against several databases. It uses gnomAD-SV and the Database of Genomic Variants for frequency and ClinGen dosage sensitivity ratings for gene-level haploinsufficiency and triplosensitivity. It also uses OMIM for gene-disease associations and RefSeq/ENSEMBL for exon-level overlap. The full rows describe the whole event, the split rows describe the event’s intersection with each gene, and the split rows are what you read.

Then filter aggressively. We keep calls that are rare (gnomAD-SV allele frequency below 1 percent), overlap at least one exon of a gene with a ClinGen haploinsufficiency score of 3 for deletions or triplosensitivity score of 3 for duplications, or exceed 500 kb regardless of gene content. Everything else goes into a second-tier file you can return to with a specific question.

For reading the results, a tabular viewer beats a spreadsheet. CNVizard is a Streamlit application built for interactive review of annotated CNV call sets, with filtering by gene panel, frequency, and call quality, and it is light enough to run on a laptop against your own files.⁴ If you prefer to write your own analysis, load the AnnotSV TSV into a pandas DataFrame and treat it as a join problem, which is what it is.

6. Visually confirm anything you intend to act on

Every CNV you care about should be looked at, because the failure modes of depth callers are visually obvious and statistically subtle. Load the region in IGV with the CRAM and a coverage track, then check three things. First, does the coverage step change at the reported breakpoints and stay changed across the interval. Second, do mate-pair orientations and insert sizes at the breakpoints match the called type, meaning larger-than-expected inserts for a deletion and reversed orientation for a tandem duplication. Third, does the B-allele frequency of heterozygous SNPs inside a called deletion collapse toward 0 and 1 as it must if one haplotype is gone.

# Static plots for a batch of regions
samplot plot -n SAMPLE -b sample.cram \
  -c chr16 -s 28800000 -e 29100000 -t DEL \
  -o chr16_28.8M_DEL.png --reference GRCh38.fa

If the coverage drop is gradual rather than stepped, you are probably looking at GC or mappability bias. If the interval sits inside a segmental duplication, the call may be a mismapping artifact, and short reads cannot settle it. Long-read sequencing or optical mapping can, and a clinical laboratory using targeted MLPA or a validated array is the right place to confirm a specific event.

7. Make the pipeline reproducible

A CNV call set is only interpretable if you can say exactly how it was produced, which means pinning tool versions, reference build, bin sizes, and filter thresholds in a workflow file rather than a shell history. Snakemake and Nextflow both do this well, and clinical variant pipelines built on Snakemake with containerized tools are a reasonable model to copy for structure, even when the biology in question differs from yours.⁵ Write the whole thing as rules with explicit inputs, outputs, and a conda or container directive per rule, keep the reference build and bin sizes in a single config YAML, and record the tool versions in the output directory. When you re-run this in a year on deeper data, you will want to compare against the old call set, and that comparison is only meaningful if both runs are specified.

Common problems

Depth callers over-call at the boundaries of low-mappability regions. Any window whose reads mostly cannot be placed uniquely gives a depth estimate that is not a copy number estimate. Apply an exclusion mask before segmentation rather than after, because a bad window in the middle of a segment drags the whole segment’s log2 ratio. The mask should be the hg38 ENCODE blacklist plus centromeres, telomeres, and segmental duplications above 98 percent identity.

Bin size determines what you can see, and there is no single right choice. At 1 kb bins you can resolve small exonic deletions but the per-bin count noise generates many false segments. At 100 kb bins the segmentation is clean but you miss everything below about 50 kb. Run both and merge, treating the small-bin calls as a lower-confidence tier. The segmentation algorithm matters too: methods based on genome alteration detection analysis were designed to be fast and flexible enough to segment thousands of samples for association work, and the choice of segmentation and its sparsity parameter changes the number of calls substantially.⁶

Do not compare your WGS calls directly against array-based CNV calls without adjusting expectations. Array platforms detect CNVs from probe intensity at fixed probe positions, so their breakpoints are coarse and their sensitivity depends on probe density. Different array algorithms also disagree with each other considerably even on the same data.⁷ A WGS call that an old array missed is not necessarily wrong, and an array call your WGS pipeline missed is often in a region your mask excluded.

Common CNVs are not a bug in your pipeline. Plenty of large duplications and deletions are ordinary polymorphisms carried by a significant fraction of people. Copy number variation has been examined as a contributor to complex traits with results that are often modest. One example is genome-wide work on age-related macular degeneration subtypes, where heritability attributable to common variation was estimated but subtype differences were hard to resolve.⁸ Rarity plus dosage sensitivity plus phenotype is the filter that matters, and the phenotype half of that is not something a pipeline can supply.

Finally, a CNV call is a measurement. It does not amount to a conclusion. Copy number changes are the mechanism behind a number of well-described conditions. The interpretation of a specific rare deletion or duplication in a specific person requires a clinical genetics evaluation, family history, and confirmation in an accredited laboratory. If your pipeline surfaces a rare deletion overlapping a haploinsufficient disease gene, the correct next step is to take the coordinates to a clinical geneticist or genetic counselor. Reading more papers about the gene is not the next step. Tools that summarize genome-wide patterns rather than individual events, such as pipelines for oncogenomic signature analysis, are a reminder that the same call set answers different questions at different scales. The individual-variant question is the one that needs a clinician.⁹

Footnotes

1. Stefano Marangoni, Federica Furia, Debora Charrance, et al. GeNePi: a graphics processing unit enhanced next-generation bioinformatics pipeline for whole-genome sequencing analysis. Briefings in Bioinformatics, 2026. https://doi.org/10.1093/bib/bbag001 ↩

2. Marc-André Legault, Simon Girard, Louis-Philippe Lemieux Perreault, et al. Comparison of Sequencing Based CNV Discovery Methods Using Monozygotic Twin Quartets. PLOS ONE, 2015. https://doi.org/10.1371/journal.pone.0122287 ↩

3. Daniall Masood, Luyao Ren, Cu Nguyen, et al. Evaluation of somatic copy number variation detection by NGS technologies and bioinformatics tools on a hyper-diploid cancer genome. Genome Biology, 2024. https://doi.org/10.1186/s13059-024-03294-8 ↩

4. Jeremias Krause, Carlos Classen, Daniela Dey, et al. CNVizard—a lightweight streamlit application for an interactive analysis of copy number variants. BMC Bioinformatics, 2024. https://doi.org/10.1186/s12859-024-06010-2 ↩

5. Jingyu Yang, Tim Beißbarth, Jürgen Dönitz. Onkopipe: A Snakemake Based DNA-Sequencing Pipeline for Clinical Variant Analysis in Precision Medicine. Studies in Health Technology and Informatics, 2023. https://doi.org/10.3233/shti230694 ↩

6. Roger Pique-Regi, Alejandro Cáceres, Juan R González. R-Gada: a fast and flexible pipeline for copy number analysis in association studies. BMC Bioinformatics, 2010. https://doi.org/10.1186/1471-2105-11-380 ↩

7. L. Winchester, C. Yau, J. Ragoussis. Comparing CNV detection methods for SNP arrays. Briefings in Functional Genomics and Proteomics, 2009. https://doi.org/10.1093/bfgp/elp017 ↩

8. Lucia Sobrin, Stephan Ripke, Yi Yu, et al. Heritability and Genome-Wide Association Study to Assess Genetic Differences between Advanced Age-related Macular Degeneration Subtypes. Ophthalmology, 2012. https://doi.org/10.1016/j.ophtha.2012.03.014 ↩

9. Marta Portasany-Rodríguez, Gonzalo Soria-Alcaide, Elena G Sánchez, et al. PULPO: pipeline of understanding large-scale patterns of oncogenomic signatures. Bioinformatics, 2026. https://doi.org/10.1093/bioinformatics/btag118 ↩