Novel Variance-Component TWAS method for studying complex human diseases with applications to Alzheimer's dementia.

Transcriptome-wide association studies (TWAS) have been widely used to integrate transcriptomic and genetic data to study complex human diseases. Within a test dataset lacking transcriptomic data, traditional two-stage TWAS methods first impute gene expression by creating a weighted sum that aggregates SNPs with their corresponding cis-eQTL effects on reference transcriptome. Traditional TWAS methods then employ a linear regression model to assess the association between imputed gene expression and test phenotype, thereby assuming the effect of a cis-eQTL SNP on test phenotype is a linear function of the eQTL's estimated effect on reference transcriptome. To increase TWAS robustness to this assumption, we propose a novel Variance-Component TWAS procedure (VC-TWAS) that assumes the effects of cis-eQTL SNPs on phenotype are random (with variance proportional to corresponding reference cis-eQTL effects) rather than fixed. VC-TWAS is applicable to both continuous and dichotomous phenotypes, as well as individual-level and summary-level GWAS data. Using simulated data, we show VC-TWAS is more powerful than traditional TWAS methods based on a two-stage Burden test, especially when eQTL genetic effects on test phenotype are no longer a linear function of their eQTL genetic effects on reference transcriptome. We further applied VC-TWAS to both individual-level (N = ~3.4K) and summary-level (N = ~54K) GWAS data to study Alzheimer's dementia (AD). With the individual-level data, we detected 13 significant risk genes including 6 known GWAS risk genes such as TOMM40 that were missed by traditional TWAS methods. With the summary-level data, we detected 57 significant risk genes considering only cis-SNPs and 71 significant genes considering both cis- and trans- SNPs, which also validated our findings with the individual-level GWAS data. Our VC-TWAS method is implemented in the TIGAR tool for public use.

Introduction

Genome-wide association studies (GWAS) have succeeded in identifying thousands of genetic loci associated with complex traits and diseases [1-3]. However, for the most part the molecular mechanisms linking these genes with these complex traits and diseases remain unexplained [4]. Studies have shown that gene expression plays a key role in the phenotypic manifestation of human diseases [5]. Many common genetic variants associated with the phenotypes manifested by human diseases are highly likely to be expression quantitative trait loci (eQTL) [6,7]. Therefore, integrating gene expression data together with genetic data from GWAS is expected to help identify novel risk genes as well as elucidate the mechanisms underlying the associations of genetic loci with complex traits and disease phenotypes.

The transcriptome-wide association study (TWAS) is an innovative strategy to leverage the enhanced power of integrating gene expression of the transcriptome together with genetic data from GWAS [8] for gene-based association studies of complex traits. Using a reference transcriptomic panel like Genotype-Tissues Expression (GTEx) [9], traditional TWAS methods first train a regression model that treat gene expression as the outcome and SNP genotype data (generally cis-SNPs nearby the test gene) as predictors, which can be viewed as a gene expression imputation model. traditional two-stage TWAS tools such as PrediXcan [10], FUSION [11], and TIGAR [8] employ different regression methods to fit such gene expression imputation models, where corresponding cis-eQTL effect sizes are derived. For example, PrediXcan method [10] implements the Elastic-Net penalized regression method [12] while TIGAR [8] implements the nonparametric Bayesian Dirichlet process regression (DPR) [13] method. Regardless of technique, these existing two-stage TWAS methods produce a set of estimates of cis-eQTL effects on reference transcriptome. Within a GWAS lacking transcriptomic data, existing two-stage TWAS tools then imputes genetically-regulated gene expression (GReX) by weighted summing cis-eQTL SNP data from the GWAS with corresponding estimates of cis-eQTL effects on reference transcriptome. Once created, existing two-stage TWAS tools test for association between the imputed GReX and traits of interest using GWAS data, where the test is based on a single variant linear regression model. TWAS methods have been deployed widely since their inception and have had substantial success in improving our understanding the genetic regulation of various complex traits [14].

Since existing two-stage TWAS tools derive GReX as the summation of eQTL SNP genotypes weighted by their corresponding effects on reference transcriptome, they essentially conduct weighted burden tests [8] to evaluate the association between GReX and outcome (referred to as Burden-TWAS in this paper). Burden-TWAS methods assume the effects of eQTL SNPs on phenotype is a linear function of their corresponding estimated effects on reference transcriptome (see Methods), which may not be true in practice. In particular, the effect of eQTL SNP on phenotype may be a non-linear, rather than a linear, function of the eQTL effect on reference transcriptome. Moreover, the assumed eQTL effects on transcriptome in the test dataset could be mis-specified due to ancestral differences between the test and reference datasets. For example, recent studies show that gene-expression prediction models trained in one population perform poorly when applied to a different population [15]. The non-linearity assumption about the effects of eQTL SNPs on the phenotype of interest does not affect the validity of TWAS but it can reduce the power of these traditional two-stage Burden-TWAS methods.

Here, we derive a Variance-Component TWAS (VC-TWAS) method that relaxes the linearity assumption of Burden-TWAS. This modification makes this method more robust to misspecification of eQTL effect size estimates derived from the reference transcriptome, thus improving its power compared to the Burden-TWAS methods. VC-TWAS is a TWAS analogue of the Sequence Kernel Association Test (SKAT) [16-19] used for gene-based association studies. Unlike Burden-TWAS methods, our VC-TWAS aggregates genetic information across test SNPs using a kernel similarity function that allows upweighting or downweighing of specific variants in the similarity score based on eQTL effect size magnitudes (which can be derived using DPR methods in TIGAR or Elastic-Net penalized regression methods implemented in PrediXcan). The test statistic can be thought of as a Variance-Component score statistic based on a mixed model where each test variant has a random effect whose variance is a linear function of the squared values of corresponding eQTL effect size.

By modeling variants with random effects, the technique relaxes the main assumption of Burden-TWAS and is robust to misspecification of eQTL effect size estimates derived from the reference transcriptome (in terms of the direction and magnitude). The variance component score test employed by VC-TWAS also enjoys efficient computation and flexibility for studying both quantitative and dichotomous phenotypes, which enables VC-TWAS to consider both cis- and trans- eQTL SNPs. VC-TWAS is further applicable to both individual-level and summary-level GWAS data.

We note that VC-TWAS is not the first TWAS method developed that relaxes the linearity assumption of Burden-TWAS. A recently derived collaborative mixed model (CoMM) [20,21] likewise accounts for the uncertainty of eQTL effect size estimates from reference transcriptome data by jointly modeling reference and test data within a linear mixed-model framework. The PMR-Egger [22] and moPMR-Egger [23] also take the approach of jointly modeling reference and test data to address the uncertainty of eQTL effect size estimates from the reference data. However, these methods that are based on the maximum likelihood inference framework and implement likelihood ratio tests, are derived only for quantitative phenotypes and could suffer computation burden for testing thousands of cis-SNPs per gene, which often happens in practice particularly when considering imputed SNP data or whole-genome sequencing data for TWAS. In this study, we show that the likelihood ratio test approach used by CoMM and PMR-Egger is computationally more expensive to run than VC-TWAS in practice when there are thousands of test SNPs such as using imputed or Whole Genome Sequencing (WGS) genotype data.

This manuscript is as organized as follows. First, we provide an overview of traditional two-stage Burden-TWAS as well as existing techniques for estimating cis-eQTL SNP effects from reference transcriptome data. We then derive VC-TWAS approach for using both individual-level and summary-level GWAS data and compare its performance with Burden-TWAS using simulated data generated under different eQTL model assumptions. We then apply VC-TWAS to individual-level GWAS data from Religious Order Study and Memory Aging Project (ROS/MAP) [24-27] and Mayo Clinic Late-Onsite Alzheimer’s disease (LOAD) [28,29] cohorts, as well as the summary-level GWAS data from the International Genomics of Alzheimer’s Project (IGAP) [3], for studying Alzheimer’s dementia (AD). By considering only cis-eQTL effect sizes, VC-TWAS identified both novel and known risk genes for AD within 2MB of the well-known major risk gene APOE of AD, including the known risk gene TOMM40 and APOE. Considering both cis- and trans- eQTL effect sizes estimated by the Bayesian Genome-wide TWAS (BGW-TWAS) method [30], VC-TWAS detected 71 risk genes for AD which complemented existing TWAS results using only cis-eQTL SNP data. After describing the results, we provide a brief discussion summarizing our findings and describing implementation of VC-TWAS into our previously developed software tool TIGAR [8] for public use.

Methods

Ethics statement

Real ROS/MAP and Mayo Clinic GWAS data analyzed in this study were generated under the improvement by the Institutional Review Board (IRB) of Rush University Medical Center, Chicago, IL and Mayo Clinic, respectively. All samples analyzed in this study were de-identified and all analyses were approved by the IRB of Emory University School of Medicine.

Traditional Two-Stage TWAS procedure

Two-stage TWAS first fits gene expression imputation models by taking genotype data as predictors and assuming the following additive genetic model for expression quantitative traits, Here, is the genotype matrix for all considered SNP genotypes (encoded as the number of minor alleles or genotype dosages of SNPs within 1MB of the target gene region), is the eQTL effect size vector, and is the profiled gene expression levels for the target gene . Given the eQTL effect size estimates from reference data, will be imputed by the following equation where is the genotype matrix for the test cohort.

