The mutational constraint spectrum quantified from variation in 141,456 humans

Abstract

Genetic variants that inactivate protein-coding genes are a powerful source of information about the phenotypic consequences of gene disruption: genes that are crucial for the function of an organism will be depleted of such variants in natural populations, whereas non-essential genes will tolerate their accumulation. However, predicted loss-of-function variants are enriched for annotation errors, and tend to be found at extremely low frequencies, so their analysis requires careful variant annotation and very large sample sizes¹. Here we describe the aggregation of 125,748 exomes and 15,708 genomes from human sequencing studies into the Genome Aggregation Database (gnomAD). We identify 443,769 high-confidence predicted loss-of-function variants in this cohort after filtering for artefacts caused by sequencing and annotation errors. Using an improved model of human mutation rates, we classify human protein-coding genes along a spectrum that represents tolerance to inactivation, validate this classification using data from model organisms and engineered human cells, and show that it can be used to improve the power of gene discovery for both common and rare diseases.

Fig. 1. Aggregation of 141,456 exome and genome sequences.

a, Uniform manifold approximation and projection (UMAP)^, plot depicting the ancestral diversity of all individuals in gnomAD, using seven principal components. Note that long-range distances in the UMAP space are not a proxy for genetic distance. b, The number of individuals by population and subpopulation in the gnomAD database. Colours representing populations in a and b are consistent. c, d, The mutability-adjusted proportion of singletons (MAPS) is shown across functional categories for SNVs in exomes (c; x axis shared with e and g) and genomes (d; x axis shared with f and h). Higher values indicate an enrichment of lower frequency variants, which suggests increased deleteriousness. e, f, The proportion of possible variants observed for each functional class for each mutational type for exomes (e) and genomes (f). CpG transitions are more saturated, except where selection (for example, pLoFs) or hypomethylation (5′ untranslated region) decreases the number of observations. g, h, The total number of variants observed in each functional class for exomes (g) and genomes (h). Error bars in c–f represent 95% confidence intervals (note that in some cases these are fully contained within the plotted point).

Fig. 2. Generating a high-confidence set of pLoF variants.

a, The percentage of variants filtered by LOFTEE grouped by ClinVar status and gnomAD frequency. Despite not using frequency information, LOFTEE removes a larger proportion of common variants, and a very low proportion of reported disease-causing variation. b, MAPS (see Fig. 1c, d) is shown by LOFTEE designation and filter. Variants filtered out by LOFTEE exhibit frequency spectra that are similar to those of missense variants; predicted splice variants outside the essential splice site are more rare, and high-confidence variants are very likely to be singletons. Only SNVs with at least 80% call rate are included here. Error bars represent 95% confidence intervals. c, d, The total number of pLoF variants (c), and proportion of genes with more than ten pLoF variants (d) observed and expected (in the absence of selection) as a function of sample size (downsampled from gnomAD). Selection reduces the number of variants observed, and variant discovery approximately follows a square-root relationship with the number of samples. At current sample sizes, we would expect to identify more than 10 pLoF variants for 72.1% of genes in the absence of selection.

Fig. 3. The functional spectrum of pLoF impact.

a, The percentage of genes in a set of curated gene lists represented in each LOEUF decile. Haploinsufficient genes are enriched among the most constrained genes, whereas recessive genes are spread in the middle of the distribution, and olfactory receptor genes are largely unconstrained. b, The occurrence of 6,735 rare LoF deletion structural variants (SVs) is correlated with LOEUF (computed from SNVs; linear regression r = 0.13; P = 9.8 × 10⁻⁶⁸). Error bars represent 95% confidence intervals from bootstrapping. c, d, Constrained genes are more likely to be lethal when heterozygously inactivated in mouse and cause cellular lethality when disrupted in human cells (c), whereas unconstrained genes are more likely to be tolerant of disruption in cellular models (d). For all panels, more constrained genes are shown on the left.

Fig. 4. Biological properties of constrained genes and transcripts.

a, The mean number of protein–protein interactions is plotted as a function of LOEUF decile: more constrained genes have more interaction partners (LOEUF linear regression r = −0.14; P = 1.7 × 10⁻⁵¹). Error bars correspond to 95% confidence intervals. b, The number of tissues where a gene is expressed (transcripts per million > 0.3), binned by LOEUF decile, is shown as a violin plot with the mean number overlaid as points: more constrained genes are more likely to be expressed in several tissues (LOEUF linear regression r = −0.31; P < 1 × 10⁻¹⁰⁰). c, For 1,740 genes in which there exists at least one constrained and one unconstrained transcript, the proportion of expression derived from the constrained transcript is plotted as a histogram.

