Main

Exome sequencing has enabled the discovery of rare coding variants, and has thus provided insights into gene function that have accelerated the pace of disease-associated gene discovery across Mendelian and common disorders1,2,3,6,9,10,11,12. Furthermore, exome sequencing has identified protective alleles that highlight drug targets that could be amenable to pharmacological intervention2,13,14,15,16,17. For example, anti-PCSK9 drug therapy is based on the observation that a loss of PCSK9 function is associated with reduced levels of cholesterol18.

Cataloguing rare coding variation can help with the implementation of precision medicine19,20. Large datasets of genetic variation that are representative of the human population are essential for the comprehensive discovery and interpretation of rare variants. Roadmaps for numerous large-scale sequencing studies have been proposed, and several efforts are now underway21,22,23,24. The Genome Aggregation Database4 (gnomAD) and Trans-Omics for Precision Medicine8 (TOPMed) initiatives have developed large public databases of genetic variation derived from approximately 200,000 individuals and 132,000 individuals, respectively. Here, we describe a harmonized collection of exonic data derived from 983,578 individuals who represent a diverse array of ancestries. We calculate continental and fine-scale ancestry-based allele frequencies across this dataset and make the data publicly available through the RGC research browser: https://rgc-research.regeneron.com/me.

Survey of variation in the RGC-ME dataset

The Regeneron Genetics Center Million Exome (RGC-ME) dataset contains the genetic variation observed in 983,578 individuals. These data span dozens of collaborations, including large biobanks and health systems. All data were generated by Regeneron Genetics Center using a single harmonized sequencing and informatics protocol. Previously published datasets, such as the UK Biobank and the Mexico City Prospective Study, were reprocessed9,25. The RGC-ME dataset comprises both outbred and founder populations spanning African (AFR), European (EUR), East Asian (EAS), Indigenous American (IAM), Middle Eastern (MEA) and South Asian (SAS) continental ancestries, and includes cohorts with relatively high rates of consanguinity. More than 190,000 of the unrelated participants (23%) are of non-EUR ancestry in the RGC-ME dataset, as compared with 35,000 in gnomAD genomes (v.3.1.2), 53,000 in gnomAD exomes (v.2.1.1), and 91,000 in TOPMed Freeze 8, indicating that RGC-ME represents a large increase in the number of individuals of non-EUR ancestry in datasets of genetic variation4,8 (Fig. 1a and Supplementary Table 1a).

Fig. 1: Variant survey and population counts in the RGC-ME dataset.
figure 1

a, Summed proportional ancestry (sum of weighted ancestry probabilities) at continental, sub-continental and regional levels for 821,979 unrelated samples (Supplementary Table 1b). All subsequent variant counts and surveys have been performed in the unrelated analysis set. UNK, unknown. b, Count of variants unique to RGC-ME (that is, variants not in gnomAD v.3.1.2 genomes, gnomAD v.2.1.1 exomes and TOPMed Freeze 8), broken down by singletons and variant functional category. M, million. c, Variant counts in different functional categories, proportion of singletons and per-individual median values. All counts were based on variants in the canonical transcript. pLOF includes frameshift, essential splice donor and acceptor (excluding splice sites in UTRs) and stop-gained variants.

Full size image

We performed a comprehensive survey of genetic variation, encompassing single-nucleotide variants (SNV) and insertion–deletion (indel) variants. To estimate population allele frequencies, we focused on 821,979 unrelated samples (referred to hereafter as the 822K unrelated set; Supplementary Table 1a). We identified 16,425,629 unique mutated genomic positions (that is, sites) in autosomal and X-chromosomal coding regions, with one unique reference–alternate allele change (that is, variant) every two bases on average. In canonical transcripts within sequencing target regions, mutations at 35.6%, 32.2% and 9.5% of all possible genomic positions that can lead to synonymous, missense and stop-gained variants, respectively, were observed. In highly methylated CpG sites, we observed 95.0% of all possible synonymous, 92.2% of missense and 78.6% of stop-gained variants. Across all mutational contexts, 21.4% and 8.4% of all possible synonymous variants and stop-gained variants, respectively, were observed (Extended Data Fig. 1). Thus, RGC-ME represents a major advance towards the comprehensive discovery of rare variants.

