MiSTi: set-based mixed effects score test for GxE interactions

Yu-Ru Su

2018-12-13

This is an illustration for using the package MiSTi to implement set-based GxE interaction test under mixed effects models. To form the input data, make sure that the columns of the data are arranged in the following order [Y E X G], where Y is an outcome variable, E is a matrix of environmental risk factors, X is a matrix of confounder(s), and G is a matrix of genotype(s). The current version of MiSTi allows multiple environmental risk factors (dimension is specified by option m), multiple confounders (dimension is specified by option d), and multiple variants (dimension is specified by option p).

Several R packages are called in MiSTi and should be installed beforehand. These include:

The main function to implement the integrative GxE analysis is MiSTi. The required inputs include

Other options to customize the analysis are available, including

The references for the methods are:

The following is an example of having simple burden score (weight being 1’s for all variants) and variance component test in MiSTi.

library(MiSTi)

# Y: binary outcome variable. a vector of length n.
# E: environmental risk factors. a vector of length n or a matrix of n rows.
# X: d confounders 
#    Either a vector of length n or a n*d matrix if d is greather than 1.
# G: genotypes of varaints. a n*p matrix.

# A special case with 1 environmental factor.
# Data generation
n = 2000
X = rbinom(n,size=1,prob=0.5)
E = rbinom(n,size=1,prob=0.5)
MAF = runif(10,min=0.001,max=0.01)
G = sapply(MAF,function(maf) rbinom(n,size=1,prob=maf))
eta = E*0.5 + G%*%rnorm(10,mean=0,sd=0.5) + (G*E)%*%rnorm(10,mean=0.5,sd=5)
Y = rbinom(n,size=1,prob=exp(eta)/(1+exp(eta)))
d = 1
p = 10

data = data.frame(Y=Y,E=E,X=X,G=G)
misti = MiSTi(data = data,
              outcome_type = "Binary",
              d = d,
              p = p,
              m = 1,
              R = 1,
              weight_method = "User"
)
misti$pvalue
misti$stat
misti$rho

# If only the burden and variance component tests are of interest:
misti.NoCombination = MiSTi(data = data,
              outcome_type = "Binary",
              d = d,
              p = p,
              m = 1,
              R = 1,
              weight_method = "User",
              combination_return = FALSE
)
misti.NoCombination$pvalue["BurdenComp"]
misti.NoCombination$pvalue["VarComp"]

# If only the Fisher's combination method is of interests:
misti.FisherCombination = MiSTi(data = data,
              outcome_type = "Binary",
              d = d,
              p = p,
              m = 1,
              R = 1,
              weight_method = "User",
              combination_return = TRUE,
              combination_preference = "Fisher"
)
misti.FisherCombination$pvalue["Fisher"]

If a functional annotation for the genetic variants are available, one can incorporate the functional annotation into the GxE test. Below is the sample code.

# FA: a vector of length p containing functional annotation for variants in consideration.

FA = runif(p)
misti.FAOnly = MiSTi(data = data,
                     outcome_type = "Binary",
                     d = d,
                     p = p,
                     m = 1,
                     R = 1,
                     weight_method = "User",
                     user_weight = FA
)

If the user would like to jointly test the effects of two weighted burden scores (e.g., a vector of 1’s and functional annotation), the sample code is shown below. In this case, the p-value of the joint BxE effects of the two burdens and their individual p-values after adjusting for the presence of each other will be returned by default.

misti.2burdens = MiSTi(data = data,
                       outcome_type = "Binary",
                       d = d,
                       p = p,
                       m = 1,
                       R = 2,
                       weight_method = "User",
                       user_weight = cbind(rep(1,p),FA)
)

Below we demonstrate sample code for situations under which there are multiple environmental factors to be tested, and their joint interaction effect with the genetic varaints are of interest. Extension to situations with functional annotations is straightforward by assigning values to the option ‘user_weight’.

E = sapply(c(0.5,0.5,0.5),function(p_e) rbinom(n,size=1,prob=p_e))
GE = do.call(cbind,sapply(1:ncol(E),simplify=FALSE,function(j) G*E[,j]))
eta = E%*%c(0,0,0) + G%*%rnorm(10,mean=0,sd=0.5) + GE%*%rnorm(ncol(GE),mean=0.5,sd=5)
Y = rbinom(n,size=1,prob=exp(eta)/(1+exp(eta)))
m = 3
data = data.frame(Y=Y,E=E,X=X,G=G)
misti.3Es = MiSTi(data = data,
                  outcome_type = "Binary",
                  d = d,
                  p = p,
                  m = m,
                  R = 1,
                  weight_method = "User"
)
misti.3Es$pvalue
misti.3Es$stat
misti.3Es$rho