The general test framework of Burden-TWAS [8,10,11] that test for association between and the phenotype of interest can be written as: where denotes the phenotype of interest, g(.) denotes a link function, is imputed gene expression levels, and ′ denotes the coefficient vector for other non-genomic covariates . Basically, Burden-TWAS tests the null hypothesis of H0: β = 0, where eQTL effect size estimates () are taken as variant weights and SNP effect sizes on phenotype () are assumed to be a linear function of [8,10,11]. As noted above this linear relationship between eQTL effect sizes on reference transcriptome and the SNP effect sizes on test phenotype may not be true when analyzing human data and limits the power of the existing two-stage Burden-TWAS methods.

Estimation of eQTL effect sizes

Different methods can be used to estimate eQTL effect sizes from Eq (1). In this study, we applied PrediXcan (Elastic-Net penalized regression) and TIGAR (nonparametric Bayesian DPR methods) [8,10] to estimate that only consider cis-eQTL. PrediXcan TWAS method [10] employs Elastic-Net penalized regression method [12] to estimate cis-eQTL effect sizes from Eq (1) (S1A Text). TIGAR8 provides a more flexible approach to nonparametrically estimate cis-eQTL effect sizes by a Bayesian DPR method [13] (S1B Text). Additionally, we also considered modeling gene expression using both cis- and trans- eQTL effect sizes estimated by the recently proposed Bayesian genome-wide TWAS (BGW-TWAS) method [30].

VC-TWAS with individual-level GWAS data

Here, we propose a powerful VC-TWAS method that is analogous to the previously proposed SKAT method for SNP-set based association studies [16]. Similar to SKAT, the general test framework of VC-TWAS can be written as for continuous quantitative traits, and for dichotomous traits of sample i. Here, is the genetic effect size vector, is the genotype matrix for all test SNPs with respect to the test gene, is the non-genomic covariate matrix, and is the error term. VC-TWAS will test H0: τ = 0, which is equivalent to testing H0: = 0. The Variance-Component score statistic used by VC-TWAS is given by where is the estimated phenotype mean under and with weight w for the jth variant.

In contrast to SKAT methods, VC-TWAS takes eQTL effect size estimates from Eq (1) as variant weights (w). That is, the variances () of SNP effect sizes on phenotype are assumed to be a linear function of eQTL effect size estimates, which is robust to both direction and magnitude of eQTL effect size estimates. Since the Variance-Component score statistic Q (Eq (4)) follows a mixture of chi-square distributions under the null hypothesis [31,32], p-value can be conveniently obtained from several approximation and exact methods like the Davies exact method [33].

We note that our VC-TWAS method is computationally more expensive than standard Burden-TWAS given the need to perform eigen-decomposition of the kernel matrix in Eq 4 to obtain an analytic p-value (S1C Text). Such eigen-decomposition has computational complexity O(m3) for considering m SNPs with non-zero eQTL effect sizes. As DPR method produces non-zero cis-eQTL effect size estimates for almost all test SNPs (with most cis-eQTL effect size estimates being close to zero [8]), we explored an alternate VC-TWAS that considered a reduced set of SNPs by filtering out those with cis-eQTL effect size estimates smaller than the median cis-eQTL effect size estimate per gene. As we show using simulated data, we can reduce up to 80% computation time while having negligible impact on performance relative to using all SNPs with non-zero cis-eQTL effect size estimates.

Compared to Burden-TWAS, VC-TWAS is based on a random effect model and test for if the variation of random effects differs from 0. That is, the variance component score statistic (Eq (4)) used by VC-TWAS does not directly model the directions of cis-eQTL effect sizes (w) and is robust for mis-specifications of cis-eQTL effect sizes.

VC-TWAS with summary-level GWAS data

Since summary-level GWAS data generally provide SNP effect sizes on phenotype and corresponding standard errors based on single variant tests (often meta-analysis), it is reasonable to assume that the test phenotypes were adjusted for the confounding covariates. That is, in the Q statistic used by VC-TWAS as in Eq (4), we can assume the phenotype mean under is 0. This leads to a simplified formula [34,35] where is the single variant score statistic of the j variant, is the estimated phenotype variance.

As derived in the previous studies [36], we can estimate phenotype variance and then the score statistics by using only GWAS summary statistics including the single variant effect size estimate and corresponding standard error for j SNP, sample size n, and a reference LD covariance matrix Σ of all test SNPs (S1C Text). That is, the numerator and denominator of the score statistic can be estimated by,

ROS/MAP data

In this study we used clinical and postmortem data from older adults participating in two ongoing community-based cohorts studies the Religious Orders Study (ROS) and Rush Memory and Aging Project (MAP) [24-27] which document risk factors and chronic conditions of aging including dementia. Participants are senior adults without known dementia, who agree to annual clinical evaluation and brain autopsy at the time of death. All participants signed an informed consent and Anatomic Gift Act, and the studies are approved by an Institutional Review Board of Rush University Medical Center, Chicago, IL. All participants in this study also signed a repository consent to allow their data to be re-purposed.

Currently, we have microarray genotype data generated for 2,093 European-decent subjects from ROS/MAP [24-27], which are further imputed to the 1000 Genome Project Phase 3 [37] in our analysis [38]. Post-mortem brain samples (gray matter of the dorsolateral prefrontal cortex) from ~30% these ROS/MAP participants with assayed genotype data are also profiled for transcriptomic data by next-generation RNA-sequencing [39], which are used as reference data to train GReX prediction models in our application studies.

Using ROS/MAP data, we conducted TWAS for clinical diagnosis of late on-site Alzheimer’s dementia (LOAD) as well as three pathology indices of AD including quantified β-amyloid load and PHFtau tangle density as well as a summary measure of AD pathology burden (S1D Text) [24,25,27].

Mayo Clinic LOAD GWAS data

Mayo Clinic LOAD GWAS data contain samples from two clinical AD Case-Control series: Mayo Clinic Jacksonville (MCJ: 353 AD cases and 331 Controls), Mayo Clinic Rochester (MCR: 245 AD cases and 701 Controls) and a neuropathological series of autopsy-confirmed subjects from the Mayo Clinic Brain Bank (MCBB: 246 AD cases and 223 non-AD Controls) [28,29]. In total, we have 844 cases with LOAD and 1,255 controls without a dementia diagnosis. Mayo Clinic LOAD GWAS data have microarray genotype data profiled for 2,099 European-decent samples that are further imputed to the 1000 Genome Project Phase 3 [37] in our analysis [38]. This cohort only profiles the phenotype of clinical diagnosis of AD.

IGAP GWAS summary statistics of AD

We used the GWAS summary statistics of AD generated from the stage 1 of the International Genomics of Alzheimer’s Project (IGAP) with individuals of all European ancestry [3]. This summary-level GWAS data was generated by meta-analysis (N = ~54K) of four previously-published GWAS datasets consisting of 17,008 Alzheimer’s disease cases and 37,154 controls (The European Alzheimer’s disease Initiative–EADI the Alzheimer Disease Genetics Consortium–ADGC The Cohorts for Heart and Aging Research in Genomic Epidemiology consortium–CHARGE The Genetic and Environmental Risk in AD consortium–GERAD).

Simulation study design

The purpose of this simulation study is to compare the performance of Burden-TWAS and VC-TWAS with variant weights estimated by PrediXcan and DPR methods, as well as validate the VC-TWAS approach of using only summary-level GWAS data. We used the real genotype data from ROS/MAP [40] participants to simulate quantitative gene expression and phenotype traits, where the genotype data were of 2,799 cis-SNPs (with MAF>5% and Hardy Weinberg p-value >10−5) of the arbitrarily chosen gene ABCA7.

Specifically, quantitative gene expression traits are generated by the following equation where denotes the genotype matrix of randomly selected true causal eQTL based on a target proportion of causal eQTL (p) within the test gene, denotes cis-eQTL effect sizes generated from with variance chosen to ensure a target gene expression heritability (), and is the error term generated from .

Phenotype data are generated based on two models to mimic two different genetic architectures of complex traits that may be encountered with human data:

Model I: The genetic effects on the phenotype of interest are completely driven by genetically regulated gene expression (GReX), where SNP effect sizes are of a linear function of their corresponding cis-eQTL effect sizes as assumed by Burden-TWAS methods [8,10,11]. Phenotype data are generated from the following equation where is the gene expression generated from Eq (7) and r = is a scalar chosen to ensure a target phenotype variation proportion due to gene expression ().

Model II: The magnitudes of SNP effect sizes on the phenotype are no longer a linear function of corresponding cis-eQTL effect sizes as in Model I, but rather are derived randomly from a distribution whose variance is a linear function of the squares of such effects. By doing this, the SNP effect sizes on phenotype are a function of cis-eQTL effect sizes but the assumption of a linear relationship between the two (assumed by Burden-TWAS) is relaxed. Phenotype data are generated from the following equation where denotes the genotype matrix of randomly selected true causal SNPs that are also true causal cis-eQTL as in Eq (7), respective SNP effect sizes are generated from β with corresponding cis-eQTL effect size w as used in Eq (7) and r chosen to ensure a phenotype heritability .

We considered scenarios with various proportions of causal cis-eQTL/SNPs p = (0.001, 0.01, 0.1, 0.2) for the test gene, and various combinations of expression heritability () and phenotype variance proportion/heritability () that were chosen to ensure TWAS power falling within the range of (25%, 90%). The values of p and were taken as for simulating phenotypes from Model I, while taken as for simulating phenotypes from Model II.

We used 499 ROS/MAP samples as training data that were also used as training data in our application studies, and randomly selected 1,232 ROS/MAP samples as test data. To show the power performance with respect to test sample size, we considered different test sample sizes (400, 800, 1232) under the scenario with p = 0.2.