Among coding variation in canonical transcripts, 1,115,116 pLOF variants were identified, which include those causing a premature stop, affecting essential splice donor and acceptor sites or causing frameshifts (Fig. 1c). Of these pLOF variants, 53.3% were observed as singletons; that is, only observed in one individual. In addition, 4,645,092 synonymous (35.7% as singletons) and 10,444,562 missense (40.0% as singletons) variants in canonical transcripts were detected. A total of 48% of coding variants in canonical transcripts were unique to RGC-ME and absent in other large-scale datasets4,8 (Fig. 1b). Each sample had a median of 137 pLOF, 8,652 missense and 10,184 synonymous variants (Fig. 1c). AFR individuals had, on average, 18.6% more variants across all functional categories compared with individuals of other ancestries (Extended Data Fig. 2), as expected on the basis of the ‘Out of Africa’ model of human population history26.

Constrained genes

Population-scale sequencing allows the quantification of pLOF variation in genes, which is key to understanding the relationship between genes and diseases. Several gene constraint metrics have been developed to estimate the pLOF tolerance of genes27. Here, we estimated pLOF depletion using the cumulative frequency of pLOF variants in a gene to derive a selection coefficient, shet, that quantifies fitness loss due to heterozygous pLOF variation28. We estimated the indispensability of 16,710 protein-coding genes on the basis of the observed number of rare pLOF variants per gene with a cumulative alternate allele frequency (AAF) of less than 0.1% compared with the expected number based on gene-specific mutation rates (Supplementary Table 2).

The mean shet value in the RGC-ME dataset for canonical transcripts was 0.073 (95% highest posterior density (HPD)): [0.043, 0.12] (median shet = 0.021) (Fig. 2a), which suggests that, on average, a pLOF would result in 7.3% lower evolutionary fitness relative to the reference allele. This estimate is comparable with the mean shet value of 0.073 [0.029, 0.18] (median = 0.028) computed using the same method on the ExAC dataset29 (n ≈ 60,000). Our sample size (n ≈ 822,000) helped to accurately quantify rare pLOF variants and compute more precise constraint scores than the ExAC values. This finding is best illustrated in known haploinsufficient genes, which are expected to be more constrained and thus have larger shet values relative to all genes (Extended Data Fig. 3a). Compared with values that were computed with ExAC data, shet values for haploinsufficient genes in the RGC-ME dataset were significantly higher (\({\Delta \bar{s}}_{\text{het}}=0.045\), P = 0.002) and had smaller 95% HPD ranges despite those larger means (∆Var(shet) = −0.026, P = 4.5 × 10−21). Estimates for all genes were more precise in 822K samples compared with a randomly downsampled set of 60,000 samples from the RGC-ME dataset (Extended Data Fig. 3b), in which mean and median 95% HPD ranges were 6.2- and 4.0-fold larger, respectively.

Fig. 2: Estimates of gene-level constraint, representing shet, from the RGC-ME dataset.
figure 2

a, Mean shet probability density for 16,710 canonical transcripts with 95% CIs calculated with 10,000 bootstrapped samples from the means of individual genes. b, Odds ratios (points) and 95% CIs (short horizontal lines; computed using standard error) for genes with shet cut-off > 0.073 (deemed highly constrained genes) to be included in each gene category listed on the y axis compared with genes below the cut-off. Genes defined as ‘human knockouts’ are those with carriers of rare, biallelic pLOF variants observed in the RGC-ME dataset. A total of 16,710 canonical transcripts were included in each category, which contained at minimum 234 ‘true’ genes (that category being haploinsufficient genes). HGMD, Human Gene Mutation Database.

Full size image

The shet value is higher in genes that are associated with Mendelian diseases28,30 (Extended Data Fig. 3a), and can differentiate groups of genes under varying degrees of selection (Fig. 2b). We used shet to identify constrained genes by comparing the shet scores of known high-constraint genes (haploinsufficient, autosomal dominant and developmental-specific autosomal dominant) with those of low-constraint genes (haplosufficient and genes with rare biallelic pLOF variants from the RGC-ME dataset) (Extended Data Fig. 4a). Among 1,476 genes in the ‘high-constraint’ and 3,893 genes in the ‘low-constraint’ groups, 89.1% of genes with a shet score greater than the mean (0.073) and 66.6% of genes with a shet score greater than the median (0.021) belonged to the high group (Supplementary Table 2). These thresholds served as cut-offs for mean and lower bound (2.5% HPD), respectively, to identify highly constrained genes with fitness deficits on a par with dominant disease-causing genes that also reflect uncertainty in the mean.

We compared shet to other published LOF constraint measures, such as LOEUF4, and an alternate method for estimating shet based on approximate Bayesian computing31, which we refer to as shet-ABC (Supplementary Figs. 2–4). Spearman rank correlations between shet from RGC-ME and these estimates were high (−0.768 with LOEUF; 0.778 with shet-ABC). However, shet derived from RGC-ME had higher sensitivity and specificity in differentiating between constrained and unconstrained genes, compared with LOEUF and shet-ABC (Extended Data Fig. 4a).