Fig. 5. Disease applications of constraint.

a, The rate ratio is defined by the rate of de novo variants (number per patient) in 5,305 cases of intellectual disability/developmental delay (ID/DD) divided by the rate in 2,179 controls. pLoF variants in the most constrained decile of the genome are approximately 11-fold more likely to be found in cases compared to controls. Error bars represent 95% confidence intervals. b, Marginal enrichment in per-SNV heritability explained by common (minor allele frequency > 5%) variants within 100-kb of genes in each LOEUF decile, estimated by linkage disequilibrium (LD) score regression. Enrichment is compared to the average SNV genome-wide. The results reported here are from random effects meta-analysis of 276 independent traits (subsetted from the 658 traits with UK Biobank or large-scale consortium GWAS results). Error bars represent 95% confidence intervals. c, Conditional enrichment in per-SNV common variant heritability tested using regression of linkage disequilibrium score in each of 658 common disease and trait GWAS results. P values evaluate whether per-SNV heritability is proportional to the LOEUF of the nearest gene, conditional on 75 existing functional, linkage disequilibrium, and minor-allele-frequency-related genomic annotations. Colours alternate by broad phenotype category.

Extended Data Fig. 1. Overview of the sample quality control workflow.

a, Exome (square) and genome (circle) samples underwent quality control in the following stages: hard filtering (step 1), relatedness inference (step 2), ancestry inference (step 3), platform inference (step 4, for exomes only), and population- and platform-specific outlier filtering (step 5). See Supplementary Information for further details. Except for samples failing hard filters (dotted outline), all quality control analyses were applied to all samples, regardless of the presence or absence of other quality control flags (such as relatedness, lack of release permissions, or outlier status; red diagonal bar). Assignment of ancestry labels is represented by fill colour and accompanying three-letter ancestry group abbreviation. Assignment of platform labels is represented by outline colour and a numbered label for exomes (corresponding to imputed platforms) and a PCR ± label for genomes. The final set of samples included in the gnomAD release (125,748 exomes and 15,708 genomes) was defined to be the set of unrelated samples with release permissions, no hard filter flags, and no population- and platform-specific outlier metrics (step 6). b, In exomes, the chromosomal sex of samples was inferred based on the inbreeding coefficient on chromosome X and the coverage of chromosome Y into male (green), female (amber), ambiguous sex (pink), and sex chromosome aneuploid (blue). c, The top two principal components from PCA-HDBSCAN analysis of exome capture regions. Sequencing platforms were inferred for exome samples based on principal component analysis of biallelic variant call rates over all known exome capture regions, and samples were assigned a cluster label (0–15, or unknown) using HDBSCAN. d, We performed platform- and population-specific outlier filtering for several quality-control metrics. The distribution of the number of deletions in samples from south Asian individuals across platforms is shown. Distributions (and accordingly, median and median absolute deviations) for these metrics varied widely both by population and sequencing platform (numbered on the y axis). Outliers (black dots) were defined as samples with values outside four median absolute deviations (shown by dotted vertical lines) from the median of a given metric.

Extended Data Fig. 2. Variant calling performance for common variants.

a–h, Precision-recall curves are shown for variant calls in two samples with independent gold-standard data, NA12878 (a–d) and a synthetic diploid mixture (e–h). The random forest (blue) approach described here is compared to the current state-of-the-art GATK variant quality score recalibration (orange) for exome SNVs (a, e) and indels (b, f), and genome SNVs (c, g) and indels (d, h). Note that the indels presented in f and h exclude 1-base-pair (bp) indels as they are not well characterized in the synthetic diploid mixture gold standard sample. In all cases, at the thresholds chosen (dashed lines representing 10% and 20% of SNVs and indels filtered, respectively), random forest outperforms or is similar to variant quality score recalibration.

Extended Data Fig. 3. Variant calling performance for rare variants.