We estimated cis-eQTL effect sizes from training data by using PrediXcan and TIGAR DPR methods and then conducted Burden-TWAS and VC-TWAS with individual-level and summary-level test GWAS data. We note that Burden-TWAS using PrediXcan weights is equivalent to the PrediXcan method [10] while Burden-TWAS using DPR weights is equivalent to the TIGAR method [8]. We also compared Burden-TWAS and VC-TWAS with the CoMM method under simulation settings with p = 0.2 (we did not consider additional simulation settings because of the computational demands of CoMM).

For each scenario, we repeated simulations for 1,000 times and obtained the power as the proportion of simulations that had test p-value <2.5×10−6 (genome-wide significance threshold for gene-based test). Additionally, we simulated phenotype under the null hypothesis Y~N(0,1) for 106 times and evaluated type I errors of Burden-TWAS and VC-TWAS, using variant weights derived from PrediXcan and DPR methods.

For each scenario with VC-TWAS, we considered both the original form of the test as well as the alternate computationally-efficient form that considered only the filtered set of variants with cis-eQTL effect size estimates greater than the median effect size value in each simulation (about 50% SNPs).

Results

Simulation results

We compared the performance of VC-TWAS and Burden-TWAS using PrediXcan weights (cis-eQTL effect size estimates by Elastic-Net penalized regression) and DPR weights (cis-eQTL effect size estimates by DPR) under various scenarios. We also evaluated the performance of VC-TWAS and Burden-TWAS using filtered DPR weights as described in Methods. Under the scenario with p = 0.2, we considered different test sample sizes (400, 800, 1232), and compared VC-TWAS with CoMM using both individual-level and summary level test data with 1232 test samples.

First, we compared TWAS power for studying phenotypes simulated from Model I that assumed SNP effect sizes on phenotypes were of a linear function of their corresponding cis-eQTL effect sizes. When p = (0.001, 0.01) with sparse true causal signals (S1A Fig), VC-TWAS had comparable power with Burden-TWAS and both TWAS methods using PrediXcan weights achieved higher power compared to using DPR weights. When p>0.01, we observed that both TWAS methods using DPR weights achieved higher power compared to using PrediXcan weights. These results are consistent with the previous TIGAR paper [8]. This is because DPR method is preferred for modeling quantitative gene expression traits when a gene harbors a considerable proportion of true cis-eQTL with relatively smaller effect sizes, e.g., scenarios with pcausal>0.01 in our simulation studies. As expected, with pcausal>0.01 (S1A Fig, Fig 1A), Burden-TWAS methods outperformed VC-TWAS under Model I which meets the assumptions by Burden-TWAS methods. Across all considered scenarios, all TWAS methods using filtered DPR weights had similar performance as using complete DPR.

Fig 1

Power comparison of Burden-TWAS and VC-TWAS methods under simulation scenarios with 20% true causal eQTL for gene expression (i.e., p = 0.2) in the test gene.

Phenotypes simulated from Model I (A, C, E) and Model II (B, D, F) were considered. Various gene expression heritability and various types of SNP weights were considered, including those derived from PrediXcan method, DPR method, and filtered DPR weights in panel (A, B). Various test sample sizes (400, 800, 1232) were considered in panel (B, D). The VC-TWAS approach of using only summary-level GWAS data were validated in panel (E, F).

Second, we compared TWAS power for studying phenotypes simulated from Model II that assumes variances of SNP effect sizes on phenotype were of a linear function of the squared values of their corresponding cis-eQTL effect sizes. As shown in Figs 1B and S1B, we found that the VC-TWAS obtained higher power than Burden-TWAS across all scenarios. Especially, with p≥0.01, the power by VC-TWAS using DPR weights was double of that using PrediXcan weights on average (88.99% vs. 38.45%), while the power of Burden-TWAS using DPR weights was comparable with using PrediXcan weights except when pcausal = 0.001. When p = 0.001 and (S1B Fig), both TWAS approaches using PrediXcan weights performed better than using DPR weights. Again, across all considered scenarios, TWAS using filtered DPR weights had similar performance to that using complete DPR.

Third, we compared VC-TWAS and Burden-TWAS with different test sample sizes (400, 800, 1232) under the scenarios with (). As shown in Figs 1C and 1D, the power of VC-TWAS and Burden-TWAS increased with respect to the test sample size in both simulation models. Interestingly, for model II, the magnitude of the power difference between VC-TWAS and Burden-TWAS became more pronounced with increasing sample size.

Fourth, to validate the VC-TWAS approach for using only summary-level GWAS data, we compared VC-TWAS with DPR weights by using individual-level and summary-level GWAS data, under the scenarios with p = 0.2. As shown in Figs 1E and 1F, VC-TWAS using summary-level GWAS data (VC-TWAS-SS) performed equivalent as using individual-level data under both simulation models. Additionally, VC-TWAS using filtered DPR weights and summary-level GWAS data (VC-TWAS_f) still had similar performance to that using complete DPR weights.

Fifth, we also compared VC-TWAS with CoMM [20,21] using individual-level and CoMM-SS using summary-level GWAS data under the scenarios with p = 0.2 (S2 Fig). Although CoMM and CoMM-SS outperformed VC-TWAS, CoMM costed up to 600x more computation time and CoMM-SS costed up to 250x more computation time than VC-TWAS for testing a gene with ~5K test SNPs (S1 Table). Take typical genes with ~2K-5K test SNPs as examples, VC-TWAS costed ~20s versus ~38,037s by CoMM using individual-level GWAS data, and VC-TWAS costed ~3s versus ~373s by CoMM-SS using summary-level GWAS data, by a single-thread computation with 4 cores (32GB memory) on a 2.10GHZ CPU (16-node Intel Xeon node). As the number of the SNPs in a gene increases, the difference in computation run times between VC-TWAS and CoMM/CoMM-SS become even more pronounced.

Last but not least, to evaluate type I error of Burden-TWAS and VC-TWAS, we conducted 106 times simulations under the null hypothesis where phenotypes were not associated with genetic data of the test gene. Without loss of generality, we used gene expression data simulated with and generated phenotypes randomly from a N(0,1) distribution. We evaluated type I errors (Table 1) with multiple significant levels (10−2, 10−4, 2.5×10−6), demonstrating that both TWAS approaches had type I errors well controlled across all considered significance levels. We also presented the quantile-quantile (QQ) plots of p-values by all methods in S4 Fig.

Table 1

Type I errors under null simulation studies for Burden-TWAS and VC-TWAS with p = 0.2, = 0.1, using variant weights given by DPR, filtered DPR, and PrediXcan.

Significance Level	Burden-TWAS			VC-TWAS
	DPR	Filtered DPR	PrediXcan	DPR	Filtered DPR	PrediXcan
1.00×10−2	9.82×10−3	9.86×10−3	9.27×10−3	9.43×10−3	9.46×10−3	9.23×10−3
1.00×10−4	8.64×10−5	8.44×10−5	9.95×10−5	8.64×10−5	9.05×10−5	8.24×10−5
2.50×10−6	2.00×10−6	2.00×10−6	2.00×10−6	2.00×10−6	1.00×10−6	6.00×10−6

To summarize, VC-TWAS performed similarly to Burden-TWAS for studying phenotypes simulated from Model I with sparse true causal eQTL, but outperformed Burden-TWAS for studying phenotypes simulated from model II. This is because the genetic architecture assumed under Model I is the one assumed by Burden-TWAS with linear relationship between SNP effect sizes on phenotype and cis-eQTL effect sizes. Whereas, Model II assumes a genetic architecture that is a function of cis-eQTL effect sizes but not explicitly a linear relationship (leading to more general models of effect). Generally, TWAS methods using DPR weights achieved higher power than using PrediXcan weights when pcausal≥0.01, which is consistent with previous studies [8]. Additionally, VC-TWAS using filtered DPR weights achieved similar power as using complete DPR weights, while saving up to 80% of computation time.

Application studies of AD with individual-level GWAS data

We applied VC-TWAS to the individual-level GWAS data of ROS/MAP and Mayo Clinic LOAD cohorts, using SNP weights (i.e., cis-eQTL effect sizes) generated by PrediXcan and filtered DPR methods with 499 ROS/MAP training samples that had both transcriptomic and genetic data profiled [8]. As suggested by previous studies [10,41], our TWAS results included genes with 5-fold cross validation (CV) R2>0.005 for predicting quantitative gene expression traits by either PrediXcan or DPR. We obtained VC-TWAS p-values for 5,710 genes using PrediXcan weights and 12,650 genes using filtered DPR weights. Here, we roughly chose the threshold 10−4 to filter DPR weights in our VC-TWAS such that on average the variance component test considered ~50% SNPs from the test gene region. Specifically, the median number of SNPs considered by VC-TWAS per gene is 2,872 for using filtered DPR weights and 6,632 for using complete DPR weights (S3 Fig).

Leveraging the clinical and postmortem AD data available in ROS/MAP, we were able to apply VC-TWAS to four clinical and pathologic AD phenotypes (S1D Text). We examined final cognitive status diagnosis of 1,436 decedents (AD Dementia (N = 609) versus No Dementia (N = 827)), as well as three postmortem AD phenotypes including: continuous outcomes of β-amyloid load (N = 1,294), PHFtau tangle density (tangles, N = 1,303), and a global AD pathology (N = 1,329). In the VC-TWAS of all four AD phenotypes, we adjusted for covariates of age, smoking status, sex, study group (ROS or MAP), education, and the top three principal components of ancestry.