Improving shet estimates is most valuable for genes with few expected pLOF mutations4, particularly shorter genes. The RGC-ME dataset allowed shet to be estimated more precisely for the smallest quantiles of gene coding sequence (CDS) length (Extended Data Fig. 3c), and derived more informative constraint metrics using an allele-frequency-based approach and a larger sample. We derived constraint scores for 923 genes that had 5 or fewer expected pLOF variants, deemed underpowered for similar analyses with LOEUF32. These 923 underpowered genes were significantly shorter, with a mean CDS length of 573 base pairs compared to 1,797 base pairs for genes with more than 5 expected pLOF variants. Eighty-six genes were highly constrained, with a mean shet value greater than 0.073 and a lower bound greater than 0.021 (Extended Data Fig. 4b), and are promising candidates for efforts to discover new disease-associated genes. Thirty per cent (26 of 86) have been linked to human diseases or shown to be essential in mice or cell lines (Supplementary Table 2). These include well-studied genes with known importance in cellular function, such as the transcription factor TWIST1 (ref. 33), DNA- and RNA-binding protein BANF1 (refs. 34,35) and transactivator CITED2 (ref. 36).

Overall, 3,988 highly constrained genes had shet values greater than 0.073 and a lower bound greater than 0.021. Although 1,153 of these lack known associations with human diseases or lethal mouse knockout phenotypes, they are likely to have high functional importance. These constrained genes might lack disease associations because the loss of even a single copy is incompatible with life or causes reduced reproductive success without clinical disease37.

Constrained coding regions

Identifying sub-genic regions that are intolerant of mutations can reveal functionally important regions that would otherwise be missed when constraint scores are aggregated at the gene level. Models of local coding constraint are powerful tools for identifying protein domains with crucial functions and for variant prioritization38,39,40. In addition to gene constraint derived from pLOF variation, we also identified regions depleted of missense variation using the missense tolerance ratio (MTR)40,41, defined as the ratio of the observed to the expected proportion of missense variants adjusted by synonymous variation in a defined codon window. Using the 822K unrelated samples, we calculated the MTR for each amino acid along the CDS within sliding windows of 21 and 31 amino acids (MTR scores available on figshare; see ‘Data availability’) and characterized continuous segments of missense constrained regions (Supplementary Table 3).

Compared with benign missense variants, ClinVar pathogenic missense (two stars or more) variants were highly enriched in the top percentile of exome-wide MTR scores (odds ratio = 100.0 and 89.8, computed with 21- and 31-codon windows, respectively; Fig. 3a). Our sample size, which is nearly four times larger than that used in previous MTR estimates41, resulted in improved discrimination between pathogenic and benign variants for top-10-percentile MTR scores in which we observed significant enrichment (Fig. 3a). This larger sample size enabled us to identify 24% more missense variants in the top-1-percentile constrained MTR scores (512,499 versus 413,147) compared with a subsampled set of 225,000, after adjusting for a false discovery rate (FDR) lower than 0.1. In addition, the increased power derived from 822K samples resulted in higher resolution for distinguishing pathogenic from benign variants for MTR computed with 21-codon windows, albeit at the expense of having fewer scored missense variants overall (295,958 constrained missense variants).

Fig. 3: Missense regional constraint captured by MTR.
figure 3