a–j, The x axes show the cumulative ranked percentile for our random forest (blue) model and, as a comparison, for the current state-of-the-art GATK variant quality score recalibration (orange). That is, the point at 10 shows the performance of the 10% best-scored data; the point at 50 shows the performance 50% best-scored data. a–d, The number of transmitted singletons (singletons in the unrelated individuals that are transmitted to an offspring) on chromosome 20 for exome SNVs (a) and indels (b), and genome SNVs (c) and indels (d). Chromosome 20 was not used for training our random forest model. We expect most of these to be real variants because we observe Mendelian transmission of an allele that was sequenced independently in a parent and child. e–h, The number of bi-allelic de novo calls per child (4,568 exomes, 212 genomes) outside of low-complexity regions. The expectation is that there is approximately 1.6 de novo SNV (e) and 0.1 de novo indels per exome (f), and 65 de novo SNVs (g) and 5 de novo indels (h) per genome. i, j, The number of independently validated de novo mutations, available for a subset of 331 exome samples for which de novo mutations were validated as part of other studies. In all cases, at the thresholds chosen (dashed lines representing 10% and 20% of SNVs and indels filtered, respectively), random forest outperforms or is similar to variant quality score recalibration.

Extended Data Fig. 4. Variant discovery at large sample sizes.

a, b, The total number of variants observed (a) and the proportion of possible variants observed (b) as a function of sample size, broken down by variant class. At large sample sizes, CpG transitions become saturated, as previously described. Colours are consistent in a and b. c, This results in a decrease of the transition/transversion (Ti/Tv) ratio. d, When broken down by functional class, we observe the effects of selection, in which synonymous variants have the highest proportion observed, followed by missense and pLoF variants. e, f, The number of additional pLoF variants introduced into the cohort as a function of sample size on a log (e) and linear (f) scale. Here, gnomAD (black) refers to a uniform sampling from the population distribution of the full cohort of exome-sequenced individuals.

Extended Data Fig. 5. Using LOFTEE to create a high-confidence set of pLoF variation.

a, Schematic of LOFTEE filters. LOFTEE filters out putative stop-gained, essential splice, and frameshift variants based on sequence and transcript context, as well as flagging exonic features such as conservation (not shown). For instance, variants that are not predicted to disrupt splicing based on retention of a strong splice site, or rescue of a nearby splice site. Additional filters not shown include: ANC_ALLELE (the alternative allele is the ancestral allele), NON_ACCEPTOR_DISRUPTING and DONOR_RESCUE (opposite to those already shown). b, To tune the END_TRUNC filter, we retained variants that pass the 50-bp rule (are more than 50 bp before the 3′-most splice site). The overall MAPS score for variants that fail this rule is shown in grey. For the remaining 39,072 variants, we computed the sum of the genomic evolutionary rate profiling (GERP) score of bases deleted by the variant. At 40 bins of this score, we compute the MAPS score for those variants retained at this threshold (red) compared to variants removed at this threshold (blue), and plot this as a function of the proportion of variants filtered at this threshold. We chose the 50% point as it retains variants with a MAPS score of 0.14, while removing variants with a MAPS score of 0.06. Error bars represent 95% confidence intervals. c, Density plot of aggregate pLoF frequency computed from high-confidence pLoF variants discovered using LOFTEE.

Extended Data Fig. 6. Computing the depletion of variation of functional categories.

a, The distribution of mean methylation values across 37 tissues and across every CpG dinucleotide in the genome. We divided the genome into 3 levels (low methylation, missing or < 0.2; medium, 0.2–0.6; and high, >0.6) and computed all ensuing metrics based on these categories. b, Comparison of estimates of the mutation rate with previous estimates. For transversions and non-CpG transitions, we observe a strong correlation (linear regression r = 0.98; P = 2.6 × 10⁻⁶⁵). For CpG transitions, the new estimates are calculated separately for the three levels of methylation and track with these levels. Colours and shapes are consistent in b–d. c, For c–e, only synonymous variants are considered. The proportion of possible variants observed for each context is correlated with the mutation rate. We compute two fit lines, one for CpG transitions, and one for other contexts to calibrate our estimates. d, Calibration of each context to compute a predicted proportion observed after fitting the two models in c, which is used to calculate an expected number of variants at high coverage. e, With an expectation computed from high coverage regions, the observed/expected ratio follows a logarithmic trend with the median coverage below 40×, which is used to correct low coverage bases in the final expectation model. f–h, For each transcript, the observed number of variants is plotted against the expected number from the model described above, for synonymous (f), missense (g), and pLoF (h) variants, and the linear regression coefficient is shown. Note that the expectation does not include selection, and so, pLoF and, to a lesser extent, missense variants exhibit lower observed values than expected.