With Mayo Clinic cohort, we conducted VC-TWAS for AD clinical diagnosis with 844 cases diagnosed with LOAD and 1,255 controls showed no signal of dementia, while adjusting for covariates age, sex, and top three principal components of ancestry. Since only the phenotype of AD clinical diagnosis was profiled by both ROS/MAP and Mayo Clinic cohorts (under different diagnosis criteria) and different sets of covariates were adjusted in VC-TWAS, we conducted meta-analysis with VC-TWAS summary statistics for each study by using Fisher’s method (meta VC-TWAS) [42] to increase power with a larger sample size.

By meta VC-TWAS, we detected 13 significant risk genes with FDR < 0.05 that were located within ~2MB region around the well-known AD risk gene APOE on chromosome 19 (Fig 2A; Table 2). Seven of those significant genes were known risk genes of AD by previous GWAS (CLASRP, TOMM40, MARK4, CLPTM1, CEACAM19, RELB) [43,44] and Burden-TWAS using DPR weights (TRAPPC6A) [8].

Fig 2

Manhattan plots of meta VC-TWAS for AD clinical diagnosis (A) and VC-TWAS of global AD pathology (B) with filtered DPR weights. Genes with FDR < 0.05 are colored in (A), with red for novel risk genes and blue for known AD risk genes. Genes with FDR < 0.05 in meta VC-TWAS of AD clinical diagnosis and p-value < 0.0013 in VC-TWAS of global AD pathology are colored in red in (B).

Table 2

Significant genes for phenotype AD clinical diagnosis by meta VC-TWAS with filtered DPR weights.

Significant genes have FDR < 0.05 by meta-TWAS with ROS/MAP and Mayo Clinic cohorts. AD risk genes identified by previous GWAS are shaded in grey.

Gene name	CHR	Start	End	P-value	FDR
ZNF234a	19	44,645,710	44,664,462	1.11×10−16	1.40×10−12
CLASRP a	19	45,542,298	45,574,214	4.44×10−16	2.81×10−12
TRAPPC6A a	19	45,666,187	45,681,485	3.60×10−14	1.52×10−10
TOMM40 a	19	45,394,477	45,406,935	9.05×10−13	2.86×10−9
MARK4	19	45,754,550	45,808,541	4.62×10−16	1.17×10−8
PPP1R13L	19	45,882,892	45,909,607	1.82×10−10	3.84×10−7
CLPTM1	19	45,457,848	45,496,598	5.71×10−8	1.03×10−4
EML2	19	46,112,660	46,148,726	1.88×10−7	2.97×10−4
FBXO46	19	46,213,887	46,234,151	4.13×10−7	5.80×10−4
CEACAM19 a	19	45,174,724	45,187,631	3.93×10−6	4.68×10−3
GIPR	19	46,171,502	46,185,704	4.07×10−6	4.68×10−3
RELB	19	45,504,695	45,541,452	6.63×10−6	6.99×10−3
ZNF225	19	44,617,548	44,637,255	2.59×10−5	2.51×10−2

a: Genes with significant p-values <0.0013 (Bonferroni correction with respect to 13 genes and 3 phenotypes) for at least one AD pathology phenotype

As clinical AD diagnosis can be associated with AD related pathologies, in further analyses we investigated whether the genes found to be associated with clinical AD diagnosis (Table 2) were also associated with pathologic AD phenotypes. We examined the VC-TWAS p-values of these significant genes with respect to AD pathology phenotypes (β-amyloid, tangles and global AD pathology) (S2 Table; Fig 2B; S5 and S6 Figs). Interestingly, 5 out of these 13 genes had at least one VC-TWAS p-value <0.0013 (Bonferroni correction with respect to 13 genes and 3 AD pathology phenotypes) for one of the three AD pathology phenotypes. Particularly, genes ZNF234 and CLASRP have VC-TWAS p-value <0.0013 for all three AD pathology phenotypes, which are likely to be involved in the biological mechanisms of both β-amyloid and tangles. The other three genes (TRAPPC6A, TOMM40, CEACAM19) have VC-TWAS p-value <0.0013 for β-amyloid and global AD pathology, which are likely to be involved only in the biological mechanism of β-amyloid.

For example, the top significant gene ZNF234 (with FDR = 1.40×10−12) by meta VC-TWAS of AD clinical diagnosis is also the top significant gene (p-value = 2.10×10−4) by VC-TWAS of β-amyloid, the second most significant gene (p-value = 6.39×10−5) by VC-TWAS of global AD pathology, and has p-value = 1.06×10−3 by VC-TWAS of tangles. These results showed that the genetic factor of gene ZNF234 on AD could be potentially mediated through its gene expression, and the expression of this gene could be also potentially involved in the mechanisms of both AD pathology indices of β-amyloid and tangles. Besides AD, gene ZNF234 is also a known risk gene for lipid traits [45]. The genetically regulated gene expression of this gene might also affect lipid traits, thus leading to a pleiotropy phenomenon of AD and lipid traits [46]. Additionally, ZNF234 is known to be involved in the super pathway of gene expression and is annotated with the Gene Ontology term of nucleic acid binding and DNA-binding transcription factor activity [47].

Another significant gene of interest is TOMM40, which has FDR = 2.86×10−9 by meta VC-TWAS for AD dementia and VC-TWAS p-values = (4.44×10−4, 6.95×10−2, 1.91×10−4) for β-amyloid, tangles, and global AD pathology, respectively. These findings suggest that the association of this well-known AD risk gene TOMM40 [48] could be mediated through its gene expression via β-amyloid but not tangles.

For all SNPs considered by meta VC-TWAS of genes ZNF234 and TOMM40 for studying AD clinical diagnosis, we colocalized the meta GWAS results for AD clinical diagnosis with ROS/MAP and Mayo Clinic cohorts and the corresponding DPR weight (i.e., cis-eQTL effect size) magnitude. Interestingly, we found that the VC-TWAS association of these two genes were likely to be driven by SNPs around APOE/TOMM40 loci that also possessed major cis-eQTL effect size magnitudes (S7 Fig).

In addition, our VC-TWAS identified a significant gene HSPBAP1 (FDR = 0.058) for tangles (S5B Fig). As shown by previous studies, mRNA of gene HSPBAP1 was abnormally expressed in the anterior temporal neocortex of patients with intractable epilepsy [49]. Based on our VC-TWAS results, gene HSPBAP1 might not have a significant genetic effect on AD dementia, but might have a significant effect on tangle pathology (p-value = 4.57×10−6) and may account for the increasing recognition of non-cognitive AD phenotypes [50]. This provides support about that gene HSPBAP1 could be involved in the mechanism of brain pathology tangles and other neurological diseases such as intractable epilepsy [49,51].

We also applied VC-TWAS with PrediXcan weights, but no significant genes were identified with FDR <0.05 (S8, S9 and S10 Figs), neither the genes with smallest p-values were proximal to APOE. In contrast, our results by VC-TWAS with filtered DPR weights provided potential biological interpretations for several known AD risk genes via gene expression and for their associations with both clinical and pathologic AD phenotype.

Application studies of AD with summary-level GWAS data

Next, we applied VC-TWAS to the stage1 summary-level GWAS data of AD from IGAP [3], which has a much larger sample size (~54K with 17,008 AD cases and 37,154 controls). We considered filtered cis-eQTL DPR weights as used in the above applications with individual-level GWAS data, as well as the SNP weights of both cis- and trans- eQTL generated by the BGW-TWAS method [30]. LD covariance matrices from ROS/MAP individual-level GWAS data were used for implementing VC-TWAS with summary-level GWAS data.

Using filtered cis-eQTL DPR weights

By using the filtered cis-eQTL DPR weights, we identified 57 significant risk genes with FDR <0.05 by VC-TWAS (Fig 3A; S3 Table and S11A Fig), including the most well-known AD risk genes TOMM40 and APOE, along with 45 genes that were located within ~2MB region around those two risk genes. Other significant genes are distributed on chromosome 6, 8, 11 and 20 (Fig 3A). Among these 57 significant genes, 20 genes were known risk genes by previous GWAS [43,44] (S3 Table). Genes detected by VC-TWAS using the individual-level GWAS data of ROS/MAP and Mayo Clinic cohorts were all replicated (S3 Table). Burden-TWAS using IGAP summary data and filtered cis-eQTL DPR weights (S3 Table) also identified 22 out of 57 of these significant genes.

Fig 3

Manhattan plots of VC-TWAS using IGAP summary data with filtered cis-eQTL DPR weights (A) and TWAS locus zoom plots for the loci on chromosome 11 (B) and chromosome 19 (C). Significant genes with FDR < 0.05 are colored, with red for significant ones only identified by using IGAP summary data and blue for the ones replicating our VC-TWAS findings using the individual-level GWAS data of ROS/MAP and Mayo Clinic cohort. Between gene R2 in (B, C) were calculated with respect to GReX values. The R2 in locus zoom plot denoted by various colors for the dots is the squared correlation of GReX between the most significant gene and other neighborhood genes. The heatmap is based on the squared correlation matrix of GReX.

Since most significant genes are nearby genes on chromosome 11 and 19, we made analogous TWAS locus zoom plots for these two major loci in Fig 3B and 3C, where the between gene R2 were calculated with respect to the predicted GReX values using the individual-level GWAS data of ROS/MAP and Mayo Clinic cohorts. We observed that most significant genes do have highly correlated GReX values, showing that nearby significant TWAS associations are likely to be not independent. For the locus on chromosome 19 (Fig 3C), we can see that the R2 between gene MARK4 (and a few other genes shown in orange color) and the top significant gene CLASRP is greater than 0.4. The significant genes in blue colors are likely to be independent associations from the top significant gene CLASRP. For the locus on chromosome 11 (Fig 3B), significant genes TMEM132A, OSBP, STX3, and TMEM109 are highly correlated with the top significant gene PRPF19, while gene ME3 tend to be another independent association.

