Latent Growth Curve, Genome-Wide Association Study

Suppose you have longitudinal phenotypic measurements and genetic data on the same individuals. The gwsem package does not offer a latent growth curve model. Here we demonstrate how to build a latent growth curve model from scratch and perform a genome-wide association test on the slope. Hopefully this demonstration will inspire researchers to incorporate novel structural equation models into genome-wide association studies.

Design the Model

Visit the Onyx website. Download Onyx, and create the pictured model.

LGC Model

LGC Model

In this model, we have four occasions of measurement for the phenotype (t1, t2, t3, and t4). There are two latent variables, int for the intercept and slope for the slope of the trajectory. The loadings from int to the manifests are all fixed to 1. In contrast, the loadings from slope to the manifests increase. In OpenMx, it is possible to use the actual date of measurement instead of time increments of 1 unit. However, the use of non-unit time increments would force us to estimate the model using maximum likelihood, a method much slower than weighted least squares. So we recommend unit time increments. The single nucleotide polymorphism (SNP) data will appear in the snp indicator. It will be convienent to use exactly the name snp because gwsem knows which column to use by its name.

We export the model. This model does not entirely suit our purpose. After copying the code to a file, we remove the bits relating to modelData. We will also need to make a few changes to accomodate covariates. To make these changes easier, we gather up all the mxPath statements into a list.

model export

model export

