GATK-gCNV enables the discovery of rare copy number variants from exome sequencing data

Abstract

Copy number variants (CNVs) are major contributors to genetic diversity and disease. While standardized methods, such as the genome analysis toolkit (GATK), exist for detecting short variants, technical challenges have confounded uniform large-scale CNV analyses from whole-exome sequencing (WES) data. Given the profound impact of rare and de novo coding CNVs on genome organization and human disease, we developed GATK-gCNV, a flexible algorithm to discover rare CNVs from sequencing read-depth information, complete with open-source distribution via GATK. We benchmarked GATK-gCNV in 7,962 exomes from individuals in quartet families with matched genome sequencing and microarray data, finding up to 95% recall of rare coding CNVs at a resolution of more than two exons. We used GATK-gCNV to generate a reference catalog of rare coding CNVs in WES data from 197,306 individuals in the UK Biobank, and observed strong correlations between per-gene CNV rates and measures of mutational constraint, as well as rare CNV associations with multiple traits. In summary, GATK-gCNV is a tunable approach for sensitive and specific CNV discovery in WES data, with broad applications.

Fig. 1

GATK-gCNV pipeline steps. a, Coverage information is collected from genome-aligned reads over a set of predefined genomic intervals. b, The original interval list is filtered to remove coverage outliers, unmappable genomic sequence, and regions of segmental duplications. c, Samples are clustered into batches based on read-depth profile similarity and each batch is processed separately. d, Chromosomal ploidies are inferred using total read-depth of each chromosome. e, The GATK-gCNV model learns read-depth bias and noise and iteratively updates copy number state posterior probabilities until a self-consistent state is obtained; after convergence, constant copy number segments are found using the Viterbi algorithm along with segmentation quality scores. Abbreviations: CN - copy number; QS - quality score.

Fig. 2.

Calling and benchmarking of GATK-gCNV callset in a cohort of more than 7,000 samples with matching deep WGS sequencing. a, A heatmap illustration of the distinct read-count signal of the 7,981 intervals chosen for the batch creation procedure. b, After normalizing for median read count, the first three PCs are clustered to determine which samples will be processed together with GATK-gCNV, colored by the assigned batch. c, For each of the 14 batches generated, a random subset of 200 samples was chosen to generate a read-count model using cohort-mode; the remaining samples were processed in case-mode. d, The recall (and e, precision) of rare CNVs in GATK-gCNV ES CNVs compared to WGS gold-standard CNVs as a function of the number of exons the variant spans. f, The recall (and g, precision) of de novo CNVs in GATK-gCNV compared to gold-standard WGS CNVs as a function of the number of exons. h, The recall (and i, precision) of rare CNVs in GATK-gCNV, XHMM, CONIFER, cn.mops, and ExomeDepth WES CNVs compared to WGS gold-standard CNVs as a function of the number of exons the variant spans. Abbreviations: PCA - principal component analysis; WES - exome sequencing WGS - whole genome sequencing.

Fig. 3.

A high-quality rare CNV callset was generated on 200,624 exomes from the UK Biobank (UKBB) using GATK-gCNV a, The variant-size distribution of high-quality, rare CNVs in the UKBB as a function of the number of exons each variant spans. b, The distribution of the number of rare, high-quality CNVs per-sample in the UKBB. c, Using 177,158 UKBB samples with matching CMA data, we find excellent validation of high-quality GATK-gCNV WES calls using Genome STRiP Intensity Rank Sum testing. d, GD CNV rates in the UKBB GATK-gCNV WES callset were highly concordant with rates from previous reports based on UKBB CMA data. e, The number of rare deletions observed over a gene in the UKBB GATK-gCNV callset is tightly correlated with LOEUF, with grey band representing LOESS smoothing of the 95% confidence intervals on corresponding point estimates. f, The number of rare duplications observed over a gene in the UKBB GATK-gCNV callset is also strongly correlated with the pTriplo score measuring intolerance to duplications, with grey band representing LOESS smoothing of the 95% confidence intervals on corresponding point estimates. g, The number of high-confidence duplications (IED) with both breakpoints within the boundaries of a gene are also correlated with LOEUF, with grey band representing LOESS smoothing of the 95% confidence intervals on corresponding point estimates. h, 16p11.2 deletions are associated with a significant increase in normalized BMI (n=41 carried a CN=1 deletion, n=61 carried a CN=3 duplication, and 169,711 individuals copy normal; boxplot corresponding to first, second, and third quartile of data, with whiskers denoting 1.5x interquartile range). i, PDZK1 deletions are associated with a significant decrease in normalized urate levels (n=145 carried a CN=1 deletion overlapping, n=143 carried a duplication of CN=3 overlapping, and 161,773 individuals copy norma; boxplot corresponding to first, second, and third quartile of data, with whiskers denoting 1.5x interquartile range l). j, CST3 duplications are significantly associated with decreased normalized eGFR values (n=6 carried a CN=3 duplication overlapping, n=3 carried a CN=3 duplication overlapping, and 162,666 individuals copy normal; boxplot corresponding to first, second, and third quartile of data, with whiskers denoting 1.5x interquartile range), on par with eGFR of individuals with renal failure (n=5,455). Abbreviations: CNV - copy number variation; DEL - deletion; DUP - duplication; CMA - chromosomal microarray; UKBB - UK Biobank; LOEUF - loss-of-function observed over expected upper bound fraction; pTriplo - probability of triplosensitivity; IED - intragenic exonic duplication; GD - genomic disorder; WES - exome sequencing; BMI - body mass index; eGFR - estimated glomerular filtration rate.