Rare coding variant analysis for human diseases across biobanks and ancestries

Abstract

Large-scale sequencing has enabled unparalleled opportunities to investigate the role of rare coding variation in human phenotypic variability. Here, we present a pan-ancestry analysis of sequencing data from three large biobanks, including the All of Us research program. Using mixed-effects models, we performed gene-based rare variant testing for 601 diseases across 748,879 individuals, including 155,236 with ancestry dissimilar to European. We identified 363 significant associations, which highlighted core genes for the human disease phenome and identified potential novel associations, including UBR3 for cardiometabolic disease and YLPM1 for psychiatric disease. Pan-ancestry burden testing represented an inclusive and useful approach for discovery in diverse datasets, although we also highlight the importance of ancestry-specific sensitivity analyses in this setting. Finally, we found that effect sizes for rare protein-disrupting variants were concordant between samples similar to European ancestry and other genetic ancestries (β_Deming = 0.7-1.0). Our results have implications for multi-ancestry and cross-biobank approaches in sequencing association studies for human disease.

Figure 5

Figure 1. Study overview for rare variant discovery across human disease.

Three studies were included in the analysis: All of Us (AoU) with whole-genome sequence data, UK Biobank (UKB) with exome sequencing data, and Mass General Brigham Biobank (MGB) with exome sequence data. Over 600 disease Phecodes were identified using a hierarchal clustering algorithm. Disease Phecodes were analyzed using exome-wide gene-based testing of rare genetic variants using three masks (LOF, LOF+missense and ultra-rare missense) after which P-values were combined into a single P-value using the Cauchy distribution for each gene-disease pair.

Figure 2. Multi-ancestry meta-analysis of rare genetic variation across three sequenced biobanks in over 750,000 individuals identifies 363 rare variant associations.

Panel a shows a stacked bar chart with the proportion of each continental ancestry on the y-axis and dataset on the x-axis. Ancestral diversity was largest in AoU. Panel b is a violin plot with overlaid boxplot showing the prevalence of Phecodes on the y-axis and each dataset on the x-axis. Plotted Phecodes were those included in the analysis with at least 50 cases in each dataset (N=546). Panel c is a stacked bar chart showing the number of identified disease associations (Cauchy Q<0.01) on the y-axis and each dataset and x-axis, as well as the meta-analysis results. Bars are stacked by the class of mask that yielded the lowest P-value (from LOF masks, LOF+missense masks, and ultra-rare missense masks). Panel d is a multi-trait gene-based Manhattan plot highlighting results from the overall meta-analysis, each dot representing one gene-trait test, with the -log10 of the Cauchy P-value on the y-axis and different disease categories on the x-axis. For disease categories with strong associations, the top three non-redundant genes are annotated with the gene names. Panel e is a violin plot with overlaid boxplot showing the distribution of inflation factors by phenotype (λ estimated at 95 percentile) on the y-axis, and different rare variant masks on x-axis, as well as the distribution for the Cauchy combination results (on the far right). Dotted lines show the 0.75 and 1.25 cutoffs for inflation factor on the y-axis. The number of phenotypes is 601 in all violins. Similarly, panel f shows the distribution of inflation factors by gene across the different masks and for the Cauchy combination results (on the far right), where the number of genes equals 14,388, 15,529, 17,809, 16,742, 15,462, 18,238, and 18,456. Note: Cauchy P-values represent the omnibus P-value of all masks for a gene-phecode pair (unadjusted for multiple testing) after combining them using the Cauchy distribution. The Cauchy Q-values represent the Benjamini-Hochberg FDR adjustments of these Cauchy P-values. P-values for mask-phecode pairs - prior to the Cauchy combination - were derived from Z-score-based meta-analyses of score tests from logistic mixed-effects models with saddle-point-approximation. All statistical tests and P-values are two-sided. All boxplots show median (center), 25^th percentile (bottom of box), 75^th percentile (top of box), smallest/largest value within 1.5*inter quartile range from hinge (bottom/top whiskers, respectively), and data points outside of this range (dots). Abbreviations: EUR, European ancestry; AFR, African ancestry; AMR, Admixed American Ancestry; SAS, South-Asian ancestry; EAS, East-Asian ancestry; UND, undefined ancestry; LOF, loss-of-function.

Figure 3. Assessment of bias from inclusion of non-European samples among the significant associations.