Using cis- and trans- eQTL weights generated by BGW-TWAS

To provide a complementary list of significant genes by considering both cis- and trans- eQTL, we conducted VC-TWAS with the IGAP summary data using the cis- and trans- eQTL weights generated by the BGW-TWAS method [30]. We detected total 71 significant genes with FDR < 0.05 (Fig 4; S4 Table; S11B Fig), among which 6 genes were identified by GWAS [44,52-55] and 22 genes were shown to be related with AD or other neurological diseases by previous studies and (Table 3).

Fig 4

Manhattan plot of VC-TWAS using IGAP summary data with BGW cis- and trans- eQTL weights.

Significant genes with FDR < 0.05 are colored.

Table 3

Significant genes identified by VC-TWAS using IGAP summary statistics data with BGW cis- and trans- eQTL weights, which are either known GWAS loci or shown to be related with AD or other neurological diseases by previous studies.

AD risk genes identified by previous GWAS are shaded in grey.

Gene Name	CHROM	Start	End	P-value	FDR
ARHGEF2a	1	155,916,644	155,966,129	1.75×10−16	2.46×10−13
GAS5 a	1	173,833,037	173,838,020	1.33×10−22	3.11×10−19
NCOA1 a	2	24,714,782	24,993,571	2.42×10−25	8.51×10−22
CENPO a	2	25,016,004	25,045,245	4.07×10−7	1.85×10−4
SLC4A10 a	2	162,280,842	162,841,792	8.80×10−6	2.58×10−3
MAN2A1	5	109,025,066	109,205,326	6.30×10−8	3.70×10−5
ZMAT2 a	5	140,079,746	140,086,266	1.83×10−6	6.43×10−4
HLA-DRB5 a	6	32,485,119	32,498,064	7.81×10−7	3.23×10−4
AHR a	7	17,338,245	17,385,776	2.11×10−5	5.62×10−3
JAZF1	7	27,870,191	28,220,362	7.49×10−8	4.22×10−5
ACHE a	7	100,487,614	100,494,594	5.93×10−5	1.39×10−2
DLGAP2 a	8	1,449,530	1,656,642	1.18×10−5	3.38×10−3
RFX3 a	9	3,218,296	3,526,004	4.69×10−5	1.16×10−2
SLC25A28 a	10	101,370,281	101,380,535	1.55×10−7	7.79×10−5
PDCD11 a	10	105,156,404	105,206,049	2.52×10−8	1.69×10−5
ROBO4 a	11	124,753,586	124,768,396	1.13×10−6	4.30×10−4
USP30 a	12	109,460,893	109,525,831	1.50×10−10	1.41×10−7
PPP1R3E a	14	23,765,111	23,772,057	2.11×10−5	5.62×10−3
MTFMT a	15	65,294,844	65,321,977	3.10×10−5	8.09×10−3
PARP6	15	72,533,521	72,565,340	6.77×10−6	2.07×10−3
GRIN2A a	16	9,852,375	10,276,611	2.22×10−8	1.56×10−5
SETD6 a	16	58,549,382	58,554,431	3.44×10−5	8.81×10−3
NUP88 a	17	5,264,257	5,323,480	2.65×10−7	1.24×10−4
NOS2 a	17	26,083,791	26,127,555	9.55×10−5	2.07×10−2
CEACAM19	19	45,174,723	45,187,631	1.46×10−5	4.03×10−3
APOC1 a	19	45,417,920	45,422,606	2.80×10−6	9.63×10−4

a: Genes shown to be related with AD or other neurological diseases by previous studies

For example, gene ARHGEF2 on chromosome 1(with FDR = 2.46×10−13) is shown to be interact with all four members of the MARK family including MARK4 on chromosome 19 which has an emerging tole in the phosphorylation of MAPT/TAU in Alzheimer’s disease and identified by VC-TWAS using cis-eQTL DPR weights [56]. Gene GAS5 (with FDR = 3.11×10−19) was shown to have a novel role in microglial polarization and the pathogenesis of demyelinating diseases which suggested the potential therapeutic benefit of targeting GAS5 for the treatment of neurological disorders including AD [57]. Gene NCOA1 (with FDR 8.51×10−22) was shown to have pivotal roles in memory and learning [58]. SNPs in Gene PDCD11 (with FDR = 1.69×10−5) have been showed to be significantly associated with AD [59]. Gene APOC1 was detected by several GWASs [43,44], even when conditioning on APOE gene [60]. The presence of gene ACHE (with FDR = 1.39×10−2) is a common feature in AD brain [61]. We also detected genes NUP88 [62], ROBO4 [63], DLGAP2 [64], AHR [65], PPP1R3E [66], RFX3 [67] and NOS2 [68] that were shown to possess biological link with AD. Additional identified genes USP30 [69], GRIN2A [70], SLC25A28 [71], HLA-DRB5 [72], ZMAT2 [73], SLC4A10 [74], MTFMT [75] and SETD6 [76] were shown to be related with neurological diseases by previous studies.

Overall, we show that VC-TWAS using BGW weights provide a complementary list of risk genes to those previously identified by considering additional trans- eQTL information.

Discussion

In this paper, we propose a novel variance-component TWAS (VC-TWAS) method that leverages eQTL effect sizes from reference transcriptome but does not assume a linear relationship between SNP effect sizes on phenotypes and cis-eQTL effect sizes. VC-TWAS is applicable to both quantitative and dichotomous outcomes and can further handle both individual-level and summary-level GWAS data. By implementing this VC-TWAS with cis-eQTL effect sizes estimated by DPR method [8,13], we created a powerful test statistic that had good performance in simulation studies and obtained biologically meaningful TWAS results for both clinical and pathologic AD phenotypes. In particular, with the individual-level GWAS data of ~3.4K samples, we detected 13 TWAS genes for AD dementia, including the well-known GWAS risk gene TOMM40 and previously identified TWAS gene TRAPPC6A [8]. Moreover, 6 out of these 13 genes were identified by previous GWAS [43]. The pleiotropy effects of 5 of these genes with respect to AD dementia and indices of AD pathology demonstrated the possible biological mechanisms linking AD risk genes via β-amyloid and tangles with AD dementia.

By applying VC-TWAS with summary-level GWAS data of AD with a much larger sample size, we not only validated our findings with the individual-level GWAS data but also detected additional novel risk genes. In particular, by applying VC-TWAS with both cis- and trans- eQTL effect sizes estimated by the BGW-TWAS method, we identified a list of significant AD risk genes that complement risk genes identified by considering only cis-eQTL information.

To help users conduct our VC-TWAS method conveniently and efficiently, we added this method into our previously developed tool––Transcriptome Integrated Genetic Association Resource (TIGAR) [8]. The user has the option of estimating cis-eQTL effect sizes by either the PrediXcan method (i.e., Elastic-Net) [12] or nonparametric Bayesian DPR method [13], or using eQTL weights generated by previous studies. The VC-TWAS function works for using individual-level GWAS data to study continuous or dichotomous phenotypes, as well as using summary-level GWAS data.

Since the variance component test statistic used by VC-TWAS involves calculating and performing an eigen-decomposition of a genotypic kernel matrix, efficient computation is required (even when filtering variants to include only those variants with relatively larger magnitude of eQTL estimates) for obtaining the corresponding p-values in practice. Our TIGAR tool implements multi-threaded computation to take advantage of high-performance cloud computing clusters and enable practical computation for testing genome-wide genes with respect to reference transcriptome data of multiple tissue types. A typical VC-TWAS with ~20K genes can be accomplished within ~80 hours by using a single thread with 1 CPU core and 32 GB memory.

Of course, current TWAS methods including our VC-TWAS still have their limitations. Because of genetic and transcriptomic heterogeneities across different ethnicities, one may have difficulty in translating eQTL effect size estimates across cohorts with different ethnicities [77]. That being said, VC-TWAS is likely more robust to this phenomenon than Burden-TWAS given the former relaxes the assumption of the latter of a linear relationship between the SNP effect size and eQTL effect in the transcriptome. Nevertheless, reference panels with diverse ethnicities and multiple tissue types are certainly needed to enhance TWAS to study complex diseases across different ancestral groups. We should note though that both ROS/MAP and Mayo Clinic cohorts that we considered in this work have similar (European) origins.

Compared to the alternative CoMM and PMR-Egger methods that jointly model both reference and test data, VC-TWAS might be less powerful when the eQTL effect sizes are homogeneous in both reference and test cohorts. However, CoMM and PMR-Egger methods are derived for quantitative phenotypes and implement likelihood ratio tests under a maximum likelihood reference framework that are computationally expensive when thousands of cis-SNPs need to be tested per gene. Although one may apply CoMM and PMR-Egger methods to dichotomous phenotypes by taking cases as 1’s and controls as 0’s, this may lead to inflated type I errors when population stratification leads to violation of the constant-residual variance assumption. The recently published generalized linear mixed model association test (GMMAT) method paper has shown this by both simulation and real studies and suggested that a logistic mixed model would be more appropriate for analyzing dichotomous traits [78]. Thus, VC-TWAS would be preferred because of computational efficiency when imputed or whole genome sequencing genotype data are considered for the reference panel, or when dichotomous traits are analyzed. We provide a summary table for pros and cons of existing popular TWAS methods in S5 Table.

