Meta-SAIGE: Scalable and Accurate Meta-Analysis for Rare Variants

Meta-SAIGE: Scalable and Accurate Meta-Analysis for Rare Variants | medRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-P4HH5NV'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search Meta-SAIGE: Scalable and Accurate Meta-Analysis for Rare Variants View ORCID Profile Eunjae Park , View ORCID Profile Kisung Nam , Seokho Jeong , Karl Keat , Dokyoon Kim , Vikas Bansal , Wei Zhou , Seunggeun Lee doi: https://doi.org/10.1101/2024.09.17.24313855 Eunjae Park 1 Graduate School Data Science, Seoul National University , Seoul, Republic of Korea 2 Bioinformatics and Systems Biology, University of California San Diego , La Jolla, CA, United States of America Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Eunjae Park Kisung Nam 1 Graduate School Data Science, Seoul National University , Seoul, Republic of Korea Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Kisung Nam Seokho Jeong 1 Graduate School Data Science, Seoul National University , Seoul, Republic of Korea Find this author on Google Scholar Find this author on PubMed Search for this author on this site Karl Keat 3 Informatics in Biostatistics and Epidemiology, University of Pennsylvania , Philadelphia, Pennsylvania, United States of America Find this author on Google Scholar Find this author on PubMed Search for this author on this site Dokyoon Kim 3 Informatics in Biostatistics and Epidemiology, University of Pennsylvania , Philadelphia, Pennsylvania, United States of America Find this author on Google Scholar Find this author on PubMed Search for this author on this site Vikas Bansal 2 Bioinformatics and Systems Biology, University of California San Diego , La Jolla, CA, United States of America Find this author on Google Scholar Find this author on PubMed Search for this author on this site Wei Zhou 4 Psychiatric and Neurodevelopmental Genetics Unit, Center for Genomic Medicine , Boston, Massachusetts, United States of America 5 Stanley Center for Psychiatric Research, Broad Institute of MIT and Harvard , Cambridge, Massachusetts, United States of America 6 Program in Medical and Population Genetics, Broad Institute of MIT and Harvard , Cambridge, Massachusetts, United States of America 7 Institute for Molecular Medicine Finland, University of Helsinki , Helsinki, Finland Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: wzhou{at}broadinstitute.org lee7801{at}snu.ac.kr Seunggeun Lee 1 Graduate School Data Science, Seoul National University , Seoul, Republic of Korea 8 Genomic Medicine Institute, Seoul National University , Seoul, Republic of Korea Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: wzhou{at}broadinstitute.org lee7801{at}snu.ac.kr Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract A meta-analysis is a practical approach to increasing the power of rare variant tests by combining summary statistics from multiple cohorts. However, existing methods for rare variant meta-analysis often fail to correctly control type I error rates when analyzing low-prevalence binary traits and are computationally intensive when analyzing many phenotypes. This paper introduces Meta-SAIGE, a novel approach for rare variant meta-analyses that addresses these challenges. Meta-SAIGE reduces type I error inflation through precise estimation of the distribution of test statistics and allows the reuse of the linkage disequilibrium (LD) matrix across phenotypes, significantly improving computational efficiency for phenome-wide analyses. Simulation studies using UK Biobank whole-exome sequencing (WES) genotypes demonstrate that Meta-SAIGE effectively controls type I error rates and yields power similar to that of pooled individual-level data through SAIGE-GENE+. A meta-analysis of UK Biobank and All of Us WES data for 83 low prevalence disease phenotypes identified 237 associations. Notably, 80 of these associations were not significant in either dataset alone, underscoring the power of our meta-analysis. Introduction The advent of mega biobanks and consortia has significantly enhanced our ability to sequence large cohorts, facilitating the identification of rare genetic variants linked to human diseases and traits. Due to low minor allele frequencies, single-variant tests are underpowered for identifying rare variant associations [ 1 – 3 ]. In response, gene or region-based association tests, which aggregate the effects of multiple rare variants, have been proposed. Different statistical tests such as the Burden test [ 4 ], SKAT [ 5 ], and SKAT-O [ 6 ] have been pivotal in the analysis of rare variant associations. More recently, STAAR [ 7 ], SAIGE-GENE [ 8 ] and SAIGE-GENE+ [ 9 ] have been introduced, offering efficient implementation of rare variant tests for large-scale biobank datasets, while accommodating sample relatedness and case–control imbalances. Meta-analysis is a practical approach for identifying genetic associations that may not be detectable in individual cohorts but become statistically significant when summary statistics from multiple studies are combined [ 10 ]. Meta-analysis is particularly relevant for rare variant association tests, where the low frequencies of the variants can limit test power in individual studies. Substantial efforts have been made to facilitate the meta-analysis of rare variants association studies. For example, international consortia such as the Biobank Rare Variant Analysis (BRaVa) have been established. As more biobanks carry out sequencing studies, rare-variant meta-analysis is likely to become a primary approach for identifying novel rare variant associations for human diseases and traits. Several methodologies for rare variant meta-analyses have been developed, notably RareMetal [ 11 , 12 ] and MetaSKAT [ 1 ]. More recently, MetaSTAAR [ 13 ] was introduced, enhancing the analysis by integrating variant functional annotations and accommodating sample relatedness. However, MetaSTAAR can exhibit significantly inflated type I error rates under imbalanced case–control ratios (Supplementary Note 1 and Supplementary Figure 1), a prevalent challenge in biobank-based disease phenotype studies. Additionally, MetaSTAAR’s requirement for constructing separate linkage disequilibrium (LD) matrices for each phenotype substantially increases the computational load when analyzing multiple phenotypes. An alternative approach is the Fisher method [ 14 ], which combines gene or region test p-values from each cohort, but it often has substantially lower power than joint analysis with individual-level data. This paper proposes a new rare variant meta-analysis method, Meta-SAIGE, which extends SAIGE-GENE+ [ 9 ] to meta-analysis. First, Burden, SKAT, and SKAT-O test p values are calculated using per-variant score statistics and LD matrices from each cohort. Next, to address type I error inflation in the presence of case–control imbalance, Meta-SAIGE employs two-level saddlepoint approximation (SPA) [ 15 , 16 ], including SPA on score statistics of each cohort, and a genotype count-based SPA for combined score statistics from multiple cohorts [ 17 ]. In addition, Meta-SAIGE allows the use of a single sparse LD matrix across all phenotypes, significantly reducing computational costs when conducting phenome-wide analyses involving hundreds or thousands of phenotypes. Simulation studies using UK Biobank (UKB) whole-exome sequencing (WES) data [ 18 ] demonstrate that Meta-SAIGE effectively controls type I error rates and is computationally faster and more efficient than MetaSTAAR. We conducted a meta-analysis of UKB and All of Us [ 19 ] exome-wide rare-variant gene-based tests for data for 83 disease phenotypes and identified 237 associations at the exome-wide significance level. Among these associations, 80 were not significant in either dataset alone, underscoring the power of our meta-analysis. Results Overview of Meta-SAIGE Meta-SAIGE contains three steps: (1) preparing per-variant level association summaries and a sparse LD matrix for each cohort, (2) combining summary statistics from all the studies into a single superset, and (3) running gene-based tests, as illustrated in Figure 1 . Download figure Open in new tab Fig. 1. Overview of Meta-SAIGE (A) Summary statistics generated by each cohort for Meta-SAIGE. (B) Meta-analysis incorporating per-variant summary statistics (score statistics S and its variance V ) and sparse LD matrices. The variance V is adjusted using SPA. Multiple summary statistics are merged into a superset, with additional variance adjustments via the GC-based SPA method. Ultrarare variants are collapsed, followed by gene-based tests using various functional annotations and minor allele frequency (MAF) cutoffs. P-values from these tests are combined using the Cauchy combination method. The first step uses SAIGE [ 16 ] to derive per-variant score statistics, S , accompanied by their variance and association p-values. SAIGE conducts the score test based on a generalized linear mixed model that can produce accurate single variant p-values by adjusting for case–control imbalance and sample relatedness ( Figure 1. A ). For binary phenotypes, the computation of p-values involves two possible methods, saddlepoint approximation (SPA) and efficient resampling [ 20 ], depending on the minor allele count (MAC). Additionally, the sparse LD matrix, Ω, is generated in each cohort, which is the pairwise cross product of dosages across genetic variants. In contrast to MetaSTAAR, which uses a phenotype variance weighted LD matrix, the sparse LD matrix used by Meta-SAIGE is not phenotype specific, allowing us to reuse the same LD matrix across different phenotypes. When meta-analyzing M variants from K cohorts, the storage requirement for a single phenotype for the sparse LD matrix is O ( MFK ), where F is the number of variants with the nonzero cross-product on average, and for single variant summary statistics, it is O ( MK ). When the meta-analyses are performed for P different phenotypes, our approach requires O ( MFK + MKP ) storage, while MetaSTAAR requires O ( MFKP + MKP ) storage. Next, score statistics from multiple cohorts are consolidated into a single statistic ( Figure 1. B ). For binary traits, the variance of each score statistic is recalculated by inverting the p-value by SAIGE-GENE+, which is a method that we previously developed for rare-variant association test. To further improve type I error control, we applied the Genotype-Count (GC) based saddlepoint approximation [ 17 ], which was previously designed to improve type I error control in the meta-analysis setup. The covariance matrix of score statistics is calculated in a sandwich form, Cov ( S ) = V 1 / 2 Cor ( G ) V 1 / 2 . Here Cor ( G ) is the correlation matrix of genotypes and is calculated with the sparse LD matrix Ω, and V is the diagonal matrix of the variance of S , which is obtained by inverting the SPA-GC-adjsuted p-values. With the combined per-variant statistics and covariance matrix, Meta-SAIGE performs rare-variant tests in the same manner as SAIGE-GENE+. It conducts Burden, SKAT, and SKAT-O set-based tests utilizing various functional annotations and MAF cutoffs. The ultrarare variants (e.g. MAC < 10) are collapsed to enhance type I error control and power, while reducing computation cost as demonstrated in SAIGE-GENE+. The Cauchy combination method [ 21 , 22 ] is then used to combine p-values corresponding to different functional annotations and MAF cutoffs for each testing gene or region. Type 1 Error and Power Evaluation To evaluate the type I error rates of Meta-SAIGE, simulations were conducted using the UKB WES data of 160,000 white British participants, which were partitioned into three non-overlapping cohorts of varying sizes. Null phenotypes were generated for each cohort, with prevalence of 5% and 1%. The simulation process was replicated 60 times to yield approximately one million tests. The empirical type I error rates are presented in Figure 2 and Supplementary Table 1. Without any adjustment, mirroring MetaSTAAR, there was a significant inflation in type I error rates. For example, with a disease prevalence of 1% and a sample size ratio of 1:1:1, the type I error rate for the noadjustment method at α = 2.5 × 10 − 6 was α = 2.12 × 10 − 4 , nearly 100 times higher than the nominal level (Supplementary Table 1). The SPA applied to score statistics of each cohort (SPA adjustment) notably decreased type I error rates, albeit with some remaining inflation. Meta-SAIGE, which further applies a GC-based SPA, outperformed other methods, demonstrating relatively well controlled type I error rates. Based on the same UKB genotype data, we further conducted simulation to assess the power of Meta-SAIGE in detecting rare-variant associations ( Figure 3 and Supplementary Table 2). We examined various scenarios of different effect sizes of rare variants and benchmarked Meta-SAIGE’s performance against a joint analysis of individual-level data via SAIGE-GENE+ and the weighted Fisher’s method [ 14 ], which aggregates SAIGE-GENE+ p-values weighted by sample size. Meta-SAIGE consistently demonstrated statistical power on par with that of joint analyses across all evaluated scenarios. In contrast, the weighted Fisher’s method yielded significantly lower power, underscoring the potential of our approach to enhance the detection of rare variant associations. Download figure Open in new tab Fig. 2. QQ-plots for null phenotypes with prevalence of 5% and 1% with different sample size ratios (1:1:1 and 4:3:2). (A) QQ-plot for a 1:1:1 sample size ratio with prevalence=5%. (B) QQ-plot for a 4:3:2 sample size ratio with prevalence=5%. (C) QQ-plot for a 1:1:1 sample size ratio with prevalence=1%. (D) QQ-plot for a 4:3:2 sample size ratio with prevalence=1%. Download figure Open in new tab Fig. 3. Power evaluation results for the binary phenotype with (A) 5% prevalence and (B) 1% prevalence. The scenarios are annotated with three indices, ( i, j, k ), which indicate the proportion of causal variant ( i ), magnitude of the effect sizes ( j ), direction of the effect of causal variants ( j ). Details can be found in Supplementary Table 2. The Y-axis represents the power ratio, with a value of one indicating equal power compared to the joint test. Computational Cost Evaluation We conducted a comparative analysis of the computational costs of Meta-SAIGE and MetaSTAAR. Because MetaSTAAR software lacks the ability to read genotype data directly, we utilized seqminer[ 23 ] for genotype data preprocessing. Figure 4 shows the computation time and memory usage for analyzing the largest gene, TTN (number of rare variants = 17,577, number of ultrarare variants = 15,098), across three key meta-analysis tasks: per-variant statistics generation, LD matrix generation, and the meta-analysis itself. The first two tasks should be performed independently for each cohort. Download figure Open in new tab Fig. 4. Computational cost evaluation results on TTN (A) Computation time for input file generation, including null model fitting and summary statistics generation. (B) Computation time for LD matrix generation. (C) Computation time for the meta-analysis step. Our results indicate that Meta-SAIGE consistently results in lower computational demands in terms of both CPU time and memory usage. For example, with a sample size of 179, 176, Meta-SAIGE was significantly faster than MetaSTAAR, with speed improvements of 8.1, 6.5, and 24.7 times for pervariant statistics generation, LD matrix generation, and meta-analysis, respectively. The difference in computational efficiency between Meta-SAIGE and MetaSTAAR diminished when smaller genes, such as MSH6 (number of rare variants = 990), were analyzed (Supplementary Figure 2). As previously mentioned, MetaSTAAR requires phenotype-specific LD information. Consequently, when the P phenotypes were analyzed, the LD matrix ( Figure 4. B ) must be generated P times for each cohort. The storage sizes of the LD matrix for whole-exome sequencing (WES) using Meta-SAIGE were 3.0 GB for UKB and 2.8 GB for All of Us. However, the projected LD matrix storage values for MetaSTAAR were 8.47 GB for UKB and 7.91 GB for All of Us. When 100 phenotypes were analyzed, the storage cost of LD matrices required by MetaSTAAR is projected to be 847 GB for UKB and 791 GB for All of Us, which is about 282 times higher than Meta-SAIGE. Therefore, by using single LD matrix for all phenotypes,Meta-SAIGE significantly reduces the storage costs. Meta-analysis of UKB and All of Us Datasets To illustrate the performance of Meta-SAIGE in real data, we applied Meta-SAIGE to conduct a wholeexome meta-analysis for 83 disease phenotypes using UKB and All of Us European Ancestry samples. Our analysis included 459K European Ancestry samples from UKB and 116K European Ancestry samples from All of Us. Traits were selected to be compatible with both UKB and All of Us with sufficient case sizes while incorporating traits from previous WES phenome-wide analyses[ 9 ]. Case–control ratios varied from 1:2 (essential hypertension, PheCODE=401.1) to 1:305 (hemorrhage from gastrointestinal ulcers, PheCODE=531.1) across all analyzed phenotypes. A detailed list of phenotypes is provided in Supplementary Table 3. We applied three MAF cutoffs (1%, 0.1%, 0.01%) and three functional annotations (LoF, LoF+Missense, LoF+Missense+Synonymous), resulting in a total of nine different test masks. Figure 5 A and Supplementary Table 4 present all significant associations at the exome-wide significance level of α = 2.5 × 10 − 6 among, with 20,000 protein coding genes tested. Across all the 83 phenotypes, 237 gene–phenotype pairs were identified as significant. The estimated false discovery rate (FDR) was 0.0175 (FDR = 83 [number of traits] × 20,000 [number of genes] × 2.5 × 10 − 6 [exome-wide significance threshold]/237 [number of significant genes]). When a more stringent criterion adjusted for the number of phenotypes at the top of the exome-wide significance level (p-value < 2.5 × 10 − 6 / 83 = 3.0 × 10 − 8 ) was applied, 150 gene–phenotype associations remained significant. Download figure Open in new tab Fig. 5. UKB + AllofUs meta-analysis results (A) Manhattan plot displaying the number of significantly associated traits. (B) Manhattan plot for type 2 diabetes (PheCode 250.2). (C) Manhattan plot for colorectal cancer (PheCode 153). In (B) and (C), p-values calculated using UK Biobank (UKB) data alone are shown in blue, while p-values from Meta-SAIGE combining UKB and All of Us data are shown in red. Among the 237 associations, 139 and 46 of them were significant in the UKB data and All of US data only, respectively, underscoring the power of the meta-analysis. Among all significant genes, TTN and PKD1 were the most pleiotropic one, being associated with 7 phenotypes, followed by CHEK2, TET2 , and SRSF2 , each associated with 5 phenotypes ( Figure 5. A ). To illustrate the type I error control and improvement in p-values through meta-analysis, we generated Manhattan plots for two phenotypes: type 2 diabetes (T2D) and colorectal cancer ( Figure 5 B ). A total of five genes were significant for colorectal cancer, and 8 genes were significant for T2D. Our results clearly demonstrate that the meta-analysis can substantially improve the power to detect disease associated genes. For example, Adenomatous polyposis coli ( APC ) is a tumor suppressor gene and defects in this gene increase the risk of colorectal cancer [ 24 ]. In our analysis, the association of APC with colorectal cancer was not significant in UKB or All of US data alone (UKB p-value=4.40 × 10 − 3 and All of US p-value=1.32 × 10 − 5 ), and the meta-analysis p-value was 5.39 × 10 − 12 , which was much more significant than in individual biobanks. Discussion This study introduces Meta-SAIGE, a method for meta-analyzing rare variant association studies. By employing SPA and GC-based adjustment methods, Meta-SAIGE accurately calculates p-values, effectively controlling type I error rates while preserving true signals. As demonstrated in the type I error simulation, substantial reductions in type I error rates were observed compared to the unadjusted approach. Additionally, Meta-SAIGE improves computation efficiency by allowing the use of a single LD matrix for different phenotypes, facilitating more cost-effective phenome-wide analyses as demonstrated through a meta-analysis of UKB and All of Us WES data. Our simulation results showed that the meta-analysis based on Meta-SAIGE cohort-level summary statistics by Meta-SAIGE has similar study power to the joint analysis using SAIGE-GENE+ on individual-level data. In real data analysis of UKB and All of Us WES data, highly concordant p-values were observed between the two approaches for three example phenotypes (Supplementary Note B) and both approaches revealed identical significant genes. Many biobanks are now collecting samples from diverse ancestries to mitigate the Eurocentric bias prevalent in genetic association studies, emphasizing the need for multi-ancestry meta-analysis. Although our study focused on a single ancestry, Meta-SAIGE can also be applied to multi-ancestry samples. A key challenge in multi-ancestry analysis is collapsing ultra-rare variants, as these variants can differ significantly between ancestry groups. To address this, we added a step for ancestry-specific collapsing (see Methods section). For T2D and colorectal cancer, we conducted a multi-ancestry meta-analysis of all European ancestry samples (EUR, n=572K) from UKB and All of Us, as well as African ancestry (n=46K) and American (n=36K) samples from All of Us. Only variants with a minor allele frequency (MAF) ≤ 0.01 across all ancestry groups were included for meta-analysis (Supplementary Figure 4). Supplementary Table 6 lists the significant genes identified by either European-only Meta-analysis (EUR-Meta) or Meta-analysis of Cross Ancestry (Meta-CrossAncestry). Meta-CrossAncestry had generally lower p-values than EUR-Meta and identified two additional associations for colorectal cancer: RNU1-31P and RPL29P11. RPL29P11, Ribosomal Protein L29 Pseudogene 11, is pseudogene of the ribosomal protein L29, which is shown to up-regulate colorectal cancer [ 25 , 26 ], demonstrating the enhanced power of multi-ancestry meta-analysis. However, a more comprehensive evaluation of its performance across broader and more varied ancestry groups is required, which we leave for future investigation. Overall, we demonstrated that Meta-SAIGE can effectively handle large-scale rare-variant metaanalyses of binary phenotypes, adjusting for phenotype imbalance. As many biobanks are currently generating large-scale sequencing data, our method can facilitate more powerful meta-analysis studies for identifying rare disease-associated genetic variants and elucidating the genetic architecture of underlying complex traits. Methods Workflow of Meta-SAIGE Meta-SAIGE comprises three steps: (1) preparation of single-variant association summaries for each cohort, (2) consolidation of the summaries from all studies into a singular comprehensive dataset, and (3) implementation of gene-based association testing, as illustrated in Figure 1 . Summary Statistics from Each Cohort The first step can be performed using SAIGE software to compute per-variant score statistics ( S ), variance ( V ), p-value, and a sparse LD-matrix (Ω). For binary phenotypes, case–control imbalance is accounted for using the SPA test applied to a logistic regression model as detailed by Dey et al.[ 27 ]. Here, Y i represents the phenotype of the i th sample, X i represents the nongenetic covariates and b represents a random effect in the logistic mixed model. Note that for continuous traits, alinear mixed model is used instead. With K cohorts, the score statistic for each variant j in the k th cohort is computed as where n k is the sample size of the k th cohort, is the covariate adjusted genotype, and is the estimate of P r ( Y i = 1| x i ) in the k th cohort under H 0 . The variance of s j,k (i.e. v j,k ) is calculated as where W k is a diagonal matrix with being the i th diagonal element. Because phenotype imbalance can lead to an inflation of type 1 error, saddlepoint approximation (SPA) [ 15 ] has been proposed to calculate accurate p values. The SPA uses a cumulant generating function (CGF), which can be defined as follows: and its second derivative can be used to calculate the p value of each variant. Now, the SPA-adjusted variance, , is calculated by inverting the SPA p value with the given score statistic S j,k . The sparse LD matrix (Ω) can also be generated as Combining Summary Statistics Into a Single Super-Set In Meta-SAIGE, summary statistics from multiple studies are consolidated into a single table, namely, a superset. The superset is subsequently used to perform gene-based tests. In this step, S j and V j from each study are combined as Suppose that a region with the first m variants ( j = 1, …, m ) is tested. To obtain the covariance matrix of , we first calculate the correlation matrix of S, Cor , using the sparse LD matrix Ω k and the MAF of each variant obtained from each respective cohort. Here MAC k represents the vector of minor allele counts in cohort k. Note that Ψ is the covariance matrix of S * and diag (Ψ) is a diagonal matrix of Ψ. Then, S * follows a multivariate normal distribution, given that Cor is the correlation among genetic markers where V * is a diagonal matrix with the diagonal element being . We note that the unadjusted approach uses the unadjusted variance V j ,1 to calculate Gene-Based Test with the Super-Set Prior to running the region-based test, ultrarare variants are collapsed. Ultrarare variants were determined using the overall MAC calculated from all cohorts, specifically when the overall MAC ≤ 10. The collapsing process aligns with the principles of the burden test, aggregating summary statistics for the selected variants. Subsequently, gene-based tests are performed. The test statistics Q for Burden and SKAT are as follows: where w j represents an MAF-based weighting vector calculated from the beta density function with Beta (1, 25). The SKAT-O test is subsequently performed by combining the Q Burden and the Q SKAT . For each gene, multiple tests were conducted with various MAF cutoffs and functional annotations. P-values are combined using the Cauchy combination method[ 21 ]. In this study, MAF cutoffs of 1%, 0.1%, and 0.01% and functional annotations of, LoF (loss of function), LoF+Missense, and LoF+Missense+synonymous were used. Further Adjustment with a GC-based Meta-Analysis Approach The SPA-adjustment previously outlined assumes that score statistics approximate a normal distribution following variance adjustment using SPA. However, for significantly imbalanced binary phenotypes, this adjustment alone does not suffice to maintain type I error rates within acceptable bounds. The cause of this issue lies in the discrete nature of study-specific statistics, for which mere aggregation proves inadequate. To address this, we previously introduced a genotype-count (GC)-based procedure that uses a GC-based SPA as a reference distribution instead of a normal distribution for meta-analysis [ 17 ]. The GC-based method estimates the cumulant generating function of the aggregated score statistics utilizing genotype counts to better approximate distribution, thereby enhancing the accuracy of the p-value estimation. This procedure is specifically applied to calculate the p value of when the single variant p value, derived from variance adjustment, falls below 0.05. We applied it to calculate the p value of when the single variant p value calculated from the variance adjustment was smaller than 0.05. is subsequently determined by inverting the p-value obtained through the genotype-count (GC) based approach. Simulation studies demonstrate that this GC-based meta-analysis approach markedly diminishes type I error by providing a more accurate estimate of the score distribution. Ancestry-Specific Collapsing Given that ultra-rare variants are population-specific, Meta-SAIGE can perform ancestry-specific ultrarare variant collapsing for multi-ancestry meta-analysis. For cohorts sharing the same ancestry, ancestryspecific collapsing is conducted prior to creating the super-set. These ancestry-specifically collapsed variants are treated as unique pseudo-variants. Type 1 Error and Power Evaluation For both type I error and power simulations, we used WES genotype data from 166,960 white British samples in UKB. For type I error evaluation, null phenotypes were generated using the same procedure described in [ 9 ]. The logistic mixed model was used to generate the null phenotypes, where α 0 is the intercept term determined by the phenotype prevalence (5% and 1%); X i 1 and X i 2 are simulated binary and continuous covariates from Bernoulli(0.5) and N (0, 1), respectively; is the standardized genotype value of the ith individual and jth marker, and γ j is the genetic effect size from N (0, 1 /L ). These markers represent polygenic effects. Only even chromosome markers were used for null phenotype generation if the analysis was conducted on odd chromosomes and vice versa. A total of L = 30, 000 LD pruned markers were used. Null phenotype generation was replicated 60 times to yield approximately one million tests. The samples were divided into three cohorts with sample sizes designed to reflect specific ratios. Related samples (pairwise values in GRM ≥ 0.05) were assigned to the same cohort. In the first scenario, samples were distributed relatively equally across three cohorts, assigning (55,655, 55,654, and 55,652) samples, respectively. In the second scenario, a 4:3:2 ratio was used, with values of (74,205, 55,654, 37,102). To assess the statistical power, we simulated binary phenotypes using ten specific genes. Specifically, a phenotype was generated from the following model, where π i 0 is the probability for the i th individual being a case given the covariates and genotype, is the standardized genotype value of the j th variant of the i th individual, β j is the genetic effect size following N (0, 1 /L ), g ikq is the genotype value of the q th rare variant in gene k of the i th individual, and β kq is the genetic effect size of the q th rare variant in gene k according to each scenario described below. Like in the type I error simulation, two covariates X i 1 and X i 2 were generated from Bernoulli(0.5) and N (0, 1), respectively, and the intercept α 0 was determined to set the overall disease prevalences at 1% and 5%. Our investigation encompassed eight distinct scenarios, each characterized by the proportions of causal variants, effect size magnitudes, and effect directions of rare variants. Specifically, we considered two different configurations of causal variant proportions for non-ultra-rare (MAC > 10) variants: (1) 20% of LoF, 10% of missense, and 2% of synonymous and (2) 30% of LoF, 20% of missense and 2% of synonymous. Given that functionally deleterious variants are more likely to be rare, we posited that the proportions of causal variants among ultrarare variants (MAC ≤ 10) were threefold greater than those among less rare variants across various functional annotations. We explored three different settings for absolute effect sizes: (1) |0.5 log 10 ( MAF )| for LoF and |0.25 log 10 ( MAF )| for missense and synonymous and (2) |0.7 log 10 ( MAF )| for LoF and |0.35 log 10 ( MAF )| for missense and synonymous. We also examined two alternative configurations for effect directions among causal variants: (1) all causal variants exhibited the same directional effect, and (2) 100% LoF, 80% missense, and 50% synonymous variants increased the risk of disease, whereas the remaining causal variants decreased the risk of disease. For each simulation scenario, we conducted 100 independent repetitions, resulting in a total of 1,000 gene-level p-values. We partitioned 166,960 White British samples evenly distributed across three cohorts, assigning 55,655, 55,654, and 55,652 samples, respectively. UKB and All of Us WES Data Meta-analysis We performed a meta-analysis on data from the UKB and All of Us datasets for 83 binary traits among individuals of European ancestry (Supplementary Table 3). Binary traits were phenotyped based on Phecodes [ 28 ], which were derived from ICD records. Per-variant summary statistics and LD matrix calculations and meta-analyses were performed on the UKB Research Analysis Platform and the All of Us Workbench cloud system. WES data from 459,369 white UK Biobank participants were used to create association summaries. Hospital inpatient data coded in the ICD-10 were transformed into Phecodes, and participants were categorized into case–control groups accordingly. PLINK binary files of whole-exome sequences were used for the analysis. Genomic and electronic health record (EHR) data from All of Us were retrieved from the v7 Curated Data Repository (CDR) on the controlled tier. The multiancestry grouping of All of Us participants was derived from the Global Diversity Array (GDA), where principal components were used for ancestry labels [ 19 ]. Whole genome sequencing data from All of Us in BGEN format were used to create per-variant summary statistics. EHR data based on ICD-9 and ICD-10 codes were used to map Phecodes. Samples without EHR data were excluded from the analysis. For both the UKB and All of Us datasets, we included age, sex, batch information, and 10 ancestry principal components as covariates. Summary statistics from the UK Biobank and All of Us were generated with the maximum minor allele frequency of 0.01. To create the linkage disequilibrium (LD) matrix, variants pooled from both cohorts were annotated using VEP105 + LOFTEE GRCh 38 v. 1.04 ( https://github.com/BRaVa-genetics/vep105_loftee ) into three functional categories: LoF, missense, and synonymous variants. Additionally, we performed a meta-analysis on multiple ancestry groups across UKB and All of US on Type 2 diabetes and colorectal cancer. Multi-ancestry grouping of All of US participant is derived from Global Diversity Array(GDA), where principal components were generated and projected onto the ancestry label provided by GnomAD. Ancestry groups of over 10,000 samples from All of Us, European(EUR, n = 116,566), African(AFR, n = 45,772), and American(AMR, n=36,152), were used for this analysis. Per-variant summary statistics and LD matrix of UKB-White, All of Us EUR, AFR, and AMR groups were generated using SAIGE and applied to Meta-SAIGE. Data Availability Meta-SAIGE is implemented as an open-source R package available at https://github.com/leelabsg/META_SAIGE. Single variant summary statistics and LD matrix can be generated from SAIGE available at https://github.com/saigegit/SAIGE. UKB and All of US European meta-analysis results are available at https://meta-saige.leelabsg.org/. Computational Cost Evaluation We evaluated the computational cost of MetaSTAAR and Meta-SAIGE across different sample sizes: 10,000, 50,000, 100,000, and 179,176 samples. Each sample size was further divided into three sub-cohorts. For comparison, we primarily used the TTN gene, the largest gene, which represents a computational bottleneck for the analysis. We included LoF, missense, and synonymous variants with a minor allele frequency (MAF) ≤ 0.001, resulting in a total of 17,577 variants. Additionally, we analyzed a smaller gene, MSH6, which had 990 variants. We measured CPU time and memory usage across three distinct stages of the analysis: (1) input file generation, including fitting the null model and conducting single variant association tests; (2) linkage disequilibrium (LD) matrix generation; and (3) meta-analysis. Code and Data Availability Meta-SAIGE is implemented as an open-source R package available at https://github.com/leelabsg/META_SAIGE . Single variant summary statistics and LD matrix can be generated from SAIGE available at https://github.com/saigegit/SAIGE . UKB and All of US European meta-analysis results are available at https://meta-saige.leelabsg.org/ . URLS Meta-SAIGE: https://github.com/leelabsg/META_SAIGE SAIGE: https://github.com/saigegit/SAIGE UKB+AllofUs meta-analysis PheWeb: https://meta-saige.leelabsg.org/ SPAtest: https://github.com/leeshawn/SPAtest VEP105 + LOFTEE GRCh 38 v. 1.04 : https://github.com/BRaVa-genetics/vep105_loftee Inclusions and Ethics Statement Our study adhered to ethical conduct and inclusive practices throughout. All collaborators meeting Nature Portfolio journals’ authorship criteria are included as authors. We utilized data from All of Us and UKB, respecting their established data access protocols. Our research complied with relevant ethical regulations for secondary data analysis and posed no additional risks to participants. We implemented rigorous data management and analysis protocols. We maintained responsible practices in data interpretation, acknowledging the strengths and limitations of the datasets used. Author Contributions E.P. and S.L. implemented the Meta-SAIGE software. E.P. conducted the evaluation of type 1 error inflation. K.N. carried out the power evaluation. S.J. conducted meta-analysis using data from the UK Biobank and All of Us, with input from D.K. and K.K. Helpful advice was provided by D.K. and V.B. The manuscript was written by E.P., K.N., S.J., W.Z., and S.L. Acknowledgements This research was supported by the Brain Pool Plus (BP+, Brain Pool+) Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT (2020H1D3A2A03100666). We thank the helpful discussions from Dr. Duncan Palmer and the “Biobank Rare Variant Analysis” (BRaVa) consortium. UK Biobank data were accessed under the accession number UKB: 45227. References [1]. ↵ Lee , S. , Abecasis , G. R. , Boehnke , M. & Lin , X. Rare-variant association analysis: study designs and statistical tests . The American Journal of Human Genetics 95 , 5 – 23 ( 2014 ). OpenUrl CrossRef PubMed [2]. Auer , P. L. & Lettre , G. Rare variant association studies: considerations, challenges and opportunities . Genome medicine 7 , 1 – 11 ( 2015 ). OpenUrl [3]. ↵ Tang , Z.-Z. & Lin , D.-Y. Meta-analysis for discovering rare-variant associations: statistical methods and software programs . The American Journal of Human Genetics 97 , 35 – 53 ( 2015 ). OpenUrl CrossRef PubMed [4]. ↵ Morgenthaler , S. & Thilly , W. G. A strategy to discover genes that carry multi-allelic or monoallelic risk for common diseases: a cohort allelic sums test (cast) . Mutation Research/Fundamental and Molecular Mechanisms of Mutagenesis 615 , 28 – 56 ( 2007 ). OpenUrl [5]. ↵ Wu , M. C. et al. Rare-variant association testing for sequencing data with the sequence kernel association test . The American Journal of Human Genetics 89 , 82 – 93 ( 2011 ). OpenUrl CrossRef PubMed [6]. ↵ Lee , S. et al. Optimal unified approach for rare-variant association testing with application to small-sample case-control whole-exome sequencing studies . The American Journal of Human Genetics 91 , 224 – 237 ( 2012 ). OpenUrl CrossRef PubMed [7]. ↵ Li , X. et al. Dynamic incorporation of multiple in silico functional annotations empowers rare variant association analysis of large whole-genome sequencing studies at scale . Nature genetics 52 , 969 – 983 ( 2020 ). OpenUrl CrossRef PubMed [8]. ↵ Zhou , W. et al. Scalable generalized linear mixed model for region-based association tests in large biobanks and cohorts . Nature genetics 52 , 634 – 639 ( 2020 ). OpenUrl CrossRef PubMed [9]. ↵ Zhou , W. et al. Saige-gene+ improves the efficiency and accuracy of set-based rare variant association tests . Nature genetics 54 , 1466 – 1469 ( 2022 ). OpenUrl CrossRef [10]. ↵ Evangelou , E. & Ioannidis , J. P. Meta-analysis methods for genome-wide association studies and beyond . Nature Reviews Genetics 14 , 379 – 389 ( 2013 ). OpenUrl CrossRef PubMed [11]. ↵ Feng , S. , Liu , D. , Zhan , X. , Wing , M. K. & Abecasis , G. R. Raremetal: fast and powerful meta-analysis for rare variants . Bioinformatics 30 , 2828 – 2829 ( 2014 ). OpenUrl CrossRef PubMed [12]. ↵ Liu , D. J. et al. Meta-analysis of gene-level tests for rare variant association . Nature genetics 46 , 200 – 204 ( 2014 ). OpenUrl CrossRef PubMed [13]. ↵ Li , X. et al. Powerful, scalable and resource-efficient meta-analysis of rare variant associations in large whole genome sequencing studies . Nature genetics 55 , 154 – 164 ( 2023 ). OpenUrl [14]. ↵ Yoon , S. , Baik , B. , Park , T. & Nam , D. Powerful p-value combination methods to detect incomplete association . Scientific reports 11 , 6980 ( 2021 ). OpenUrl [15]. ↵ Daniels , H. E. Saddlepoint approximations in statistics . The Annals of Mathematical Statistics 631–650 ( 1954 ). [16]. ↵ Zhou , W. et al. Efficiently controlling for case-control imbalance and sample relatedness in large-scale genetic association studies . Nature genetics 50 , 1335 – 1341 ( 2018 ). OpenUrl CrossRef PubMed [17]. ↵ Dey , R. et al. Robust meta-analysis of biobank-based genome-wide association studies with unbalanced binary phenotypes . Genetic epidemiology 43 , 462 – 476 ( 2019 ). OpenUrl [18]. ↵ Backman , J. D. et al. Exome sequencing and analysis of 454,787 uk biobank participants . Nature 599 , 628 – 634 ( 2021 ). OpenUrl CrossRef PubMed [19]. ↵ Bick , A. G. et al. Genomic data in the all of us research program . Nature ( 2024 ). [20]. ↵ Lee , S. , Fuchsberger , C. , Kim , S. & Scott , L. An efficient resampling method for calibrating single and gene-based rare variant association analysis in case–control studies . Biostatistics 17 , 1 – 15 ( 2016 ). OpenUrl CrossRef PubMed [21]. ↵ Liu , Y. & Xie , J. Cauchy combination test: a powerful test with analytic p-value calculation under arbitrary dependency structures . Journal of the American Statistical Association 115 , 393 – 402 ( 2020 ). OpenUrl CrossRef [22]. ↵ Liu , Y. et al. Acat: a fast and powerful p value combination method for rare-variant analysis in sequencing studies . The American Journal of Human Genetics 104 , 410 – 421 ( 2019 ). OpenUrl CrossRef PubMed [23]. ↵ Yang , L. , Jiang , S. , Jiang , B. , Liu , D. J. & Zhan , X. Seqminer2: an efficient tool to query and retrieve genotypes for statistical genetics analyses from biobank scale sequence dataset . Bioinformatics 36 , 4951 – 4954 ( 2020 ). OpenUrl [24]. ↵ Zhang , L. & Shay , J. W. Multiple roles of apc and its therapeutic implications in colorectal cancer . JNCI: Journal of the National Cancer Institute 109 , djw332 ( 2017 ). OpenUrl [25]. ↵ El Khoury , W. & Nasr , Z. Deregulation of ribosomal proteins in human cancers . Bioscience Reports 41 , BSR20211577 ( 2021 ). OpenUrl CrossRef PubMed [26]. ↵ Liu , J.-J. , Huang , B. H. , Zhang , J. , Carson , D. D. & Hooi , S. C. Repression of hip/rpl29 expression induces differentiation in colon cancer cells . Journal of cellular physiology 207 , 287 – 292 ( 2006 ). OpenUrl CrossRef PubMed Web of Science [27]. ↵ Dey , R. , Schmidt , E. M. , Abecasis , G. R. & Lee , S. A fast and accurate algorithm to test for binary phenotypes and its application to phewas . The American Journal of Human Genetics 101 , 37 – 49 ( 2017 ). OpenUrl CrossRef PubMed [28]. ↵ Denny , J. C. et al. Phewas: demonstrating the feasibility of a phenome-wide scan to discover gene–disease associations . Bioinformatics 26 , 1205 – 1210 ( 2010 ). OpenUrl CrossRef PubMed Web of Science View the discussion thread. Back to top Previous Next Posted September 19, 2024. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about medRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Meta-SAIGE: Scalable and Accurate Meta-Analysis for Rare Variants Message Subject (Your Name) has forwarded a page to you from medRxiv Message Body (Your Name) thought you would like to see this page from the medRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Meta-SAIGE: Scalable and Accurate Meta-Analysis for Rare Variants Eunjae Park , Kisung Nam , Seokho Jeong , Karl Keat , Dokyoon Kim , Vikas Bansal , Wei Zhou , Seunggeun Lee medRxiv 2024.09.17.24313855; doi: https://doi.org/10.1101/2024.09.17.24313855 Share This Article: Copy Citation Tools Meta-SAIGE: Scalable and Accurate Meta-Analysis for Rare Variants Eunjae Park , Kisung Nam , Seokho Jeong , Karl Keat , Dokyoon Kim , Vikas Bansal , Wei Zhou , Seunggeun Lee medRxiv 2024.09.17.24313855; doi: https://doi.org/10.1101/2024.09.17.24313855 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Genetic and Genomic Medicine Subject Areas All Articles Addiction Medicine (574) Allergy and Immunology (865) Anesthesia (304) Cardiovascular Medicine (4462) Dentistry and Oral Medicine (445) Dermatology (383) Emergency Medicine (611) Endocrinology (including Diabetes Mellitus and Metabolic Disease) (1517) Epidemiology (15251) Forensic Medicine (31) Gastroenterology (1132) Genetic and Genomic Medicine (6621) Geriatric Medicine (669) Health Economics (1002) Health Informatics (4564) Health Policy (1372) Health Systems and Quality Improvement (1617) Hematology (544) HIV/AIDS (1272) Infectious Diseases (except HIV/AIDS) (15938) Intensive Care and Critical Care Medicine (1107) Medical Education (624) Medical Ethics (147) Nephrology (670) Neurology (6643) Nursing (346) Nutrition (1001) Obstetrics and Gynecology (1149) Occupational and Environmental Health (957) Oncology (3350) Ophthalmology (981) Orthopedics (369) Otolaryngology (421) Pain Medicine (436) Palliative Medicine (130) Pathology (665) Pediatrics (1698) Pharmacology and Therapeutics (694) Primary Care Research (714) Psychiatry and Clinical Psychology (5465) Public and Global Health (9259) Radiology and Imaging (2212) Rehabilitation Medicine and Physical Therapy (1372) Respiratory Medicine (1198) Rheumatology (598) Sexual and Reproductive Health (716) Sports Medicine (533) Surgery (715) Toxicology (100) Transplantation (289) Urology (265) (function(){function c(){var b=a.contentDocument||a.contentWindow.document;if(b){var d=b.createElement('script');d.innerHTML="window.__CF$cv$params={r:'a03cedac1f68c13d',t:'MTc4MDEzNTYxMA=='};var a=document.createElement('script');a.src='/cdn-cgi/challenge-platform/scripts/jsd/main.js';document.getElementsByTagName('head')[0].appendChild(a);";b.getElementsByTagName('head')[0].appendChild(d)}}if(document.body){var a=document.createElement('iframe');a.height=1;a.width=1;a.style.position='absolute';a.style.top=0;a.style.left=0;a.style.border='none';a.style.visibility='hidden';document.body.appendChild(a);if('loading'!==document.readyState)c();else if(window.addEventListener)document.addEventListener('DOMContentLoaded',c);else{var e=document.onreadystatechange||function(){};document.onreadystatechange=function(b){e(b);'loading'!==document.readyState&&(document.onreadystatechange=e,c())}}}})();