Extended Data Fig. 7. Genomic properties of constrained genes.

a, b, Histogram of the observed/expected ratio of pLoF variation (a) and LOEUF (b). Most genes have fewer observed variants than expected (median observed/expected = 0.48), and the genes with no observed pLoFs are distinguished between confidently constrained genes and noise by LOEUF. c, A 2D density plot of the number of observed versus expected pLoF variants. The boundaries of each decile are plotted as gradients (that is, the most constrained decile is below the lowest red line). d, The LOEUF of a gene is correlated with its coding sequence length (beta = −1.07 × 10⁻⁴; P < 10⁻¹⁰⁰): thus, for all downstream statistical tests, we adjust for gene length or remove genes with fewer than 10 expected pLoFs. e, Observed/expected ratios of various functional classes across genes within each LOEUF decile. The most constrained decile has approximately 6% of the expected pLoFs, while synonymous variants are not depleted and missense variants exhibit modest depletion. f, The percentage of each LOEUF decile that was described in ExAC as constrained, or pLI > 0.9. g, The percentage of each LOEUF decile that have at least one homozygous pLoF variant. h, Box plots of the aggregate pLoF frequency for each LOEUF decile. Centre line denotes the median; box limits denote upper and lower quartiles; whiskers denote 1.5× the interquartile range; points denote outliers). In e–g, error bars represent 95% confidence intervals (note that in some cases these are fully contained within the plotted point).

Extended Data Fig. 8. Biological properties of constrained genes.

a, The percentage of genes in each functional category from Pharos (see Supplementary Information) is broken down by the LOEUF decile. b, The mean number of tissues in which a transcript is expressed, binned by transcript-based LOEUF decile, is shown for all transcripts and canonical transcripts. c, The percentage of genes in which the most expressed transcript is also the most constrained is plotted in red, which is enriched compared to a permuted set (blue). d, For 927 genes with expected pLoF ≥10 in both the African/African American and European population subsets (n = 8,128), the LOEUF scores are highly correlated (linear regression r = 0.78, P < 10⁻¹⁰⁰), with a lower mean score observed in the African/African American population (0.49 versus 0.62; two-sided t-test P = 4.1 × 10⁻¹⁴), which has a higher effective population size. e, The mean LOEUF score for 865 genes with expected pLoF ≥ 10 in all populations (n = 8,128). Error bars represent 95% confidence intervals.

Extended Data Fig. 9. Applications of constraint metrics to rare variant analysis of disease.

a, Proportion of each LOEUF decile found in OMIM. b, Proportion of disease-associated genes discovered by whole-exome/genome sequencing (WES/WGS) compared to conventional (typically linkage) methods, plotted by LOEUF decile. The former are more constrained (LOEUF 0.674 versus 0.806, two-sided t-test P = 1.2 × 10⁻¹⁶), which suggests that these techniques are more effective for picking up genes with a de novo mechanism of disease, compared to recessive genes identified by linkage methods. c, Similar to Fig. 5a, the rate ratio is defined by the rate of de novo variants (number per patient) in autism cases divided by the rate in controls. pLoF variants in the most constrained decile of the genome are approximately fourfold more likely to be found in cases compared to controls. d, The mean odds ratio of a logistic regression of schizophrenia is plotted for each LOEUF decile. Error bars in a–d correspond to 95% confidence intervals.

Extended Data Fig. 10. Applications of constraint metrics to common variant analysis of disease.

a, The ${\hat{\mathrm{\tau }}}^{⁎}$ coefficient (see Supplementary Information) for each LOEUF decile across 276 independent traits. Unlike the enrichment measure reported in Fig. 5, ${\hat{\mathrm{\tau }}}^{⁎}$is adjusted for 74 baseline genomics annotations. Positive values of ${\hat{\mathrm{\tau }}}^{⁎}$ indicate greater per-SNP heritability than would be expected based on the other annotations in the baseline model, whereas negative values indicate depleted per-SNP heritability compared to that baseline expectation. b, Enrichment coefficient for each LOEUF decile using different window sizes to define which SNPs to include upstream and downstream of each gene. c, Enrichment coefficient for each LOEUF decile across traits after controlling for brain expression and gene size. Results are consistent with those shown in Fig. 5, which indicates that brain gene expression and gene size do not fully explain the enrichment of heritability observed in constrained genes. Error bars represent 95% confidence intervals.