a, Odds ratio (OR) (points) and 95% CIs (error bars, two-sided Fisher’s exact test) of ClinVar pathogenic versus benign variants in MTR ranking regions across the whole exome. Comparisons include MTR calculated using the 822K unrelated samples from the RGC-ME dataset on 31-codon (blue) and 21-codon (pink) windows, and MTR calculated using a random subset of 225,000 samples from the larger 822K samples using a 31-amino-acid sliding window (yellow). MTR values include variants for which the FDR is less than 0.1. A total of 21,047 benign ClinVar variants and 12,872 pathogenic ClinVar variants (two stars or more (CV2+)) were included. b, MTR ranking distribution of different protein functional regions and variant groups (in order: DNA-binding sites (n = 2,787), ClinVar pathogenic sites (n = 10,673), active sites (n = 2,787), transmembrane region (n = 2,787), localized to extracellular (n = 25,665), localized to cytoplasm (n = 32,994) and ClinVar benign sites (n = 20,739)). Box plot shows median and 25–75% interquartile range. The whisker minima and maxima represent the smallest and largest data points within 1.5× the interquartile range from the lower quartile and upper quartile, respectively. Every functional category was significantly more constrained than the category to its right with two-sided Wilcoxon rank-sum test (least significant p-value = 2 × 10−4, Bonferroni correction). c, Distribution of the proportion of genes located in exome-wide top-15-percentile MTR regions against the heterozygous selection coefficient, shet. Genes with a significant proportion in the most constrained 15-percentile MTR region are coloured in pink and yellow (P < 0.05, Bonferroni-corrected, one-sided binomial tests with π0 = 0.15), stratified by LOF constraint (shet = 0.073). Pink dots highlight genes that are tolerant of LOF, but have some regions depleted of missense variation. Blue lines indicate the density of genes with shet scores from 0 to 1 (right margin) and genes with a proportion of MTR in the top 15 percentile exome wide (above plot). d, MTR track of an oncogene, KRAS, a missense-specific constrained gene, along with the domain structure of the protein. The blue MTR-constrained region is defined by top-15-percentile exome-wide MTR rank. The N-terminal region containing amino acids 1–80 is depleted of missense variation, even though KRAS is tolerant of heterozygous LOF variation (shet = 0.002).

Full size image

Deleterious variants are expected to have lower allele frequencies than neutral variants, owing to negative selection. We can infer the functional importance of different classes of variation by comparing the proportion of singletons in each class. We computed the deleteriousness of variants using an updated mutability-adjusted proportion of singletons (MAPS) metric5,32 and derived an MTR score threshold at which their MAPS score corresponds to that of missense variants that were predicted to be deleterious by five out of five prediction algorithms in dbNSFP (v.3.2; see Supplementary Information); that is, 5/5 missense variants. Variants with MTR values in the top-15-percentile exome-wide threshold (MTR < 0.841) were predicted to be as deleterious as 5/5 missense variants (Extended Data Fig. 5a). For 31-codon windows, 1.24% (129,990) of all missense variants (excluding known ClinVar pathogenic variants) observed in the RGC-ME dataset had significant MTR scores in the top 15 percentile. These missense variants in the top 15 percentile of exome-wide MTR are potentially deleterious and could be suitable for prioritization in projects aiming to discover disease-associated genes.

MTR is a useful metric of regional constraint that may capture functionally important segments within genes. We defined MTR-constrained regions as continuous regions within a protein that have variants with MTR values in the top-15-percentile threshold (Supplementary Information and Extended Data Fig. 6a). We identified 41,114 missense constrained regions in 12,349 genes (Supplementary Table 3). Our findings overlap with results from a previous study38 that estimated the regional observed-to-expected missense ratio (ɣ) from around 60,000 ExAC samples (Extended Data Fig. 6b) to derive a composite missense deleterious score called MPC. We refer to the MPC-derived constrained regions as MPC segments, and compared these with MTR-constrained regions. MTR-constrained regions had a median length of 22 residues [14–35, quartile 1–quartile 3], compared with 358 [208, 579] in MPC segments. Overall, we identified 8.59 times more MTR-constrained regions than MPC segments (ɣ ≤ 0.612, top 15 percentile) across 2,832 transcripts with data from both methods (Extended Data Fig. 6c).

We examined the distribution of de novo missense variants in MTR-constrained regions and observed a significant enrichment (P = 2.61 × 10−10) of variants identified in individuals with neurodevelopmental disorders (Extended Data Fig. 7a,b). Case variants were 1.85 times [1.50, 2.31 (95% CI)] more likely to occur in constrained regions, compared with controls. As expected, well-supported (two stars or more) ClinVar pathogenic missense variants were also highly enriched (P ≈ 0) in MTR-constrained regions. Pathogenic variants were 8.82 times [8.17, 9.53 (95% CI)] more likely to occur in missense constrained regions than were benign variants.

Missense constrained sites were found in key functional regions, such as DNA-binding regions and active sites (Fig. 3b). Among membrane proteins, transmembrane regions ranked higher in MTR-constrained regions than did cytoplasmic and extracellular domains. We also compared the overlap of MTR-constrained and functional regions by computing Jaccard indices. Ubiquitin-conjugating (UBC) core domains and DNA-binding regions had the highest overlap with constrained regions (Jaccard index = 0.52 and 0.18, respectively), suggesting that, among UBC enzymes, more than half of the union set between MTR-constrained regions and core domains overlapped. Other enriched functional regions included protein kinases and nuclear receptor ligand-binding domains (Supplementary Table 4).