Nevertheless, our novel technique allows a more flexible framework to account for the unknown genetic architectures underlying the relationships between SNPs and the phenotype of interest with efficient computation. Using simulated and real AD-related data, we show our method VC-TWAS provides the public a useful tool for illustrating the genetic etiology of complex diseases by providing a list of risk genes whose effects on phenotypes might be mediated through transcriptomes.

Web resources

VC-TWAS, https://github.com/yanglab-emory/VC_TWAS

TIGAR, https://github.com/yanglab-emory/TIGAR

PrediXcan, https://github.com/hakyim/PrediXcan

RADC Research Resource Sharing Hub, http://www.radc.rush.edu/

ROS/MAP data, https://www.synapse.org/#!Synapse:syn3219045

MayoClinicLOAD data, https://www.synapse.org/#!Synapse:syn2910256

GWAS catalog, https://www.ebi.ac.uk/gwas/

BGW weights of brain frontal cortex tissue, https://www.synapse.org/#!Synapse:syn22316792

Details about PrediXcan’s and TIGAR’s approach of estimating cis-eQTL effect sizes (A, B), VC-TWAS approach with summary-level GWAS data (C), and ROS/MAP data (D).

(DOCX)

Click here for additional data file.

TWAS power comparison for VC-TWAS and Burden-TWAS with phenotypes simulated from Model I (A) and Model II (B). Various types of SNP weights were considered, including those derived from PrediXcan method, DPR method, and filtered DPR weights. In Model I, the combinations of causal probability and phenotype heritability are . In Model II, the combinations of causal probability and phenotype heritability are .

(PDF)

Click here for additional data file.

TWAS power comparison for VC-TWAS and CoMM with phenotypes simulated from Model I (A) and Model II (B) using individual-level and summary-level data under the scenarios with p = 0.2.

(PDF)

Click here for additional data file.

Box plot of the number of test SNPs considered by VC-TWAS of all genome-wide genes in the application studies of AD, with complete DPR weights and filtered DPR weights derived from the ROS/MAP training data.

(PDF)

Click here for additional data file.

Q-Q plots for VC-TWAS and Burden-TWAS with DPR weights, filtered DPR weights, and PrediXcan weights under null hypothesis, where quantitative gene expression traits were generated with p = 0.2 and = 0.1.

(PDF)

Click here for additional data file.

Manhattan plots of VC-TWAS results with filtered DPR weights for studying quantitative AD pathology of β-Amyloid (A) and tangles (B). Genes with FDA < 0.05 by meta VC-TWAS for studying AD clinical diagnosis are colored in red in (A) and top significant gene for studying tangles phenotype with FDR = 0.058 is colored in red in (B).

(PDF)

Click here for additional data file.

Q-Q plots of VC-TWAS results with filtered DPR weights for studying β-amyloid, tangles, and global AD pathology with ROS/MAP cohort, as well as meta VC-TWAS results with filtered DPR weights for studying AD clinical diagnosis with ROS/MAP and Mayo Clinic cohorts.

(PDF)

Click here for additional data file.

Locus zoom plots of GWAS results and the magnitude (i.e., absolute value) of cis-eQTL effect size estimates by DPR method for SNPs that were considered by VC-TWAS of genes TOMM40 (A, B) and ZNF2334 (C, D). Filtered test SNPs with the cis-eQTL effect size magnitude > 10−4 were plotted here. SNPs with GWAS p-value <5×10−8 were colored in red in (B,D), top significant SNPs by GWAS in (A,C) were shown as the blue triangle in (B,D).

(PDF)

Click here for additional data file.

Manhattan plots of VC-TWAS results with PrediXcan weights for studying AD clinical diagnosis (A) and global AD pathology (B).

(PDF)

Click here for additional data file.

Manhattan plots of VC-TWAS results with PrediXcan weights for studying quantitative AD pathology of β-Amyloid (A) and tangles (B).

(PDF)

Click here for additional data file.

Q-Q plots of VC-TWAS results with PrediXcan weights for studying β-amyloid, tangles, and global AD pathology with ROS/MAO cohort, as well as meta VC-TWAS results with PrediXcan weights for studying AD clinical diagnosis with ROS/MAP and Mayo Clinic cohorts.

(PDF)

Click here for additional data file.

Q-Q plots of VC-TWAS results with cis-eQTL DPR filtered weights and BGW weights on IGAP GWAS summary statistics.

(PDF)

Click here for additional data file.

Average computation time per gene (in the unit of second) by all methods with individual and summary-level GWAS data by a single thread with 4 cores (32GB memory) from a 2.10GHZ CPU 16-core Intel Xeon computation node, for example genes that have test SNP numbers in respective range.

(DOCX)

Click here for additional data file.

Genes with VC-TWAS p-value <0.0013 with respect to at least one AD pathology phenotype and FDR <0.05 by meta VC-TWAS of AD clinical diagnosis.

(DOCX)

Click here for additional data file.

Significant genes identified by VC-TWAS using IGAP summary statistics with filtered cis-eQTL DPR weights.

Significant genes were identified with FDR < 0.05. AD risk genes identified by previous GWAS are shaded in grey.

(DOCX)

Click here for additional data file.

Novel significant genes identified by VC-TWAS using summary statistics with BGW weights on IGAP summary statistics.

(DOCX)

Click here for additional data file.

Pros and cons of existing popular TWAS methods.

(DOCX)

Click here for additional data file.

15 Sep 2020

Dear Dr Yang,

Thank you very much for submitting your Research Article entitled 'Novel Variance-Component TWAS method for studying complex human diseases with applications to Alzheimer's dementia' to PLOS Genetics. Your manuscript was fully evaluated at the editorial level and by independent peer reviewers. The reviewers appreciated the attention to an important problem, but raised some substantial concerns about the current manuscript. Based on the reviews, we will not be able to accept this version of the manuscript, but we would be willing to review again a much-revised version. We cannot, of course, promise publication at that time.

The manuscript proposed a variance-component-based TWAS method and applied it to study Alzheimer's disease.  Considering that there are many recent TWAS methods proposed in the literature, one theme of the reviewers is that additional comparisons with more recent methods are needed. In addition, one reviewer raised some important points regarding the application to ROSMAP and Mayo Clinic data which must be addressed in the revision.

Should you decide to revise the manuscript for further consideration here, your revisions should address the specific points made by each reviewer. We will also require a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript.

If you decide to revise the manuscript for further consideration at PLOS Genetics, please aim to resubmit within the next 60 days, unless it will take extra time to address the concerns of the reviewers, in which case we would appreciate an expected resubmission date by email to plosgenetics@plos.org.

If present, accompanying reviewer attachments are included with this email; please notify the journal office if any appear to be missing. They will also be available for download from the link below. You can use this link to log into the system when you are ready to submit a revised version, having first consulted our Submission Checklist.

Please be aware that our data availability policy requires that all numerical data underlying graphs or summary statistics are included with the submission, and you will need to provide this upon resubmission if not already present. In addition, we do not permit the inclusion of phrases such as "data not shown" or "unpublished results" in manuscripts. All points should be backed up by data provided with the submission.

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool.  PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

PLOS has incorporated Similarity Check, powered by iThenticate, into its journal-wide submission system in order to screen submitted content for originality before publication. Each PLOS journal undertakes screening on a proportion of submitted articles. You will be contacted if needed following the screening process.

To resubmit, use the link below and 'Revise Submission' in the 'Submissions Needing Revision' folder.

[LINK]

We are sorry that we cannot be more positive about your manuscript at this stage. Please do not hesitate to contact us if you have any concerns or questions.

Yours sincerely,

Lin Chen, Ph.D.

Guest Editor

PLOS Genetics

David Balding

Section Editor: Methods

PLOS Genetics

The manuscript proposed a variance-component-based TWAS method and have applied it to study Alzheimer's disease. Three reviewers have read the work and made some suggestions. Considering that there are many recent TWAS methods proposed in the literature, additional comparisons with more recent methods would be helpful. Moreover, one reviewer raised some valid points regarding the application to ROSMAP and Mayo Clinic data. This should be clarified and addressed in the revision.

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Please see the attachment for comments

Reviewer #2: This paper proposes a novel Variance-Component TWAS procedure (VC-TWAS) that assumes the effects of cis-eQTL SNPs on phenotype are random (with variance proportional to corresponding cis-eQTL effects in reference dataset) rather than fixed. They further applied VC-TWAS using cis-eQTL effect sizes estimated by a nonparametric Bayesian method to study Alzheimer’s dementia (AD) related phenotypes and detected 13 genes significantly associated with AD, including 6 known GWAS risk loci such as TOMM40. It is an interesting problem that authors propose variance component tests to perform TWAS. Here are my major comments.

1. First, I think authors should compare other alternative methods in TWAS as there are several methods that may outperform standard TWAS, e.g., [1] and [2]. Authors need to numerically compare VC-TWAS these methods as well. Particularly, [1] connects CoMM with other SKAT and PrediXcan, showing how the improvement comes from (Section 2.4 in [1]). Since the test for equation (8) is closely related to equation (14) in [1], authors should cite the relevant reference for their works.

2. In the simulation studies, what are the sample sizes for equation (9)? As sample sizes for eqtl studies increase, the difference between PrediXcan and other refined ones becomes similar. In addition, Table 1 that summarizes the type I error is not enough especially for genome-wide scan. qq-plots for null simulations are more appropriate for type I error control.

3. In abstract, it is too strong to state that VC-TWAS gain > 2X power. This can be only done by simulation studies with some favorable settings. I am wondering why authors did not apply alternative methods to the real data analysis. The comparison in the real data analysis will be more interesting.

