gwsem
Mehak Chopra
2024-06-07
library(OpenMx)## To take full advantage of multiple cores, use:
## mxOption(key='Number of Threads', value=parallel::detectCores()) #now
## Sys.setenv(OMP_NUM_THREADS=parallel::detectCores()) #before library(OpenMx)install.packages("gwsem")## Installing package into '/home/mchopra/R/x86_64-pc-linux-gnu-library/4.2'
## (as 'lib' is unspecified)## Warning: package 'gwsem' is not available for this version of R
##
## A version of this package for your version of R might be available elsewhere,
## see the ideas at
## https://cran.r-project.org/doc/manuals/r-patched/R-admin.html#Installing-packageslibrary(gwsem)##
## Attaching package: 'gwsem'## The following object is masked from 'package:base':
##
## signiflocation <- 'https://jpritikin.github.io/gwsem/gwsemGxEexample'gxeData <- read.table(file.path(location, "gxeData.txt"), header=TRUE)
gxeData[1:5, 1:5]## phe mod pc1 pc2 pc3
## 1 1.03273603 -0.3487804 0.7691604 -0.28744876 -0.04961699
## 2 0.79232013 -1.0257007 -1.0626999 0.03717712 0.46603262
## 3 0.08797211 -0.9074263 -1.0965822 -1.34888323 0.94100891
## 4 0.66951518 0.7746067 -0.1920381 -0.76406130 0.95443043
## 5 0.45793053 -0.9498979 -0.7118469 0.44507673 -1.65936741#includes the dependent variable, the moderating/interacting variable and 5 PCs)
#phe is the dependent varibale we are tryong to predict or understand.
#mod is the moderating/interacting variable. It could be age or sex depending on the situation.
#In detail, if you were investigating the effect of a treatment (IV) on a health outcome (DV), and you believed that the effect of the treatment might vary depending on the age of the participants, then age would be the moderating/interacting variable. Similarly, if you were examining the effect of a counseling intervention (IV) on reducing anxiety (DV), and you believed that the effect might be different for males and females, then gender would be the moderating/interacting variable.If want to chop the indicators up into binary or ordinal variables then this is the time to do it. Once you are satisfied the data are in the appropriate shape we can build a one factor GWAS model with the following command:
gxe <- buildItem(phenoData = gxeData, depVar = "phe",
covariates = c("mod", "pc1", "pc2", "pc3", "pc4", "pc5"),
fitfun = "WLS", exogenous = T, gxe = "mod")
gxe## MxModel 'OneItem'
## type : RAM
## $matrices : 'A', 'S', 'F', and 'M'
## $algebras : 'snp_modAlg'
## $penalties : NULL
## $constraints : NULL
## $intervals : NULL
## $latentVars : 'snp', 'snp_mod', 'mod', 'pc1', 'pc2', 'pc3', 'pc4', and 'pc5'
## $manifestVars : 'phe'
## $data : 6000 x 9
## $data means : NA
## $data type: 'raw'
## $submodels : NULL
## $expectation : MxExpectationRAM
## $fitfunction : MxFitFunctionWLS
## $compute : NULL
## $independent : FALSE
## $options :
## $output : FALSEgxeFit <- mxRun(gxe)## Running OneItem with 10 parameterssummary(gxeFit)## Summary of OneItem
##
## free parameters:
## name matrix row col Estimate Std.Error lbound ubound
## 1 snp_to_phe A phe snp -0.006375336 0.01849254
## 2 snp_mod_to_phe A phe snp_mod 0.021851971 0.01800641
## 3 mod_to_phe A phe mod 0.578280571 0.02171921
## 4 pc1_to_phe A phe pc1 0.003691416 0.01289730
## 5 pc2_to_phe A phe pc2 0.011251630 0.01324260
## 6 pc3_to_phe A phe pc3 -0.006366428 0.01280235
## 7 pc4_to_phe A phe pc4 -0.001726108 0.01329703
## 8 pc5_to_phe A phe pc5 0.004182387 0.01349019
## 9 phe_res S phe phe 1.035285847 0.01899422 0.001
## 10 pheMean M 1 phe 0.517041155 0.02247526
##
## Model Statistics:
## | Parameters | Degrees of Freedom | Fit (r'Wr units)
## Model: 10 0 NA
## Saturated: 2 8 0
## Independence: 2 8 NA
## Number of observations/statistics: 6000/10
##
## chi-square: χ² ( df=0 ) = 0, p = 1
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2024-07-17 15:44:37
## Wall clock time: 0.03088212 secs
## optimizer: SLSQP
## OpenMx version number: 2.21.11
## Need help? See help(mxSummary)library(curl)## Using libcurl 7.81.0 with GnuTLS/3.7.3#pgen file
curl_download(file.path(location, 'example.pgen'),
file.path(tempdir(),'example.pgen'))
#pvar file
curl_download(file.path(location, 'example.pvar'),
file.path(tempdir(),'example.pvar'))
getwd()## [1] "/home/mchopra/Documents/PhD-Year1/gwsem/gwsem"GWAS(model = gxe, # the model object
snpData = file.path(tempdir(), 'example.pgen'),
out=file.path(tempdir(), "gxe.log")) ## Running OneItem with 10 parameters## Done. See '/tmp/RtmpqueyVK/gxe.log' for resultsgxeResult <- read.delim(file.path(tempdir(), "gxe.log"))
gxeResult[1:5, 1:5]## OpenMxEvals iterations timestamp MxComputeLoop1 snp_to_phe
## 1 421 20 Jul 17 2024 03:44:40 PM 1 -0.03312635
## 2 819 39 Jul 17 2024 03:44:40 PM 2 0.03210928
## 3 1239 59 Jul 17 2024 03:44:40 PM 3 -0.02560457
## 4 1639 78 Jul 17 2024 03:44:40 PM 4 0.02024164
## 5 1995 95 Jul 17 2024 03:44:40 PM 5 0.01552086lr <- loadResults(path = file.path(tempdir(), "gxe.log"), focus = c("snp_to_phe", "snp_mod_to_phe"))
succinctCond <- signif(lr, "snp_to_phe")
succinctInt <- signif(lr, "snp_mod_to_phe")head(succinctCond[order(succinctCond$Z, decreasing = T),])## MxComputeLoop1 CHR BP SNP A1 A2 statusCode catch1
## <int> <int> <int> <char> <char> <char> <char> <lgcl>
## 1: 1884 1 1883 snp1883 A B OK NA
## 2: 1762 1 1761 snp1761 B A OK NA
## 3: 1868 1 1867 snp1867 B A OK NA
## 4: 48 1 47 snp47 B A OK NA
## 5: 1101 1 1100 snp1100 B A OK NA
## 6: 319 1 318 snp318 B A OK NA
## snp_to_phe snp_mod_to_phe Vsnp_to_phe:snp_to_phe
## <num> <num> <num>
## 1: 0.18781103 0.131664012 0.0004968010
## 2: 0.14568877 0.090479043 0.0004853252
## 3: 0.12841686 0.004175709 0.0004595031
## 4: 0.11862335 0.029088623 0.0004997281
## 5: 0.08049730 -0.004278833 0.0005016595
## 6: 0.07681582 -0.025794074 0.0005022668
## Vsnp_mod_to_phe:snp_mod_to_phe Vsnp_mod_to_phe:snp_to_phe Z
## <num> <num> <num>
## 1: 0.0004906852 5.606327e-06 8.426164
## 2: 0.0004858684 -8.634754e-06 6.613170
## 3: 0.0004469476 -1.380366e-05 5.990703
## 4: 0.0004998344 -1.980735e-05 5.306441
## 5: 0.0004887322 -8.432265e-06 3.593989
## 6: 0.0005077531 -9.825706e-06 3.427547
## P
## <num>
## 1: 3.571827e-17
## 2: 3.761773e-11
## 3: 2.089354e-09
## 4: 1.117866e-07
## 5: 3.256532e-04
## 6: 6.090600e-04plot(succinctInt) #snp_mod_to_phenotype # To plot the interaction coefficient
plot(succinctCond) ##snp_to_phenotype ## To plot p-values of the conditional effect