A total of 4,064 genes contained regions depleted in missense variation with a significant proportion of their coding sequence in the top 15 percentile of MTR (binomial test with π0 = 0.15, P < 0.05 after multiple testing correction; Supplementary Table 5). To identify genes with signatures of missense-only constraint, we assessed the LOF-constraint metric, shet, of these highly missense constrained genes (Fig. 3c). Among the 4,064 genes, 1,482 either were not LOF constrained or lacked shet estimates. These genes had significantly shorter CDS lengths than those of the 1,424 LOF-specific constrained genes (P = 2.9 × 10−40, Wilcoxon test; Extended Data Fig. 7c). Estimating region-level LOF constraint is difficult owing to strong selection against pLOF variants, which leads to a paucity of pLOF variation. MTR serves as a complementary lens for identifying, first, functionally important regions at a higher resolution than gene-level LOF constraint, and, second, regions within genes that are depleted of missense variation but tolerant of LOF variation. For example, KRAS, a well-known oncogene, is LOF tolerant (shet = 0.002, LOEUF = 1.24); however, the first 80 amino acids (42%) of the protein sequence were ranked in the top 1 percentile of exome-wide MTR (Fig. 3d). This region includes the P-loop, switch 1 and switch 2 functional domains, which form crucial binding interfaces for effector proteins42, and these results therefore highlight the importance of regional constraint metrics.

Understanding ‘human knockouts’

Identifying genes with biallelic pLOF variants provides an opportunity to understand gene function directly through the phenotypic characterization of individuals who have such variants—effectively, naturally occurring ‘human knockouts’. The RGC-ME dataset includes founder populations and cohorts with high rates of consanguinity, contributing to a comprehensive collection of homozygous loss-of-function variation25,43,44,45. Overall, we identified 4,686 genes comprising 8,576 rare (AAF < 1%) homozygous pLOF variants in 64,852 individuals (Supplementary Table 6). Furthermore, we identified 1,205 genes with carriers of rare (AAF < 1%) heterozygous pLOF variants in trans; that is, compound heterozygotes, 162 of which lacked homozygous pLOFs. In total, 4,848 genes were discovered with carriers of biallelic pLOF variants in which both alleles of a gene were affected by pLOF variation and could be described as putative gene knockouts (pKOs). Of these, 1,751 (1,650 from homozygous pLOFs only) have not to our knowledge been previously reported. Biallelic pLOF variants in RGC-ME are rare; 64.3% of homozygous pLOF variants and 37.4% of pKOs were detected in one participant (Fig. 4a). As expected, cohorts with higher rates of consanguinity were enriched in homozygous pLOF variants, compared with outbred populations, despite smaller sample sizes (Fig. 4b,c and Extended Data Table 1).

Fig. 4: Rare biallelic pLOF variants and ‘human knockouts’ in the RGC-ME dataset.
figure 4

a, Distribution of the number of individuals per pKO on the log10 scale. Carriers of homozygous pLOFs and compound heterozygous variants were included in this analysis. b, Breakdown of the number of unique pKO genes observed in the RGC-ME dataset by ancestry. Both sets of rare biallelic variants—homozygous pLOFs and compound heterozygous—were included in this analysis. See Extended Data Table 1a for a breakdown by ancestry of each type. c, Projected accrual of pKO genes using homozygous pLOF variant data at hypothetical cohort sizes for each ancestry in 983,578 related individuals. Curves reflect the accrual of the expected number of genes with at least one, at least five and at least ten carriers, respectively, of a rare, homozygous pLOF. Asterisk denotes the inclusion of cohorts with a high rate of consanguinity.

Full size image

pKOs were significantly less constrained, with a lower shet (on average −0.074 [−0.077, −0.071 (95% CI)], t-test) relative to all other genes. Only 2.67% of pKOs had an shet value greater than 0.073, as compared with 21.6% of all human genes, and 47.2% of pKOs were in the lowest quintile of shet scores exome-wide (shet < 7.07 × 10−3). A caveat is that shet, like most gene-specific measures of constraint, is designed to capture the effect of heterozygous LOF46. Although genes containing biallelic pLOF variants are under less heterozygous selective pressure, existing sample sizes are inadequate47 to directly compute selection on homozygous variation. pKOs are overrepresented in drug and xenobiotic metabolism pathways (Supplementary Fig. 6).

