Estimating the selective effects of heterozygous protein-truncating variants from human exome data

Abstract

The evolutionary cost of gene loss is a central question in genetics and has been investigated in model organisms and human cell lines. In humans, tolerance of the loss of one or both functional copies of a gene is related to the gene's causal role in disease. However, estimates of the selection and dominance coefficients in humans have been elusive. Here we analyze exome sequence data from 60,706 individuals to make genome-wide estimates of selection against heterozygous loss of gene function. Using this distribution of selection coefficients for heterozygous protein-truncating variants (PTVs), we provide corresponding Bayesian estimates for individual genes. We find that genes under the strongest selection are enriched in embryonic lethal mouse knockouts, Mendelian disease-associated genes, and regulators of transcription. Screening by essentiality, we find a large set of genes under strong selection that are likely to have crucial functions but have not yet been thoroughly characterized.

Figure 1

Inferred distribution of fitness effects for heterozygous loss of gene function. Estimates of parameters (α̂, β̂) from maximum likelihood fit to the observed distribution of PTV counts across 15,998 genes in terciles of mutation rate, assuming s_het ~ IG(α, β). Shaded areas show 95% CI obtained from 100 bootstrapping replicates, intended to quantify the influence of sampling noise in the data set on parameter inference, with fixed estimates of local mutation rate.

Figure 2

Separation of disease genes and clinical cases by mode of inheritance. [a] The percentage of genes associated with exclusively autosomal dominant (AD, N=867) disorders versus autosomal recessive (AR, N=1,482) disorders as annotated by the Clinical Genomics Database (CGD) in each s_het bin. Logarithmic bins are ordered from greatest to smallest s_het values. [b] Overall, AD genes have significantly higher s_het values than AR genes [Mann-Whitney U p-value 3.14×10⁻⁶⁴]. [c] Similarly, in solved Mendelian clinical exome sequencing cases (Baylor), s_het values can help discriminate between AR and AD disease genes, as annotated by clinical geneticists. [d] A s_het value of 0.04 can be used as a simple classification threshold for AD genes with a PPV of 96%. [e] This finding is replicated in a separately ascertained sample from UCLA. Box plots range from 25^th–75^th percentile values and whiskers include 1.5 times the interquartile range.

Figure 3

Enrichments of s_het in known haploinsufficient disease genes of high confidence (ClinGen Dosage Sensitivity Project). In (N=127) autosomal genes, we annotate the s_het scores of genes associated with each disease category and classification. Higher s_het values are associated with [a] earlier age of onset (Mann-Whitney U p=1.46 ×10⁻²), [b] a larger fraction of de novo variants (p=8×10⁻⁵), [c] high or unspecified penetrance (p=1.79 ×10⁻²) and [d] increased phenotypic severity (p=4.87×10⁻³). Box plots range from 25^th–75^th percentile values and whiskers include 1.5 times the interquartile range. [e] Genes with the 10% highest s_het values are also similarly enriched with more severe clinical annotations.

Figure 4

Distribution of s_het values for phenotypes in known disease genes and clinical cases. We plot the distribution of selective effects for different disorder groups, providing information about the breadth and severity of selection associated with each group. [a] We include known Mendelian disease genes (Clinical Genomic Database) annotated as either Autosomal Recessive or Autosomal Dominant and [b] clinical exome sequencing cases. We contrast these with [c] all tolerated knockouts in a consanguineous cohort (PROMIS) and [d] the distribution of selective effects in all scored genes. Logarithmic bins are ordered from greatest to smallest s_het values.

Figure 5

High-throughput screens of gene essentiality in mice and cell assays, as a percentage of all genes in each s_het bin. [a] Proportion of orthologous mouse knockout genes by phenotype, from a neutrally-ascertained set of genes generated by the International Mouse Phenotyping Consortium (IMCP). Logarithmic bins are ordered from greatest to smallest s_het values. [b] ICMP mice are separated into viable (N=1,057), sub-viable (N=211) and lethal knockouts (N=477), and lethal knockouts have significantly higher s_het values than viable [Mann-Whitney U p-value 2.95×10⁻²⁸]. [c] Cell-essential genes as reported by Wang et al. from genome-wide KBM-7 tumor cell CRISPR assay (N=1,740) have significantly higher s_het values [p-value 5.13×10⁻¹⁶] [d] as do genes that were characterized as essential in a gene trap assay (N= 1,081) [p-value = 4.90×10⁻¹⁸]. In the CRISPR assay, all genes with adjusted p-values < 0.05 and negative assay scores are included, and genes with gene trap scores < 0.4 or lower are included. Box plots range from 25^th–75^th percentile values and whiskers include 1.5 times the interquartile range.

Figure 6

Protein pathways and protein-protein interactions, as a percentage of the associated developmental genes in each s_het bin. [a] In key developmental pathways in KEGG, we find that genes with higher s_het values are enriched in genes important to development. [b] We plot the distribution of the number of protein-protein interactions for each gene, as determined by a genome-wide mass spectrometry assay versus s_het value. [c] We find that s_het values are positively correlated with the number of observed interactors for each gene. Box plots range from 25^th–75^th percentile values and whiskers include 1.5 times the interquartile range.