Fast and powerful genome wide association of dense genetic data with high dimensional imaging phenotypes

Genome wide association (GWA) analysis of brain imaging phenotypes can advance our understanding of the genetic basis of normal and disorder-related variation in the brain. GWA approaches typically use linear mixed effect models to account for non-independence amongst subjects due to factors, such as family relatedness and population structure. The use of these models with high-dimensional imaging phenotypes presents enormous challenges in terms of computational intensity and the need to account multiple testing in both the imaging and genetic domain. Here we present a method that makes mixed models practical with high-dimensional traits by a combination of a transformation applied to the data and model, and the use of a non-iterative variance component estimator. With such speed enhancements permutation tests are feasible, which allows inference on powerful spatial tests like the cluster size statistic.

In fact, the Python version of FaST-LMM implements the simplified REML for fixed effect testing with an LRT, integrating out nuisance fixed effects as in REML and performs ML estimation on the SNP effect that is being tested only. I hope you find these comments useful.
Reviewer #2 (Remarks to the Author): The authors present a method for linear mixed modelling for neuroimaging phenotypes that scales well compared to existing approaches.
The authors claim two major contributions. First, they claim their use of score tests is novel, although the approach they use is very similar to their previous work (reference 40). The use of Wald, Likelihood Ratio and Score tests for a linear mixed model are implemented in the GEMMA software by Zhou and Stephens (see the software manual for this work).
The second claim of novelty that they reduce complexity by eigen-rotating the model is not true. The simplified ML function is not due to reference 40, or to reference 33 that is quoted in that paper. The use of the eigenvectors of the GRM to diagonalize the covariance structure of the model is due to references 24 and 26. This has been a standard method now for several years. The authors have not correctly credited this earlier work and have overstated the use of this approach as novel.
The authors use permutation to account for the correlation between the tested phenotypes. There are methods that can be used to estimate an effective number of tests (see for example https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3325408/) which could be applied here to estimate the correction for the 42 ROIs analysed. How does applying this method compare to using permutation?
A standard approach in LMM fitting has been to fit a model without any SNPs as fixed effects to estimate the variance components. These are then kept fixed when testing each SNP. This was first suggested in the EMMAX approach. Any covariates can be regressed out in advance so that a residual phenotype is used. I would like to see how the authors approach compares to this approach on the simulated and real datasets.
How does the permutation scheme used compare to that used in reference 46. If it is the same that should be acknowledged. It is not, then a comparison should be included in the paper.
We thank the editor and the reviewers for their constructive comments. We have made significant revisions in response, marking our changes to the manuscript in blue.
Below we provide point-by-point responses, typesetting the reviewer/editor's comments in italics.
Editor's comment: We hope you will find the referees' comments useful as you decide how to proceed. We would like to point out that for further consideration in Nature Communications, we would have to see more extensive benchmarking as requested by Reviewer #2 and demonstration of superiority (not just in speed but also accuracy) over these other methods.
In addition to our method being able to scale to the large scale phenotypes found with imaging, we have shown that we have more accurate standard errors (less biased) than FaST-LMM (FIGURE 8 Supplementary material) and can obtain exact, permutation-based familywise error (FWE) for spatial inferences (voxel-wise and cluster-wise), something no other method is able to provide.

Reviewer #1
1. In order for the speedups to work, all traits need to be observed for all individuals. As this assumption does not hold in many other multiple-phenotype scenarios, the assumption should be stated in the paper explicitly.
Response: We have addressed this point and mentioned the assumption explicitly in discussion (p8 last paragraph) and section 1.4 of online method (p16).   Response: We apologize that we did not state original contributions more clearly. We have made changes to make our contribution clear (p2). Our proposed method reduces the complexity of LMM in GWA of high dimensional imaging data as follows: a. Use of non-iterative estimator for the random effects variance instead of fully converged estimator based on REML; and b. Vectorisation of score test computation, feasible due to the simplified REML function, vastly reducing the computational complexity of association testing. The two in combination accelerate the model fitting process sufficiently to allow permutation inference, important for FWE control of spatial inference tools like the cluster size test.

In
To be clear, we did not intend to imply covariance matrix diagonalization or use of a score statistic as our major contribution. While, the non-iterative random effect estimator was introduced in our previous work, that was based on the simplified ML function, and here we have used the simplified REML function; in particular, we have shown that this improves the bias and mean square error of the estimated variance. Response: We thank the reviewer for this reference. We have calculated a FWE threshold using the method proposed by the reviewer (p6). Application of a Bonferroni correction with effective number of tests yields a slightly more stringent threshold, but finds the same 8 association statistics as found with permutation. While this result would seem to undermine the need for permutation, we have noted that (1) this use of the effective number of independent tests depends on an arbitrary GWAS threshold that depends on the chip used, (2) the approach is not applicable to cluster-wise inference and (3) doesn't scale to voxel-wise inference which are the standard spatial statistics in imaging.

3.
A standard approach in LMM fitting has been to fit a model without any SNPs as fixed effects to estimate the variance components. These are then kept fixed when testing each SNP. This was first suggested in the EMMAX approach. Any covariates can be regressed out in advance so that a residual phenotype is used. I would like to see how the authors approach compares to this approach on the simulated and real datasets.
Response: The reviewer raises a valid point, as a model based on precomputed residualised phenotypes would indeed save computation. However, we have found that this approach is not safe in general, and have added a simulation study to address this point (p6). First, we note that imaging association studies routinely use intracranial volume (ICV) as a nuisance covariate, and ICV is well known to have high heritability (h 2 ranges: 0.81-0.91, Peper et al 2007, 0.72-0.92, Glahn et al. 2007), hence we conducted a simulation study that evaluates controlling for a heritable fixed effect nuisance covariate in the null simulation setting when there is neither a SNP effect nor a covariate effect on the phenotype. Figure 4 compares performance of FaST-LMM, EMMAX and the score test based on the simplified REML function using non-iterative random effect estimator (NINGA). The simulation results show that the parametric Pvalue from all approaches when the nuisance covariate is included in the LMM is valid (Figure 4 left panel). However, the null distribution of P-values from LMM fitted on residualised phenotypes can be conservative (Figure 4 right panel).
4. How does the permutation scheme used compare to that used in reference 46. If it is the same that should be acknowledged. It is not, then a comparison should be included in the paper.