Among very rare doubleton variants for which we observed exactly two copies of the alternate allele, we observed a clear excess of homozygotes that is likely to be explained by population structure and background inbreeding. For example, among missense and synonymous variants, we observed 5,857 and 2,490 homozygotes among 1,580,917 and 679,335 doubleton variants, respectively, compared with a Hardy–Weinberg equilibrium (HWE) expectation of fewer than one homozygote in each case. These estimates corresponded to a background inbreeding coefficient of 0.37%. Among pLOF variants, we observed only 406 homozygotes among 129,405 doubleton variants (Supplementary Table 7). Although this number is much larger than HWE expectations, it is around 15% less than the expected 479 homozygotes calculated using an inbreeding coefficient of 0.37% (P = 0.0095, Fisher’s exact test). This suggests that a notable proportion of these homozygotes were never observed in our sample population.

Genes with biallelic inactivating mutations could reveal potential drug targets that can be disrupted with minimal side effects43. Drug targets with homozygous pLOF variants in humans are more likely to progress from phase I trials to approval44. Of 997 inhibitory preclinical targets listed in the Drug Repurposing Hub, 182 (18.3%) had at least one individual with a rare biallelic pLOF variant in the RGC-ME dataset48. In-depth phenotyping of human knockouts can help researchers to better understand the efficacy and side-effect profiles of these potential drug targets. Human knockouts provide a way to understand the consequences of lifelong deficiency of a gene49.

Annotation of splice-affecting variants

Several prediction tools50,51,52,53 have been developed to understand the effects of genetic variants on alternative splicing. Although these tools mainly assess whether a variant affects splicing, some also provide a pathogenicity metric or score threshold as a measure of deleteriousness. Predicted cryptic splice sites with SpliceAI scores greater than 0.8 have been validated at high rates using RNA sequencing and are as depleted at common allele frequencies as pLOF variants50. Here, we used human genetic data to optimize splice prediction score thresholds enriched for deleterious variants that affect splicing. We systematically quantified the deleteriousness of variants at various splice prediction score thresholds using the MAPS metric. As previously demonstrated2,4, pLOF variants had the highest MAPS scores, followed by missense, synonymous and noncoding variants, respectively (Fig. 5a).

Fig. 5: Identification of deleterious variants that are predicted to affect splicing.
figure 5

a, MAPS across different functional categories. Error bars show standard deviation around the mean proportion of singletons (points). The yellow dashed line represents the SpliceAI and MMSplice score threshold for variants that have a MAPS score equal to that of 5/5 missense variants (predicted deleterious by five algorithms). Variants with a SpliceAI score ≥ 0.35 or a MMSplice score ≥ 0.97 are predicted deleterious SAVs. Noncoding variants refers to intronic, downstream (variant located 5′ of a gene), upstream (variant located 3′ of a gene) and 5′ and 3′ UTR variants captured by exome sequencing. Coding variants are inclusive of canonical splice sites, splice region and UTR splice sites. All variants that passed quality control and were observed in unrelated individuals in the RGC-ME dataset were included in this analysis (n = 34,512,842 variants). b, Enrichment of ClinVar pathogenic variants (two stars or more) in predicted SAVs compared with corresponding variant sets filtered by either LOFTEE, 5/5 missense deleteriousness models or CADD. Points represent odds ratios and bars depict 95% CIs (two-sided Fisher’s exact test, no multiple testing correction). ‘All variants’ include 313,390 coding and noncoding variants, and ‘splice sites’ include essential and UTR splice sites; counts of variants included in each calculation are provided in the table. HC, high-confidence. c, Empirical validation of MAPS-predicted deleterious SAVs (intersection set): enrichment of predicted deleterious SAVs in experimentally validated SDVs compared with non-SDVs. Points represent odds ratios and bars depict 95% CIs (two-sided Fisher’s exact test, no multiple testing correction). A total of n = 36,636 variants, of which 346 SAVs are validated SDVs, are included.

Full size image

We used splice predictions from SpliceAI50 and MMSplice51 to group variants into predicted splice score bins, and identified the minimum threshold at which the MAPS score of the variants is equal to that of 5/5 missense variants (variants predicted to be deleterious by five out of five prediction methods). The proposed prediction score thresholds of 0.35 for SpliceAI and 0.97 for MMSplice pathogenicity (Fig. 5a and Extended Data Fig. 5b) identify predicted deleterious splice-affecting variants (SAVs).