[1] Yang C, Wan X, Lin X, et al. CoMM: a collaborative mixed model to dissecting genetic contributions to complex traits by leveraging regulatory information[J]. Bioinformatics, 2019, 35(10): 1644-1652.

[2] Nagpal S, Meng X, Epstein M P, et al. TIGAR: an improved Bayesian tool for transcriptomic data imputation enhances gene mapping of complex traits[J]. The American Journal of Human Genetics, 2019, 105(2): 258-266.

Reviewer #3: Tang et al. proposed a TWAS method following the idea of the gene-based variance component test SKAT with cis-variant weights trained with Dirichlet process regression (DPR). The idea of combining TWAS and SKAT is novel, and the applications on Alzheimer’s disease data are interesting. The manuscript is also clearly written. Below I list some comments that may further improve the paper.

Major comments:

“As DPR method produces non-zero cis-eQTL effect size estimates for almost all SNPs within a test gene region (with most cis-eQTL effect size estimates being close to zero [8]), we explored an alternate VC-TWAS that considered a reduced set of SNPs by filtering out those with cis-eQTL effect size estimates smaller than the median cis-eQTL effect size estimate.”

1) Is the median cis-eQTL effect size estimate for each gene or all genes?

2) Other than using this arbitrary threshold for filtering, a method BGW-TWAS developed earlier by the authors’ group adopts a spike-and-slab prior that is commonly used in Bayesian variable selection. I am wondering what the performance of DPR as compared to the spike-and-slab prior.

3) Additionally, BGW-TWAS considers both cis and trans- eQTLs and has trained the same datasets as this manuscript. It would be interesting to use the BGW-TWAS trained cis- and trans- eQTL weights in this VC-TWAS approach. This will help compare: i) DPR vs. the spike-and-slab prior, and ii) cis- only vs. both cis- and trans- eQTLs.

4) Given that the main contribution of this manuscript is to use different SNP weights, I am curious about the real data results using the original SKAT to see the gain by using this new weighting via TWAS.

“Mayo Clinic LOAD GWAS data contain samples from two clinical AD Case-Control series:

Mayo Clinic Jacksonville (MCJ: 353 AD cases and 331 Controls), Mayo Clinic Rochester (MCR: 291 AD cases and 787 Controls) and a neuropathological series of autopsy-confirmed subjects from the Mayo Clinic Brain Bank (MCBB: 298 AD cases and 223 non-AD Controls) [24, 25]. In total, we have 844 cases with LOAD and 1,255 controls without a dementia diagnosis.”

1) This part was from the data description on the Synapse website, but the sum of the sample size of the three cohorts does not equal to the total number listed above. Specifically, the sample size of MCR and AD cases of MCBB does not match with those reported in Allen et al. Scientific Data 2016.

2) In the real data analysis of Mayo Clinic data, the authors probably need to adjust for the cohort variable. MCJ and MCR samples are antemortem, but MCBB samples are postmortem.

It is not clear if the outcome/AD definition is the same across ROSMAP data and the three cohorts of Mayo Clinic data. Below is an excerpt on the AD definition of Mayo data copied from the Synapse website.

“All subjects from the clinical series (MCJ and MCR) were diagnosed by a Mayo Clinic neurologist; all control subjects had a Clinical Dementia Rating score of zero at the most recent time of testing; all LOAD patients had a diagnosis of probable or possible AD according to the NINCDS-ADRDA criteria. Subjects selected from the Mayo Clinic Brain Bank underwent neuropathologic evaluation by Dr. Dennis Dickson. All ADs had definite diagnosis according to the NINCDS-ADRDA criteria and had Braak scores of ≥4.0. All non–AD Controls had Braak scores of ≤2.5; many had brain pathology unrelated to AD and, in particular, many of the non-AD samples used in the brain expression data were selected based on the presence of primary tauopathy (progressive supranuclear palsy).”

1) Were all AD patients in both ROSMAP and Mayo LOAD GWAS data diagnosed by a neurologist with the same criterion?

2) Even if 1) is true, the diagnosis of subjects in Mayo data has been assisted with NINCDS-ADRDA criteria and Braak score. For instance, “all LOAD patients had a diagnosis of probable or possible AD according to the NINCDS-ADRDA criteria.” Is this true for ROSMAP data? “All MCBB ADs had definite diagnosis according to the NINCDS-ADRDA criteria and had Braak scores of ≥4.0, all non–AD Controls had Braak scores of ≤2.5” is this true for all other cohorts?

I know it may be hard to check this, but making sure each dataset has the same AD definition and matched other clinical variables would facilitate the interpretation of the meta-analysis.

3) The non-AD controls in MCBB have other pathologies. Is this true for other cohorts? Would this cofound the comparison as well?

Before meta-analysis is done, are there any overlapping signals between the two data cohorts: ROSMAP and Mayo LOAD GWAS? This would help check the reproducibility of results.

“11 of these genes were shown to possess significant pleiotropic effects on brain pathology phenotypes.”

This significance is based on a p-value threshold of 0.05 (Table 3, Figure 2B) and has not been adjusted for multiple testing of 11 genes and three phenotypes/outcomes.

Minor comments:

Typo: Page 15, first line, CPU *munities*;

Page 19, line 457, “using individual-level DATA for”? “Data” was omitted.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Genetics data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: None

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

Submitted filename: Review.pdf

Click here for additional data file.

Submitted filename: review.docx

Click here for additional data file.

18 Jan 2021

Submitted filename: ResponseToReviewers.docx

Click here for additional data file.

11 Feb 2021

Dear Dr Yang,

Thank you very much for submitting a revised version of your Research Article entitled 'Novel Variance-Component TWAS method for studying complex human diseases with applications to Alzheimer's dementia' to PLOS Genetics.

Although two reviewers are now satisfied, Reviewer 2 continues to have concerns, primarily relating to the accuracy and fairness of descriptions of and comparisons with the CoMM software.  We ask you to further revise your manuscript to address these concerns, which can involve either performing new comparisons or a further explanation of the reasons for the discrepancies that the reviewer has highlighted.  We encourage you also to make software for your analyses available as the reviewer requests.

Should you decide to revise the manuscript for further consideration here, you should address all the points made by the reviewer either through appropriate revisions or rebuttal in a response letter. We will also require a description of the changes you have made in the manuscript.

We hope that you will revise the manuscript for further consideration at PLOS Genetics, please aim to resubmit within the next 60 days, unless it will take extra time to address the concerns of the reviewers, in which case we would appreciate an expected resubmission date by email to plosgenetics@plos.org.

If present, accompanying reviewer attachments are included with this email; please notify the journal office if any appear to be missing. They will also be available for download from the link below. You can use this link to log into the system when you are ready to submit a revised version, having first consulted our Submission Checklist.

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see our guidelines.

Please be aware that our data availability policy requires that all numerical data underlying graphs or summary statistics are included with the submission, and you will need to provide this upon resubmission if not already present. In addition, we do not permit the inclusion of phrases such as "data not shown" or "unpublished results" in manuscripts. All points should be backed up by data provided with the submission.

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool.  PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

PLOS has incorporated Similarity Check, powered by iThenticate, into its journal-wide submission system in order to screen submitted content for originality before publication. Each PLOS journal undertakes screening on a proportion of submitted articles. You will be contacted if needed following the screening process.

To resubmit, use the link below and 'Revise Submission' in the 'Submissions Needing Revision' folder.

[LINK]

We are sorry that we cannot be more positive about your manuscript at this stage. Please do not hesitate to contact us if you have any concerns or questions.

Yours sincerely,

Lin Chen, Ph.D.

Guest Editor

PLOS Genetics

David Balding

Section Editor: Methods

PLOS Genetics

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Thanks for the revisions, my questions and concerns are well addressed.

Reviewer #2: Thank authors for detailed and comprehensive replies. I still have some comments.

1. In practice, many TWAS methods consider cis-SNPs to be either upstream or downstream 500kb with the range of a gene and the number of cis-SNPs ranges from tens of a few hundreds. In Table S1, the number of cis-SNPs considered ranges from 1000 to >5000. First, authors need to justify that increasing cis-SNP number can enhance power. Otherwise, I don’t see any points to consider such a large number of cis-SNPs to compare the computational efficiency. At least, authors need to show that using such a large of cis-SNPs in VC-TWAS can outperform CoMM/CoMM-S2 that uses a few hundred cis-SNPs. I guess authors are using imputed SNPs in the analysis. If so, authors need to justify that using imputed SNPs can improve power in TWAS over the genotyped SNPs. Second, to further account for horizontal pleiotropy, there are couple of transcriptome-wide Mendelian randomization methods proposed, e.g., PMR-Egger [1] that primary assumes a burden assumption. It will be also interesting to compare performance of VC-TWAS with PMR-Egger as well. In PMR-Egger, the authors compare the computational efficiency in Table 2 [1], but Table S1 is significantly different from that one. This point should be addressed seriously by the authors.

Overall, without justification for a large of cis-SNPs as Table S1, it is pretty aggressive to state as in the paragraph between lines 358 and 365 (as well between lines 109 and 111), which are absolutely misleading for readers as CoMM was designed to consider only a few hundreds genotyped cis-SNPs, especially such conclusions are based on the incorrect Table S1. Moreover, readers need to know whether they should use imputed cis-SNP or genotyped cis-SNP. If the power from genotyped cis-SNP is comparable to that of imputed cis-SNP, there is no reason to make comparisons in Table S1. I suggest the authors revise the parts I mentioned here.

