Powerful, scalable and resource-efficient meta-analysis of rare variant associations in large whole genome sequencing studies

Abstract

Meta-analysis of whole genome sequencing/whole exome sequencing (WGS/WES) studies provides an attractive solution to the problem of collecting large sample sizes for discovering rare variants associated with complex phenotypes. Existing rare variant meta-analysis approaches are not scalable to biobank-scale WGS data. Here we present MetaSTAAR, a powerful and resource-efficient rare variant meta-analysis framework for large-scale WGS/WES studies. MetaSTAAR accounts for relatedness and population structure, can analyze both quantitative and dichotomous traits and boosts the power of rare variant tests by incorporating multiple variant functional annotations. Through meta-analysis of four lipid traits in 30,138 ancestrally diverse samples from 14 studies of the Trans Omics for Precision Medicine (TOPMed) Program, we show that MetaSTAAR performs rare variant meta-analysis at scale and produces results comparable to using pooled data. Additionally, we identified several conditionally significant rare variant associations with lipid traits. We further demonstrate that MetaSTAAR is scalable to biobank-scale cohorts through meta-analysis of TOPMed WGS data and UK Biobank WES data of ~200,000 samples.

Extended Data Fig. 1 |. Quantile-quantile plots for gene-centric unconditional meta-analysis of lipid traits LDL-C, HDL-C, TG and TC using TOPMed WGS data (n = 30,138).

MetaSTAAR-O is a two-sided test. Different symbols represent the MetaSTAAR-O P values of different functional categories of individual genes (putative loss-of-function, missense, synonymous, promoter and enhancer). The promoter and enhancer of a gene are the promoter and the GeneHancer region that overlap with CAGE sites for a given gene, respectively (Methods). Four lipid traits were analyzed using MetaSTAAR-O: LDL-C, low-density lipoprotein cholesterol; HDL-C, high-density lipoprotein cholesterol; TG, triglycerides; and TC, total cholesterol.

Extended Data Fig. 2 |. Manhattan plots for gene-centric unconditional meta-analysis of lipid traits LDL-C, HDL-C, TG and TC using TOPMed WGS data (n = 30,138).

The horizontal line indicates the genome-wide MetaSTAAR-O P value threshold of $5.00\times {10}^{-7}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction ($0.05/\left(20,000\times 5\right)=5.00\times {10}^{-7}$). MetaSTAAR-O is a two-sided test. Different symbols represent the MetaSTAAR-O P values of different functional categories of individual genes (putative loss-of-function, missense, synonymous, promoter and enhancer). The promoter and enhancer of a gene are the promoter and the GeneHancer region that overlap with CAGE sites for a given gene, respectively (Methods). Four lipid traits were analyzed using MetaSTAAR-O: LDL-C, low-density lipoprotein cholesterol; HDL-C, high-density lipoprotein cholesterol; TG, triglycerides; TC, total cholesterol.

Extended Data Fig. 3 |. Scatterplots comparing gene-centric unconditional meta-analysis P values from MetaSTAAR-O with STAAR-O from the joint analysis of pooled individual-level data (STAAR-O-Pooled) of lipid traits LDL-C, HDL-C, TG and TC using TOPMed WGS data (n = 30,138).

Each dot represents a functional category of a gene with x-axis label being the $-{\text{log}}_{10}\left(P\right)$ of STAAR-O-Pooled and y-axis label being the $-{\text{log}}_{10}\left(P\right)$ of MetaSTAAR-O (n = 30,138). The horizontal and vertical lines indicate the genome-wide P value threshold of $5.00\times {10}^{-7}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction ($0.05/\left(20,000\times 5\right)=5.00\times {10}^{-7}$). Both MetaSTAAR and STAAR are two-sided tests. LDL-C, low-density lipoprotein cholesterol; HDL-C, high-density lipoprotein cholesterol; TG, triglycerides; TC, total cholesterol.

Extended Data Fig. 4 |. Scatterplot of P values comparing MetaSTAAR-O to Burden-MS, SKAT-MS and ACAT-V-MS (MS is short for MetaSTAAR) for quantitative and dichotomous traits when 15% of rare variants are causal variants.

In each simulation replicate, a 2-kb region was randomly selected as the signal region. Within each signal region, variants were randomly generated to be causal based on a multiple logistic model and on average there were 15% causal variants in the signal region. The effect sizes of causal variants were ${\beta }_j=c_0\left|{\text{log}}_{10}{\mathit{\text{MAF}}}_j\right|$. For quantitative traits, $c_0=0.07$; for dichotomous traits, $c_0=0.11$. All causal variants had positive effect sizes. Power was estimated as the proportion of the P values less than $\alpha ={10}^{-7}$ based on ${10}^4$ replicates. Burden-MS, SKAT-MS, ACAT-V-MS and MetaSTAAR-O are two-sided tests. Five studies were included in meta-analysis, each with a sample size of 10,000.

Figure 1 |. MetaSTAAR workflow.

a, Input data of MetaSTAAR for each study, including genotypes, phenotypes, covariates and sparse genetic relatedness matrix are prepared. b, Summary statistics, including individual variant score statistics, sparse weighted LD matrices and low-rank projection matrices accounting for covariate effects for each study are generated using MetaSTAARWorker. c, All rare variants in the merged variant list are functionally annotated and two types of variant sets are defined: gene-centric analysis by grouping variants into functional genomic elements for each protein-coding gene; and genetic region analysis using agnostic sliding windows. d, The MetaSTAAR-O P values for all variant sets defined in c are obtained. e, The conditional MetaSTAAR-O P values for all significant variant sets from d after adjusting for known variants are obtained and reported.