A total of 296,696 predicted deleterious coding SAVs (inclusive of canonical splice sites, splice region and untranslated region (UTR) splice sites) in the RGC-ME dataset had scores that exceeded the MAPS-derived splicing thresholds for both SpliceAI and MMSplice (referred to as the intersection set; Extended Data Fig. 8a). Of these, 43.5% (129,118) were cryptic splice sites (that is, non-canonical splice sites). Unsurprisingly, canonical splice sites and variants within the splice region comprised the largest category of predicted deleterious SAVs. Both SpliceAI and MMSplice identified around 80% of LOFTEE (loss of function transcript effect estimator; ref. 4) high-confidence splice sites and around 10% of variants within splice regions as predicted deleterious SAVs (Extended Data Fig. 8a,b). In addition, around 68% of LOFTEE low-confidence splice sites were predicted to be deleterious SAVs (94% of low-confidence splice sites were in the UTR). The impact of non-canonical splice variants on alternative splicing is often underestimated; we found that missense variants accounted for 11.3% of all predicted deleterious SAVs identified by both SpliceAI and MMSplice in the RGC-ME dataset (Extended Data Fig. 8a,b).

Predicted deleterious SAVs were enriched in well-supported ClinVar pathogenic variants (two stars or more) compared with other metrics of variant deleteriousness (Fig. 5b); for example, compared with combined annotation dependent depletion (CADD)54,55 score ≥ 20 (odds ratio = 4.5, P = 0). Missense SAVs were significantly enriched for pathogenic variants compared with 5/5 missense variants (odds ratio = 1.8, P = 3.3 × 10−8) and missense variants with CADD ≥ 20 (odds ratio = 3.6, P = 1.01 × 10−26), respectively. Notably, splice sites in the intersection set were also significantly enriched for pathogenic variants compared to LOFTEE high-confidence splice sites, indicating that the MAPS-derived metric identifies deleterious splice sites. Similar results were obtained when we evaluated the enrichment of pathogenic variants compared with benign ones (Supplementary Table 8a).

We next assessed the MAPS-derived splice prediction thresholds for variants that have been experimentally assessed for splicing effects56,57,58 (Supplementary Table 9). Predicted deleterious SAVs identified in the intersection set were significantly enriched in experimentally validated large-effect splice-disrupting variants (SDVs) compared with non-SDVs in all functional categories except the splice site category, although the odds ratio was greater than one for splice sites (Fig. 5c). Variants of unknown significance (VUSs) in ClinVar that were predicted as deleterious SAVs were also significantly enriched in experimentally validated SDVs (Fig. 5c). Of the 563 predicted deleterious SAVs assayed in the experimental data, 346 (61.5%) were SDVs and more than half were cryptic splice sites, including 13 ClinVar VUSs (Extended Data Fig. 8d).

We also derived stringent thresholds to identify SAVs by removing canonical splice sites and calibrating exclusively coding non-splice-site (nonSS) variants to a MAPS score comparable with 5/5 missense variants. These thresholds corresponded to a SpliceAI score of 0.43 and an MMSplice score of 0.97 (Extended Data Fig. 5c). Pathogenic enrichment was consistent when comparing deleterious coding nonSS and missense SAVs with corresponding variant categories filtered by CADD ≥ 20 (Supplementary Table 8b,c). Consistent results were also obtained when comparing the enrichment of deleterious SAVs in SDVs to non-SDVs after applying thresholds for coding nonSS variants (Extended Data Fig. 8c,e).

Clinical utility of rare variants

To understand the prevalence of disease-associated alleles in the general population, we identified well-supported ClinVar59 pathogenic variants (two stars or more) across 2,042 genes in 822K unrelated RGC-ME samples. We found that 40.7% of pathogenic variants (20,343/49,990) were observed in the RGC-ME dataset, of which 99.6% (20,262) had an AAF of less than 0.1% and 17.8% (3,619) were observed once. In comparison, 20% (9,821) and 29% (14,700) of pathogenic variants were observed in ExAC exomes (n ≈ 60,000) and gnomAD v.2.1.1 exomes (n ≈ 126,000), respectively (Extended Data Fig. 9a). This highlights the importance of the RGC-ME dataset’s larger sample size in identifying rare pathogenic variants. On average, individuals carry 1.58 pathogenic variants, with the majority of these individuals being heterozygous carriers of these variants. Specifically, 61.4% of the 822K unrelated individuals were heterozygous carriers of pathogenic recessive alleles in 1,143 of 2,659 known autosomal recessive genes (mean, 0.98 pathogenic alleles per person); 0.21% of the samples were homozygotes of pathogenic variants in 167 autosomal recessive genes; and 3.64% were heterozygous carriers of 353 of 1,629 total autosomal dominant genes. Pathogenic variant annotations should be interpreted cautiously owing to the incomplete penetrance of disease alleles60.