Panel a shows a scatter plot with each dot representing a gene-phecode pair that reached test-wide significance in our primary analysis (Q<0.01), with - log₁₀(P_Cauchy) from the primary analysis on the x-axis, and the -log₁₀(P_Cauchy) derived from a European ancestry sensitivity analysis on the y-axis (both log-transformed for clarity). Specific cutoffs on the y-axis are highlighted using dotted lines. Any strong deviation of P-values could indicate bias in our multi-ancestry approach, or alternatively indicate markedly lower power among European samples. No associations were abolished when restricting to European samples. There were 6 additional strongly attenuated genes (0.05>P_EUR>0.0005). Among these, several represent known gene-disease links (Supplementary Note). Panel b shows a scatter plot with the effect sizes for significant associations from the primary analysis on the x-axis, with the effect sizes from European-only sensitivity analyses on the y-axis. The effect size for the most significant mask is plotted for each gene-phecode pair, restricting to masks that had adequate allele counts in both the primary analysis and in the sensitivity analysis (cMAC≥20). Any large deviations from the dotted line (x=y) indicate bias from our multi-ancestry approach. Strikingly, no strong deviations of effect sizes were observed in this sensitivity analysis. For 8 associations there were insufficient alleles among European ancestry samples to compute an effect size, although represented well known gene-disease links (Supplementary Note). Taken together, these results show that the bias from inclusion of non-European samples was not substantial. Note: Bias is defined here as the spurious change in effect sizes / test statistics that is caused by inclusion of multiple ancestries but is not caused by true biological differences. Cauchy P-values represent the omnibus P-value of all masks for a gene-phecode pair (unadjusted for multiple testing) after combining them using the Cauchy distribution. The Cauchy Q-values represent the Benjamini-Hochberg FDR adjustments of these P-values. P-values for mask-phecode pairs (prior to the Cauchy combination) were derived from Z-score-based meta-analyses of score tests from logistic mixed-effects models with saddle-point-approximation. All statistical tests and P-values are two-sided. ORs were estimated using inverse-variance-weighted meta-analysis of two-sided Firth’s logistic regression results. Abbreviations: ALL, all ancestry individuals; EUR, European ancestry individuals; OR, odds ratio.

Figure 4. Large genetic effect sizes and pleiotropic associations identify core genes for the human disease phenome.

Panel a shows stacked bar charts for all genes from the meta-analysis that showed at least 3 associations, with the number of associations on the y-axis and gene on the x-axis. Bars are stacked by the class of the best mask (LOF, LOF+missense masks, or ultra-rare missense masks) for each gene-trait association. Panel b represents grouped boxplots showing rare variant effect size distributions per Phecode category, with log-scaled ORs on the y-axis and categories on the x-axis. The figure is restricted to gene-trait associations reaching Cauchy Q<0.01 and restricting to rare variant masks with P<2.6x10^-6. Per category, only masks with at least 7 associations are shown (and therefore some categories do not show all masks and not all categories are plotted). In box plots, the number of contributing associations from left to right equals 53, 44, 23, 15, 9, 16, 37, 30, 9, 15, 8, 22, 32, 54, 58, 18, 9, 9, and 8. All boxplots show median (center), 25^th percentile (bottom of box), 75^th percentile (top of box), smallest/largest value within 1.5*inter quartile range from hinge (bottom/top whiskers, respectively), and data points outside of this range (dots). Panel c represents a multiple jittered lollipop chart showing rare variant effect sizes for each Phecode category. The x-axis shows the log-scaled OR with each dot representing an association (restricting to gene-trait associations with Cauchy Q<0.01 and rare variant masks with P<2.6x10^-6), and Phecode categories on the y-axis. Horizontal lines start at 1 and end at the largest estimated effect size within the category. Dots are colored by class of rare variant mask. Select genes are annotated within each category to highlight large-effect size genes for the respective category. Note: In all panels, ORs were estimated using inverse-variance-weighted meta-analysis of two-sided Firth’s logistic regression results, while mask-phecode P-values were estimated from Z-score-based meta-analysis of score tests from logistic mixed-effects models with saddle-point-approximation. Cauchy P-values represent the omnibus P-value of all masks for a gene-phecode pair after combining them using the Cauchy distribution (unadjusted for multiple testing), while the Cauchy Q-values represent the Benjamini-Hochberg FDR adjustments of these P-values. All statistical tests and P-values are two-sided. Abbreviations: LOF, loss-of-function; OR, odds ratio.

Figure 5. Effect sizes of rare coding variants for disease correlate between genetic European and other genetic ancestries.

Figure 5 shows scatter plots with the effect sizes from European-ancestry analyses on the x-axis with the respective effect sizes estimated among individuals dissimilar to European ancestry on the y-axis. In each figure panel, a three-sample design was applied: Significant mask-disease pairs were identified from a EUR meta-analysis of UKB and MGB (significance determined at P<2.6x10-6), after which those mask-disease pairs were assessed within different ancestry groups from the AoU dataset. Each panel shows effect sizes (ie, log[OR]) for EUR analysis on the x-axis and effect sizes from other ancestries on the y-axis; the left panels show EUR versus all non-European samples, the middle panels show EUR vs African ancestry samples, and the right panels show EUR versus Admixed-American samples. Part a shows results for rare LOF variant masks with at least 20 carriers in both ancestry assignments, while part b shows results for ultra-rare missense0.5 variant masks with at least 20 carriers in both ancestry assignments. Linear trend lines from error-in-variable total-least-squares Deming regression are added to the plots. Statistics from Deming regression, including estimated β [95%CI] and P-values, are added in text in the top left corners. A regression coefficient (β_sens) and 95%CI is also provided in the bottom right corners, showing results from a combined sensitivity analysis where genes associated with age or leukemic outcomes are removed, and where analyses are adjusted for quantiles of effective sample size (Supplementary Table 14). Note: All ORs were estimated using Firth’s logistic regression models among unrelated participants. Deming regression was run using beta coefficients and their standard errors, making the analysis comparable to York regression with the assumption of uncorrelated errors. Standard errors were computed using Jackknife estimators. All statistical tests and P-values are two-sided. Abbreviations: LOF, loss-of-function; OR, odds ratio; EUR, European ancestry; AFR, African ancestry; AMR, Admixed American ancestry.; nonEUR, defined ancestry other than European; sens, sensitivity analysis.