Leveraging Base Pair Mammalian Constraint to Understand Genetic Variation and Human Disease

Abstract

Although thousands of genomic regions have been associated with heritable human diseases, attempts to elucidate biological mechanisms are impeded by a general inability to discern which genomic positions are functionally important. Evolutionary constraint is a powerful predictor of function that is agnostic to cell type or disease mechanism. Here, single base phyloP scores from the whole genome alignment of 240 placental mammals identified 3.5% of the human genome as significantly constrained, and likely functional. We compared these scores to large-scale genome annotation, genome-wide association studies (GWAS), copy number variation, clinical genetics findings, and cancer data sets. Evolutionarily constrained positions are enriched for variants explaining common disease heritability (more than any other functional annotation). Our results improve variant annotation but also highlight that the regulatory landscape of the human genome still needs to be further explored and linked to disease.

Fig. 1.

(A) Evolutionary constraint in multiple genomic partitions. X-axis=fraction of the genome occupied by a partition, Y-axis=fraction of partition under constraint in placental mammals (orange circles) and primates (blue triangles), grey line is the genome mean (0.035). The greatest constraint is found in CDS and key regulatory regions (5’UTR, ENCODE promoter-like elements, and 3’UTR). This figure is a subset of fig. S1 which shows more biotypes, protein-coding gene parts, and regulatory regions. (B) Whisker plots of constraint in variants from TOPMed WGS, stratified by CDS (red, 6.14 million biallelic SNPs) and non-CDS variants (orange, 549.64 million biallelic SNPs). X=six allele count (AC) bins, from singletons (AC=1, 44.8%) to common variants (allele frequency, AF ≥ 0.005, 1.4%). (C) PhyloP score density for ClinVar benign (N=231,642), ClinVar pathogenic (N=73,885), and gnomAD WGS variant positions with CADD ≥ 20 (N=3,958,488).

Fig. 2.. SNP-heritability analyses of variants at constrained positions in human complex traits and diseases.

(A) Heritability enrichment of common SNPs in the top percentiles of constraint scores in placental mammals (phyloP) and primates (phastCons). (B) Heritability enrichment as a function of the distance to a constrained base. (C) Heritability enrichment of constrained annotations in 11 blood and immune traits and 9 brain diseases (light color) versus other types of traits (dark color). Asterisks indicate significance at P < 0.05 and double Asterisks indicate significance at P < 0.05 after Bonferroni correction (0.05/4). (D) Heritability enrichment of constrained and functional annotations (left), and corresponding significance of the conditional effect while considered in a joint model with 106 annotations (right). (E) Heritability enrichment of constrained annotations intersected together and stratified by their genomic function. The dashed grey line represents heritability enrichment in coding regions (plotted for comparison purposes). (F) Squared trans-ancestry genetic correlation enrichment (left) with corresponding significance (right) for 7 annotations with significant depletion of squared trans-ancestry genetic correlations. (G) Standardized squared effect sizes as a function of allele frequency. Results are meta-analyzed across 63 independent GWAS (A, B, C, E), 31 independent traits with GWAS available in European and Japanese populations (F), and 27 independent UK Biobank traits (G). Dashed red lines represent a null enrichment of 1 (A-E) and a null squared trans-ancestry genetic correlation (f). Error bars are 95% confidence intervals. Numerical results are reported in tables S2, S3, S4, S6, S7, S8, and S11.

Fig. 3.. Leveraging constraint to move from prioritization to function.

(A,B) We report the cumulative distribution function (CDF) of posterior inclusion probability (PIP) scores using functionally-informed fine-mapping with different models of functional annotations. Distribution functions are split into subpanels by whether the fine-mapped SNP overlaps high constraint scores in mammals (A) and primates (B). One-way Komolgorov-Smirnov tests that CDF for PIP obtained from the baselineLF model (gray) are lower (above) than the CDF for PIP obtained from the baseline-LF+Zoonomia model (orange) with Bonferroni correction for N=4 categories across panels (*** p/N < 0.0001, N.S. not significant). (C,D) Examples of constrained fine-mapped variants. We report GWAS P-values (upper panel) and corresponding PIP under different functionally informed fine-mapping models (lower panel). Shape of the dots corresponds to constraint information. (E) Fine-mapped variants are not limited to the annotated genome as exemplified by rs72782676 (red dot in AF panel) in the GATA3 UNannotated Intergenic COnstraint RegioN (UNICORN) locus. (F,G) Constraint is formally linked to function via massively parallel reporter assays (MPRAs) at the (F) regional oligo and (G) base pair level for neutral, active and allele specific skewed effect. (H) For the LDLR promoter locus, MPRA effect is strongly correlated with phyloP score. Constrained (red), and unconstrained (orange) ClinVar pathogenic variants are plotted to highlight known deleterious positions.

Fig. 4.. Evolutionary constraint, protein-coding genes, and human disease.

(A) Scatterplot of protein-coding (PC) gene clustering (UMAP and DBSCAN). X- and Y-axis are the UMAP coordinates. Each point is a PC gene (N = 19,386). Five clusters are labeled: A = 56 genes whose CDS bases are in complex regions that align poorly; B = 221 genes apparently human- or primate-specific; C=669 genes with good alignment and possible human-specific functions (e.g., five HLA genes and 14 interferon alpha genes); D=15 genes, all highly constrained; and E = all other 18,425 PC genes. Coloring shows fracCdsCons, grey = least and red = most constrained with an anticlockwise gradient in mammalian constraint from upper middle to lower right. (B, C) Gene constraint deciles versus external gene sets as “lollipop plots”. Each panel has 6 subgraphs for autosomal recessive genes, ClinGen level 3 genes, essential genes from Hart, essential genes in mouse, olfactory receptors, and severe haploinsufficiency genes. X-axis = constraint decile (0 = least,, 9=most constrained, 99 = missing). Y-axis = circles are the fraction of the PC genes in a gene set in each decile. (B) Zoonomia fracCdsCons and (C) recapitulates Figure 3 from ref. (3) with LOEUF decile reversed and showing missing data. (D, H) Gene heritability enrichment for SNPs linked to genes of each decile of fracCdsCons (D) and of SNPs linked to genes of decile of constraint in different gene gene features (H). Dashed red lines represent a null enrichment of 1. Error bars are 95% confidence intervals. (E) Spearman correlation of constraint fraction between the parts of PC genes. (F, G) Fraction of CDS constraint (fracCdsCons) vs. fraction or promoter constraint (F) and fraction of distal enhancer constraint (shrinked to values <0.3) (G). Each point is a PC gene, and HOX genes (purple) and defensin beta (DEFB) genes (green) are highlighted.

Fig. 5.. Cancer driver genes identified using NCCM rates.

(A) Distribution of the rates of NCCM for medulloblastoma. (B) An example set of the candidate driver genes found either in pediatric (orange) or adult (purple) samples. Age of diagnosis (years) of the patient is indicated together with the tumor subgroup. (C) ZFHX4 locus contains 9 NCCMs drawn from 8 patients.