2. According to CoMM, “The analysis for an individual trait in the NFBC1966 dataset can be done around 2 min and around 18 min for an individual disease in the GERA dataset on a Linux platform with 2.6GHz intel Xeon CPU E5-2690 with 30 720 KB cache and 96 GB RAM (only 7.8 GB RAM used) on 24 cores.” In CoMM, the authors assume that the number of cis-SNPs is less than sample size which makes sense as the median of the number of cis-SNPs (500kb on each side) in GWAS is around 200. For this, CoMM solves the problem using the X’X, which X is genotype and the computation is in order less than O(p^3). If the number of cis-SNPs is larger than the sample size, it is optimal to solve using XX’ and the computation is at most O(n^3). This is IMPORTANT based on your assumptions for the number cis-SNPs. In Table S1, we should expect at most 73*5^3 = 9125 seconds for SNP number > 5000 as I explained above. But the reported time is 10 times that of expected. Again, the results in Table S1 here is significantly different from Table 2 reported in PMR-Egger [1]. Thus, it must be somethings wrong here and authors must correct any mistakes for this table. As PMR-Egger is an extension considering a fixed intercept (burden test), authors may compare it as well. For a fair comparison, authors need to state clearly for the assumptions that all methods made. It will be helpful to have a table of summary for pros, cons, and any assumptions among different methods.

3. Though CoMM was developed using linear models, it can still be applied to binary outcome. The differences between using linear and logistic regression would not be that much as long as the case-control ratio is too extreme which is the case in this study. At least, it does not hurt to compare. Moreover, the sentence in introduction between lines 108 and 109 should be rewritten. It is more informative to state the assumption of VC-TWAS as the binary outcome and CoMM as the continuous outcome. It is not necessary that a linear model cannot be applied to the studies with binary outcomes. Please revise the corresponding statement in the manuscript.

4. To promote reproducibility, all codes for results (tables and figures) in this manuscript should be provided, e.g. in github.

5. Provided the parallel implementation, it is computationally feasible using CoMM for large sample size as GERA (PMR-Egger even used UKB data) so it would be interesting to compare CoMM/CoMM-S2 as well as PMR-Egger with VC-TWAS in real data analysis. Even without parallelization, it would at most take 18/60*24=7.2 hours to complete analysis for a trait in GERA that is much larger than the sample used in this manuscript.

[1] Yuan, Z., Zhu, H., Zeng, P., Yang, S., Sun, S., Yang, C., ... & Zhou, X. (2020). Testing and controlling for horizontal pleiotropy with probabilistic Mendelian randomization in transcriptome-wide association studies. Nature communications, 11(1), 1-14.

Reviewer #3: none

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Genetics data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: No: To promote reproducibility, please provides all codes for Figures and Tables.

Reviewer #3: None

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

Submitted filename: review2.docx

Click here for additional data file.

4 Mar 2021

Submitted filename: Response_R2_03_03_2021.docx

Click here for additional data file.

11 Mar 2021

Dear Dr Yang,

Thank you very much for submitting your Research Article entitled 'Novel Variance-Component TWAS method for studying complex human diseases with applications to Alzheimer's dementia' to PLOS Genetics.

Reviewer 2 has only minor suggestions that we believe can be easily addressed.  In particular we think that the table of pros and cons of different methods would be useful. Please address the reviewer comments, we don't plan to send it back to reviewers again and so we hope to move quickly to acceptance once we receive your revision.

Please also:

1) Provide a list of your responses to the review comments and a description of the changes you have made in the manuscript.

2) Upload a Striking Image with a corresponding caption to accompany your manuscript if one is available (either a new image or an existing one from within your manuscript). If this image is judged to be suitable, it may be featured on our website. Images should ideally be high resolution, eye-catching, single panel square images. For examples, please browse our archive. If your image is from someone other than yourself, please ensure that the artist has read and agreed to the terms and conditions of the Creative Commons Attribution License. Note: we cannot publish copyrighted images.

We hope to receive your revised manuscript very soon.  If you anticipate any delay beyond 30 days we ask you to let us know the expected resubmission date by email to plosgenetics@plos.org.

Please notify the journal office if the reviewer attachment is missing. They will also be available for download from the link below. You can use this link to log into the system when you are ready to submit a revised version, having first consulted our Submission Checklist.

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Please be aware that our data availability policy requires that all numerical data underlying graphs or summary statistics are included with the submission, and you will need to provide this upon resubmission if not already present. In addition, we do not permit the inclusion of phrases such as "data not shown" or "unpublished results" in manuscripts. All points should be backed up by data provided with the submission.

Please review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript.

PLOS has incorporated Similarity Check, powered by iThenticate, into its journal-wide submission system in order to screen submitted content for originality before publication. Each PLOS journal undertakes screening on a proportion of submitted articles. You will be contacted if needed following the screening process.

To resubmit, you will need to go to the link below and 'Revise Submission' in the 'Submissions Needing Revision' folder.

[LINK]

Please let us know if you have any questions while making these revisions.

Yours sincerely,

Lin Chen, Ph.D.

Guest Editor

PLOS Genetics

David Balding

Section Editor: Methods

PLOS Genetics

Reviewer's Responses to Questions

Comments to the Authors:

Reviewer #2: Attached

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Genetics data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #2: No

Submitted filename: review3.docx

Click here for additional data file.

12 Mar 2021

Submitted filename: Response_R3_03_11_2021.docx

Click here for additional data file.

15 Mar 2021

Dear Dr Yang,

We are pleased to inform you that your manuscript entitled "Novel Variance-Component TWAS method for studying complex human diseases with applications to Alzheimer's dementia" has been editorially accepted for publication in PLOS Genetics. Congratulations!

Before your submission can be formally accepted and sent to production you will need to complete our formatting changes, which you will receive in a follow up email. Please be aware that it may take several days for you to receive this email; during this time no action is required by you. Please note: the accept date on your published article will reflect the date of this provisional acceptance, but your manuscript will not be scheduled for publication until the required changes have been made.

Once your paper is formally accepted, an uncorrected proof of your manuscript will be published online ahead of the final version, unless you’ve already opted out via the online submission form. If, for any reason, you do not want an earlier version of your manuscript published online or are unsure if you have already indicated as such, please let the journal staff know immediately at plosgenetics@plos.org.

In the meantime, please log into Editorial Manager at https://www.editorialmanager.com/pgenetics/, click the "Update My Information" link at the top of the page, and update your user information to ensure an efficient production and billing process. Note that PLOS requires an ORCID iD for all corresponding authors. Therefore, please ensure that you have an ORCID iD and that it is validated in Editorial Manager. To do this, go to ‘Update my Information’ (in the upper left-hand corner of the main menu), and click on the Fetch/Validate link next to the ORCID field.  This will take you to the ORCID site and allow you to create a new iD or authenticate a pre-existing iD in Editorial Manager.

If you have a press-related query, or would like to know about making your underlying data available (as you will be aware, this is required for publication), please see the end of this email. If your institution or institutions have a press office, please notify them about your upcoming article at this point, to enable them to help maximise its impact. Inform journal staff as soon as possible if you are preparing a press release for your article and need a publication date.

Thank you again for supporting open-access publishing; we are looking forward to publishing your work in PLOS Genetics!

Yours sincerely,

Lin Chen, Ph.D.

Guest Editor

PLOS Genetics

David Balding

Section Editor: Methods

PLOS Genetics

www.plosgenetics.org

Twitter: @PLOSGenetics

----------------------------------------------------

Comments from the reviewers (if applicable):

----------------------------------------------------

Data Deposition

If you have submitted a Research Article or Front Matter that has associated data that are not suitable for deposition in a subject-specific public repository (such as GenBank or ArrayExpress), one way to make that data available is to deposit it in the Dryad Digital Repository. As you may recall, we ask all authors to agree to make data available; this is one way to achieve that. A full list of recommended repositories can be found on our website.

The following link will take you to the Dryad record for your article, so you won't have to re‐enter its bibliographic information, and can upload your files directly:

http://datadryad.org/submit?journalID=pgenetics&manu=PGENETICS-D-20-01014R3

More information about depositing data in Dryad is available at http://www.datadryad.org/depositing. If you experience any difficulties in submitting your data, please contact help@datadryad.org for support.

Additionally, please be aware that our data availability policy requires that all numerical data underlying display items are included with the submission, and you will need to provide this before we can formally accept your manuscript, if not already present.

----------------------------------------------------

Press Queries

If you or your institution will be preparing press materials for this manuscript, or if you need to know your paper's publication date for media purposes, please inform the journal staff as soon as possible so that your submission can be scheduled accordingly. Your manuscript will remain under a strict press embargo until the publication date and time. This means an early version of your manuscript will not be published ahead of your final version. PLOS Genetics may also choose to issue a press release for your article. If there's anything the journal should know or you'd like more information, please get in touch via plosgenetics@plos.org.

29 Mar 2021

PGENETICS-D-20-01014R3

Novel Variance-Component TWAS method for studying complex human diseases with applications to Alzheimer's dementia

Dear Dr Yang,

We are pleased to inform you that your manuscript entitled "Novel Variance-Component TWAS method for studying complex human diseases with applications to Alzheimer's dementia" has been formally accepted for publication in PLOS Genetics! Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out or your manuscript is a front-matter piece, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Genetics and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Alice Ellingham

PLOS Genetics

On behalf of:

The PLOS Genetics Team

Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom

plosgenetics@plos.org | +44 (0) 1223-442823

plosgenetics.org | Twitter: @PLOSGenetics