library(gwsem)  # load gwsem
manifests<-c("t1","t2","t3","t4","snp")
latents<-c("int","slope")
path <- list(mxPath(from="int",to=c("t1","t2","t3","t4"), free=c(FALSE,FALSE,FALSE,FALSE), value=c(1.0,1.0,1.0,1.0) , arrows=1, label=c("int__t1","int__t2","int__t3","int__t4") ),
                 mxPath(from="slope",to=c("t1","t2","t3","t4"), free=c(FALSE,FALSE,FALSE,FALSE), value=c(0.0,1.0,2.0,3.0) , arrows=1, label=c("slope__t1","slope__t2","slope__t3","slope__t4") ),
                 mxPath(from="one",to=c("int","slope"), free=c(TRUE,TRUE), value=c(0.0,0.0) , arrows=1, label=c("const__int","const__slope") ),
                 mxPath(from="snp",to=c("slope","int"), free=c(TRUE,TRUE), value=c(1.0,0.0) , arrows=1, label=c("snp__slope","snp__int") ),
                 mxPath(from="int",to=c("int"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_int") ),
                 mxPath(from="slope",to=c("slope","int"), free=c(TRUE,TRUE), value=c(1.0,0.1) , arrows=2, label=c("VAR_slope","COV_slope_int") ),
                 mxPath(from="t1",to=c("t1"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_err") ),
                 mxPath(from="t2",to=c("t2"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_err") ),
                 mxPath(from="t3",to=c("t3"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_err") ),
                 mxPath(from="t4",to=c("t4"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_err") ),
                 mxPath(from="snp",to=c("snp"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("snp_res") ),
                 mxPath(from="one",to=c("t1","t2","t3","t4","snp"), free=F, value=0, arrows=1))
model <- mxModel("lgc", 
                 type="RAM",
                 manifestVars = manifests,
                 latentVars = latents,
                 path)

Simulate Fake Data

Even if we actually had real data, it might be useful to simulate fake data to develop expertise with our model and better understand what to expect from the real analyses. When we built the model in Onyx, we were careful to set starting values that were also appropriate for data simulation. Due to this preparation, data are easily simulated. We simulate N=500 people because this is the amount of fake SNP data available in the gwsem package (in the extdata subdirectory).

sim1 <- mxGenerateData(model, N)
head(sim1)
#>           t1         t2         t3         t4         snp
#> 1 -1.2127418  0.2550414  1.3196908  1.7728925  1.36860331
#> 2 -0.6847547  1.2415671 -0.1643415  2.3191331  0.02415591
#> 3 -0.4762370 -1.9330935 -2.8584654 -2.9570656 -0.78808290
#> 4 -1.4248013  0.4294114  0.6950322  0.6989669  1.38289021
#> 5 -1.3448042 -0.9842334  0.0289127  2.9854624  0.03050324
#> 6 -1.8998432  1.3816539  1.3972637  4.2704799  0.35110941

A placeholder snp column needs to exist in the sim1 data.frame, and will be replaced with each SNP when we run the GWAS. We add some covariates too.

for (ii in 1:5) {
  sim1[[paste0('pc', ii)]] <- rnorm(N)
}

Model Customization

To add covariates, we need to rebuild the model. To make it easy to rebuild the model, we saved all the mxPath statements in the path list (above). All these paths can be implemented in the new model by including the path list in the model construction. We add mxFitFunctionWLS to select the weighted least squares fit function and add our simulated data.

m2 <- mxModel("lgc", type="RAM",
        manifestVars = model$manifestVars,
        latentVars = c(model$latentVars, paste0('pc', 1:5)),
        path,
        mxExpectationRAM(M="M"),
        mxFitFunctionWLS(allContinuousMethod="marginals"),
        mxData(observed=sim1, type="raw", minVariance=0.1, warnNPDacov=FALSE))

At this point, m2 almost corresponds with the output of buildOneItem or similar model building function. We just need to set up the covariates.

m2 <- setupExogenousCovariates(m2, paste0('pc', 1:5), paste0('t',1:4))

The model is ready to run. If we had ordinal indicators then we could use setupThresholds to perform some additional set up.

Running the GWAS

We run the model on the first SNP as a sanity check to ensure that our model is working correctly. The GWAS function writes parameter estimates to the given out.log file, but it also returns the model back to R. We can call summary on this model to see the results from the first SNP.

tdir <- tempdir()
dir <- system.file("extdata", package = "gwsem")
snp1 <- GWAS(m2, file.path(dir,"example.pgen"), file.path(tdir, "out.log"), SNP=1)
#> Running lgc with 29 parameters
#> Done. See '/tmp/RtmpQ9nDAC/out.log' for results
summary(snp1)
#> Summary of lgc 
#>  
#> free parameters:
#>             name matrix   row   col      Estimate  Std.Error
#> 1       snp__int      A   int   snp -0.0349744341 0.06073971
#> 2     snp__slope      A slope   snp -0.0165135385 0.05826856
#> 3      pc1_to_t1      A    t1   pc1 -0.0634448796 0.06624284
#> 4      pc1_to_t2      A    t2   pc1 -0.0700728086 0.09463976
#> 5      pc1_to_t3      A    t3   pc1 -0.0824881406 0.14700949
#> 6      pc1_to_t4      A    t4   pc1  0.0801478609 0.20388143
#> 7      pc2_to_t1      A    t1   pc2  0.0262558711 0.05912795
#> 8      pc2_to_t2      A    t2   pc2  0.1656387661 0.08897793
#> 9      pc2_to_t3      A    t3   pc2  0.2549976994 0.14372096
#> 10     pc2_to_t4      A    t4   pc2  0.4455104212 0.19849419
#> 11     pc3_to_t1      A    t1   pc3  0.0669053447 0.06744135
#> 12     pc3_to_t2      A    t2   pc3  0.0726561974 0.09170514
#> 13     pc3_to_t3      A    t3   pc3  0.1605438466 0.13886029
#> 14     pc3_to_t4      A    t4   pc3  0.3595074167 0.19965554
#> 15     pc4_to_t1      A    t1   pc4 -0.1193209769 0.06529782
#> 16     pc4_to_t2      A    t2   pc4 -0.0702491986 0.09741140
#> 17     pc4_to_t3      A    t3   pc4  0.0309637039 0.15472166
#> 18     pc4_to_t4      A    t4   pc4  0.1371833370 0.22218885
#> 19     pc5_to_t1      A    t1   pc5 -0.0261296915 0.06551636
#> 20     pc5_to_t2      A    t2   pc5 -0.0629219036 0.09251596
#> 21     pc5_to_t3      A    t3   pc5  0.0498048070 0.15291992
#> 22     pc5_to_t4      A    t4   pc5  0.0006757091 0.21240210
#> 23       VAR_err      S    t1    t1  1.1558779066 0.06294267
#> 24       snp_res      S   snp   snp  0.4149100566 0.01344624
#> 25       VAR_int      S   int   int  0.8758430256 0.12115942
#> 26 COV_slope_int      S   int slope  0.0407914368 0.08860477
#> 27     VAR_slope      S slope slope  1.8843513886 0.12169536
#> 28    const__int      M     1   int  0.0865348717 0.06049005
#> 29  const__slope      M     1 slope  0.0140700095 0.06821112
#> 
#> Model Statistics: 
#>                |  Parameters  |  Degrees of Freedom  |  Fit (r'Wr units)
#>        Model:             29                   2471              181.499
#>    Saturated:             20                   2480                0.000
#> Independence:             10                   2490                   NA
#> Number of observations/statistics: 500/2500
#> 
#> chi-square:  χ² ( df=16 ) = 181.499,  p = 4.216078e-30
#> CFI: NA 
#> TLI: NA   (also known as NNFI) 
#> RMSEA:  0.1438311  [95% CI (0.1217977, 0.1666155)]
#> Prob(RMSEA <= 0.05): 3.330669e-16
#> To get additional fit indices, see help(mxRefModels)
#> timestamp: 2020-03-27 15:27:53 
#> Wall clock time: 0.5825934 secs 
#> OpenMx version number: 2.17.3.14 
#> Need help?  See help(mxSummary)

This looks alright so we proceed to run all the SNPs and retrieve the results. The loadResults function calculates Z scores and P values for the snp__slope parameter. We might also examine the snp__int parameter.

GWAS(m2, file.path(dir,"example.pgen"), file.path(tdir, "out.log"))
#> Running lgc with 29 parameters
#> 
[00:21] MxComputeGradientDescent(SLSQP) evaluations 18616 fit 1820.21 change -0.07378
[00:20] MxComputeGradientDescent(SLSQP) evaluations 32813 fit 1097.86 change -407.8  
[00:17] MxComputeGradientDescent(SLSQP) evaluations 50839 fit 1346.89 change 1255  
[00:16] MxComputeGradientDescent(SLSQP) evaluations 67855 fit 1895.49 change -0.2962
[00:16] MxComputeGradientDescent(SLSQP) evaluations 85182 fit 7984.11 change 7936   
[00:15] MxComputeGradientDescent(SLSQP) evaluations 103684 fit 8729.33 change 7464
[00:14] MxComputeGradientDescent(SLSQP) evaluations 122367 fit 1444.28 change -3.018e-05
[00:13] MxComputeGradientDescent(SLSQP) evaluations 138445 fit 1425.19 change 774       
[00:11] MxComputeGradientDescent(SLSQP) evaluations 157006 fit 6935.51 change 4334
[00:10] MxComputeGradientDescent(SLSQP) evaluations 175627 fit 8353.27 change 8140
[00:09] MxComputeGradientDescent(SLSQP) evaluations 196619 fit 2605.49 change 1328
[00:07] MxComputeGradientDescent(SLSQP) evaluations 220477 fit 107.103 change -1.13
[00:06] MxComputeGradientDescent(SLSQP) evaluations 241628 fit 526.335 change -4.954
[00:05] MxComputeGradientDescent(SLSQP) evaluations 261367 fit 14162.2 change 1.387e+04
[00:04] MxComputeGradientDescent(SLSQP) evaluations 279455 fit 8192.44 change 6165     
[00:03] MxComputeGradientDescent(SLSQP) evaluations 300557 fit 199.181 change -0.002179
[00:02] MxComputeGradientDescent(SLSQP) evaluations 321369 fit 3060.18 change 1309     
[00:01] MxComputeGradientDescent(SLSQP) evaluations 342780 fit 1277.28 change 995.3
[00:00] MxComputeGradientDescent(SLSQP) evaluations 361459 fit 1574.97 change -535.2
                                                                                    
#> Done. See '/tmp/RtmpQ9nDAC/out.log' for results
got <- loadResults(file.path(tdir, "out.log"), 'snp__slope')

Some of the models failed to fit. If we’re interested to discover why then we can load these separately,

susp <- loadSuspicious(file.path(tdir, "out.log"), 'snp__slope')
head(susp)
#>    MxComputeLoop1 CHR    BP     SNP A1 A2 statusCode
#> 1:              7   1  8000  RSID_8  G  A       <NA>
#> 2:              9   1 10000 RSID_10  G  A       <NA>
#> 3:             16   1 17000 RSID_17  G  A       <NA>
#> 4:             18   1 19000 RSID_19  G  A       <NA>
#> 5:             20   1 21000 RSID_21  G  A       <NA>
#> 6:             23   1 24000 RSID_24  G  A       <NA>
#>                                                 catch1 snp__slope snp__slopeSE
#> 1: lgc.data: 'snp' has observed variance less than 0.1          1          NaN
#> 2: lgc.data: 'snp' has observed variance less than 0.1          1          NaN
#> 3: lgc.data: 'snp' has observed variance less than 0.1          1          NaN
#> 4: lgc.data: 'snp' has observed variance less than 0.1          1          NaN
#> 5: lgc.data: 'snp' has observed variance less than 0.1          1          NaN
#> 6: lgc.data: 'snp' has observed variance less than 0.1          1          NaN
#>      Z   P
#> 1: NaN NaN
#> 2: NaN NaN
#> 3: NaN NaN
#> 4: NaN NaN
#> 5: NaN NaN
#> 6: NaN NaN

It looks like these SNPs had a variance below the threshold we specified in mxData (minVariance=0.1).

We can present a Manhattan plot consisting of the models that succeeded,

plot(got)