SMUT Package: Testing the Mediation Effect of Multiple SNPs on an Outcome Through a Mediator

2.1 Mediator model

M = α₁ + Xι^M + Gβ + ϵ₁ where M is a n by 1 vector of the mediator, α₁ is the intercept, G is a n by q matrix of the SNPs data, β is a q by 1 vector of SNPs’ effects on the mediator, X is a n by p matrix of covariates data, ι^M is a p by 1 vector of covariates effects’ on the mediator, ϵ₁ is a n by 1 vector of error terms in the mediator model.

2.2 Outcome model

Y = α₂ + Mθ + Xι^Y + Gγ + ϵ₂ where Y is a n by 1 vector of the outcome, α₂ is the intercept, M is the mediator, θ is the mediator’s effect on the outcome, X is the covariates data, ι^Y is a p by 1 vector of covariates effects’ on the outcome, G is the SNPs data, γ is a q by 1 vector of SNPs’ effects on the outcome, ϵ₂ is a n by 1 vector of error terms in the outcome model.

2.3 Testing β in the mediator model and θ in the outcome model

SMUT method adopts SKAT method (Wu et al., 2011) to test H₀^β : β = 0 versus H₁^β : β ≠ 0 in the mediator model and adopts Expectation-Maximization (EM) algorithm as well as score statistics to test H₀^θ : θ = 0 versus H₁^θ : θ ≠ 0 in the outcome model. Suppose the p value of testing β = 0 is p_β and the p value of testing θ = 0 is p_θ. SMUT method uses the hypotheses H₀ : βθ = 0 versus H₁ : βθ ≠ 0 to test the mediation effect. By leveraging IUT, the p value of the medition effect is p_IUT = max (p_β, p_θ).

2.4 Example of running R function SMUT

The example data has a genotype matrix (G) of 100 individuals with 200 SNPs.

library(SMUT)
dim(Genotype_data)
#> [1] 100 200
Genotype_data[1:3,1:4]
#>              SNP_1 SNP_2 SNP_3 SNP_4
#> Individual_1     0     2     1     0
#> Individual_2     1     1     0     0
#> Individual_3     1     1     0     1

We generate one continuous mediator and one continuous outcome based on this genotype matrix, as well as two covariates.

N_individual = nrow(Genotype_data)
N_SNPs = ncol(Genotype_data)

set.seed(1)

# generate two covariates 
X1 = rnorm(N_individual, 2, 3)
X2 = sample(c(0,1), N_individual, replace = TRUE) 
X = cbind(X1, X2)

# generate coefficients: iota_M, iota_Y, beta, theta and gamma
iota_M = c(0.3,0.5)
iota_Y = c(0.2,0.6)
beta = rnorm(N_SNPs, 1, 2)
theta = 1.2
gamma = rnorm(N_SNPs, 0.5, 2)

# generate error terms
e1 = rnorm(N_individual, 0, 1)
e2 = rnorm(N_individual, 0, 1)

# generate the mediator
mediator = 1 + X %*% iota_M + Genotype_data %*% beta + e1

# generate the outcome
outcome = 2 + mediator*theta + X %*% iota_Y + Genotype_data %*% gamma + e2

We apply SMUT function to test the mediation effect.

result_continuous = SMUT(G = Genotype_data, mediator = mediator, outcome = outcome, covariates = X)
#> Warning: 2 SNPs with either high missing rates or no-variation are excluded!
#> Warning: Genotypes of some variants are not the number of minor alleles! These
#> genotypes are flipped!
print( unlist( result_continuous ))
#>   p_value_IUT p_value_theta  p_value_beta 
#>  1.869722e-02  1.869722e-02  3.811493e-08

The warning messages are generated by the SKAT function in the R package SKAT. From the result, we can see that the p value of the mediation effect (p_value_IUT) is 0.02595523, which is the maximum of the p value of testing β (p_value_beta) and the p value of testing θ (p_value_theta).

2.5 References

Full details of the SMUT method can be found in the manuscript:

Zhong, W., Spracklen, C.N., Mohlke, K.L., Zheng, X., Fine, J. and Li, Y., 2019. Multi-SNP mediation intersection-union test. Bioinformatics, 35(22), pp.4724-4729.

3.1 Mediator model

The mediator model in the GSMUT method is the same as the one in the SMUT method. M = α₁ + Xι^M + Gβ + ϵ₁

3.2 Outcome model

3.3 Testing β in the mediator model and θ in the outcome model