The American College of Medical Genetics identified a set of genes (ACMG SF v.3.1) with clinically actionable variants that predispose individuals to diseases and for which medical interventions are available to reduce mortality and morbidity61. Among the 822K unrelated individuals, 22,846 (2.77%) had at least one ClinVar-reported (two stars or more) pathogenic missense or pLOF variant for 72 out of 76 autosomal genes on the ACMG list (Supplementary Table 10). As expected, two of the most prevalent pathogenic variations were the HFE Cys282Tyr allele (enriched in EUR, nEUR-homozygotes = 3,220 and AAFEUR = 13.8%) and the TTR Val142Ile allele (enriched in AFR, nAFR = 1,670 and AAFAFR = 3.4%).

We also tallied carriers of likely pathogenic pLOF variants (novel variants not yet reported as pathogenic in ClinVar) in 44 genes in which truncation is known to lead to disease. A total of 2,357 (0.3%) individuals in the RGC-ME dataset carried 1,407 likely pathogenic variants across 40 of these genes. In total, 3.06% of the individuals in the RGC-ME dataset were carriers of pathogenic or likely pathogenic variants. Excluding individuals with high-frequency pathogenic variants in the HFE (Cys282Tyr) and TTR (Val142Ile) genes, 2.38% of the individuals in the RGC-ME dataset carried an actionable variant (Supplementary Table 10). This number is comparable with those from other reports6,7,62 of actionable variants, which range from 2% to 4.1% for gene sets that include ACMG v.2.0 and v.3.0. As expected, pathogenic variants are rare in large-scale studies of the general population. We found that 39% and 79% of pathogenic and likely pathogenic variants, respectively, were singletons. Focusing on non-ACMG genes, we found that 1.27% of individuals were heterozygous carriers of pathogenic variants in autosomal dominant genes, and 0.21% were homozygotes of pathogenic variants in autosomal recessive genes.

Because the RGC-ME dataset includes uniformly processed exome data from a relatively large proportion of individuals from continental ancestries other than EUR, we assessed the range of allele frequencies of variants present in ClinVar across four continental populations: AFR, EUR, IAM and SAS. Approximately 34% of unique pathogenic coding variants in equalized subsamples were observed only in individuals of non-EUR ancestry, which indicates that sampling diverse populations is necessary for the comprehensive identification of rare variation. Across all unrelated individuals, on average, those of EUR ancestry had 63% more pathogenic variants that were well characterized (rated two stars or more) per sample than did individuals of AFR ancestry. Conversely, individuals of EUR ancestry had, per sample, 25.6% fewer VUSs (Extended Data Fig. 9b,c) and 18.6% fewer variants across all functional types (Extended Data Fig. 2). In individuals of AFR ancestry, a consistent pattern of significantly fewer high-confidence (two stars or more) pathogenic variants (−0.576 [−0.567, −0.585 (95% CI)], t-test) to a surplus of VUSs (42.13 [421.97, 42.28]), compared with individuals of EUR ancestry, suggests that the most well-characterized pathogenic variants were depleted in this population (Extended Data Fig. 9b). Recruiting diverse individuals to enable the identification and characterization of novel pathogenic variants might help to address this ascertainment bias. Further analyses of pathogenic coding variants and differentiated alleles between ancestries are included in Supplementary Fig. 7 and Supplementary Table 11.

Understanding VUSs is currently a bottleneck in the interpretation of variation in clinically relevant genes and a challenge in clinical management19. Although VUSs have less empirical evidence for pathogenicity, they comprise the bulk of ClinVar, with more than one million variants. Notably, VUSs in regions of low MTR may be deleterious, comprising 5,079 (0.68%) VUSs in the top 1 percentile of MTR-constrained regions and 17,500 VUSs (2%) in the top 15 percentile (Supplementary Table 12). Using the MAPS-derived splicing score thresholds, we identified more than 11,000 candidate deleterious cryptic splice sites among VUSs (1,366 synonymous variants in 822 genes and 10,407 missense variants in 3,501 genes), offering potential insights into their functional consequences for clinical prioritization and interpretation efforts.

Discussion

The RGC-ME dataset, derived from 983,578 exomes, provides a harmonized catalogue of around 20 million coding variants in individuals from a diverse array of ancestries and is publicly accessible at https://rgc-research.regeneron.com/me/home. Cataloguing variation at scale provides an opportunity to accurately estimate the frequency of rare variants—allowing us to precisely compute gene and regional constraint metrics, expand the compendium of rare human knockouts, annotate deleterious cryptic splice sites, characterize variant frequencies across different ancestries and assess the population prevalence of pathogenic variation. RGC-ME will be an invaluable resource for interpreting rare variants and is a step towards the realization of precision medicine.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.