GSMUT method adopts SKAT method to test H₀^β : β = 0 versus H₁^β : β ≠ 0 in the mediator model and adopts likelihood ratio statistics to test H₀^θ : θ = 0 versus H₁^θ : θ ≠ 0 in the outcome model. Suppose the p value of testing β = 0 is p_β and the p value of testing θ = 0 is p_θ. GSMUT method uses the hypotheses H₀ : βθ = 0 versus H₁ : βθ ≠ 0 to test the mediation effect. By leveraging IUT, the p value of the medition effect is p_IUT = max (p_β, p_θ).

In the manuscript, GSMUT is referred to as SMUT_GLM for a continuous, binary or count outcome and as SMUT_PH for a time-to-event outcome.

3.4 Example of running R function GSMUT

The example data has a genotype matrix (G) of 100 individuals with 200 SNPs.

library(SMUT)
dim(Genotype_data)
#> [1] 100 200
Genotype_data[1:3,1:4]
#>              SNP_1 SNP_2 SNP_3 SNP_4
#> Individual_1     0     2     1     0
#> Individual_2     1     1     0     0
#> Individual_3     1     1     0     1

We generate one continuous mediator, one binary outcome and one time-to-event outcome based on this genotype matrix, as well as two covariates.

N_individual = nrow(Genotype_data)
N_SNPs = ncol(Genotype_data)

set.seed(1)
# generate two covariates 
X1 = rnorm(N_individual, 2, 3)
X2 = sample(c(0,1), N_individual, replace = TRUE) 
X = cbind(X1, X2)

# generate coefficients: iota_M, iota_Y, beta, theta and gamma
iota_M = c(0.3,0.5)
iota_Y = c(0.2,-0.6)
beta = rnorm(N_SNPs, 0, 0.5)
theta = 1
gamma = rnorm(N_SNPs, 0, 0.3)

# generate error terms
e1 = rnorm(N_individual, 0, 1)

# generate the mediator
mediator = 1 + X %*% iota_M + Genotype_data %*% beta + e1

# generate the binary outcome 
eta = 2 + mediator*theta + X %*% iota_Y + Genotype_data %*% gamma
pi = 1/(1+exp( -(eta ) ))
binary_outcome = rbinom(length(pi),size=1,prob=pi) 

# generate the time-to-event outcome based on Weibull baseline hazard
v = runif(N_individual)
lambda=0.01; rho=1; rateC=0.01
Tlat = (- log(v) / (lambda * exp( eta  )))^(1 / rho)
# censoring times
C = rexp(N_individual, rate=rateC)
# follow-up times and event indicators
time = pmin(Tlat, C)
status = as.numeric(Tlat <= C)
survival_outcome = cbind(time,status)      
colnames(survival_outcome) = c("time","status")

We apply GSMUT function to test the mediation effect.

result_binary = GSMUT(G = Genotype_data, mediator = mediator, outcome = binary_outcome, covariates = X, outcome_type = "binary")
#> Warning: 2 SNPs with either high missing rates or no-variation are excluded!
#> Warning: Genotypes of some variants are not the number of minor alleles! These
#> genotypes are flipped!
print( unlist( result_binary ) )
#>   p_value_IUT p_value_theta     theta_hat  p_value_beta 
#>  5.780220e-04  5.361588e-05  9.119566e-01  5.780220e-04

result_survival = GSMUT(G = Genotype_data, mediator = mediator, outcome = survival_outcome, covariates = X, outcome_type = "survival")
#> Warning: 2 SNPs with either high missing rates or no-variation are excluded!
#> Warning: Genotypes of some variants are not the number of minor alleles! These
#> genotypes are flipped!
print( unlist( result_survival ))
#>   p_value_IUT p_value_theta     theta_hat  p_value_beta 
#>   0.025367552   0.025367552   0.159486545   0.000578022

The warning messages are generated by the SKAT function in the R package SKAT. From the result, we can see that the p value of the mediation effect (p_value_IUT) for the binary outcome is 0.000421673; p value of the mediation effect for the time-to-event outcome is 0.02470066. Here the p value of the mediation effect (p_value_IUT) is the maximum of the p value of testing β (p_value_beta) and the p value of testing θ (p_value_theta). And theta_hat is the point estimate of the θ in the outcome model.

3.5 Options for the parameter outcome_type in the GSMUT function

3.6 References

Full details of the GSMUT method can be found in the manuscript:

Zhong, W., Darville, T., Zheng, X., Fine, J. and Li, Y., 2019. Generalized Multi-SNP Mediation Intersection-Union Test. Submitted