Compute Bayesian credible interval for mediating effect

Data and packages

library(medinome)
library(limma)
library(qvalue)
library(MCMCpack)

Heart vs. Kidney data

Take from Lauren’s data comparing kidney versus heart tissue samples.

df <- get(load("../data/example-kidney-v-heart-human.rda"))
Y=df$exprs_pair
X=df$tissue
M=df$methyl_pair

# scale gene expression to mean 0 and standard deviation 1 for each gene
Y <- t(scale(t(Y)))
M <- t(scale(t(M)))

# apply limma to scaled expression to select DE genes
design_1 <- model.matrix(~X)
Y_counts <- 2^Y
Y_voom <- voom(Y_counts, design=design_1, normalize.method = "none")
model_1 <- lmFit(Y_voom, design_1)
model_1 <- eBayes(model_1)

# qvaule package for FDR control
qval <- qvalue(model_1$p.value[,2])

# select DE genes with q-value  < .05
ii.de <- which(qval$qvalues < .05)
Y.de <- Y[ii.de, ]
M.de <- M[ii.de, ]

# create a permutated dataset
M.permute <- do.call(rbind, lapply(1:nrow(M.de), function(g) {
  m.perm <- sample(M.de[g,])
  return(m.perm)
  }))

# fitting mediation test
fit.voom <- mediate.test.voom(Y=Y.de, X=X, M=M.de)

fit.voom.perm <- mediate.test.voom(Y=Y.de, X=X, M=M.permute)

Function to compute Bayesian credible interval.

ci <- vector("list", nrow(Y.de))
type_1 <- 0.05
prior_mn_model_3 <- rep(0,3)
prior_mn_model_2 <- rep(0,2)

for (g in 1:nrow(Y.de)) {

  prior_sd_model_2 <- c(fit.voom$sigma_gamma2[g], fit.voom$sigma_alpha[g])
  prior_sd_model_3 <- c(fit.voom$sigma_gamma3[g], fit.voom$sigma_tau_prime[g],
                        fit.voom$sigma_beta[g])
  
  M_g <- M.de[g,]
  Y_g <- Y.de[g,]

  model_2 <- MCMCregress(M_g ~ X, burnin = 1000 , mcmc = 1000, 
              seed = 7, b0= prior_mn_model_2, B0 = prior_sd_model_2^(-1))
              # sigma.mu = fit.voom$sigma_model2[g]^2,
              # sigma.var = var(fit.voom$sigma_model2^2))
  
  model_3 <- MCMCregress(Y_g ~ X + M_g, burnin = 1000 , mcmc = 1000, 
              seed = 7, b0= prior_mn_model_3, B0 = prior_sd_model_3^(-1))
              # sigma.mu = fit.voom$sigma_model3[g]^2,
              # sigma.var = var(fit.voom$sigma_model3^2))

  ab <- model_2[,2]*model_3[,3]

  ci.lower <- quantile(ab,probs=type_1/2,type=4)
  ci.higher <- quantile(ab,probs=(1-type_1/2),type=4)

  ci[[g]] <- data.frame(ci.lower, ci.higher)
}
ci <- do.call(rbind, ci)

res <- sapply(1:g, function(g) {fit.voom$ab[g] < ci[g,1] | fit.voom$ab[g] > ci[g,2]})
mean(res)

[1] 0.03683763

sum(res)

[1] 130

ci.perm <- vector("list", nrow(Y.de))
type_1 <- 0.05
prior_mn_model_3 <- rep(0,3)
prior_mn_model_2 <- rep(0,2)

for (g in 1:nrow(Y.de)) {

  prior_sd_model_2 <- c(fit.voom.perm$sigma_gamma2[g], fit.voom.perm$sigma_alpha[g])
  prior_sd_model_3 <- c(fit.voom.perm$sigma_gamma3[g], 
                        fit.voom.perm$sigma_tau_prime[g],
                        fit.voom.perm$sigma_beta[g])
  
  M_g <- M.permute[g,]
  Y_g <- Y.de[g,]

  model_2 <- MCMCregress(M_g ~ X, burnin = 1000 , mcmc = 1000, 
              seed = 7, b0= prior_mn_model_2, B0 = prior_sd_model_2^(-1))
              # sigma.mu = fit.voom.perm$sigma_model2^2,
              # sigma.var = var(fit.voom.perm$sigma_model2^2))
  
  model_3 <- MCMCregress(Y_g ~ X + M_g, burnin = 1000 , mcmc = 1000, 
              seed = 7, b0= prior_mn_model_3, B0 = prior_sd_model_3^(-1))
              # sigma.mu = fit.voom.perm$sigma_model3^2,
              # sigma.var = var(fit.voom.perm$sigma_model3^2))

  ab <- model_2[,2]*model_3[,3]

  ci.lower <- quantile(ab,probs=type_1/2,type=4)
  ci.higher <- quantile(ab,probs=(1-type_1/2),type=4)

  ci.perm[[g]] <- data.frame(ci.lower, ci.higher)
}
ci.perm <- do.call(rbind, ci.perm)

res.perm <- sapply(1:g, function(g) {fit.voom.perm$ab[g] < ci.perm[g,1] | fit.voom.perm$ab[g] > ci.perm[g,2]})
mean(res.perm)

[1] 0.01558515

Conclusions

• Preliminary results show that only 3 percent of the genes, there’s significant effect of mediation. Using permuted methylation data, we found significant effect in only 1 percent of the genes. Though this trend is consistent with the expectation, I have little confidence in this approach given the low percentage of genes found signficant in the test of mediation coupled with the potentiality for confounding between the relations between methylation and expression.

Session information

sessionInfo()

R version 3.4.1 (2017-06-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ltm_1.1-0        polycor_0.7-9    msm_1.6.6        assertthat_0.2.0
 [5] MCMCpack_1.4-2   MASS_7.3-49      coda_0.19-1      qvalue_2.10.0   
 [9] limma_3.34.9     medinome_0.0.1   ggplot2_2.2.1   

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.15       pillar_1.2.1       compiler_3.4.1    
 [4] git2r_0.21.0       plyr_1.8.4         tools_3.4.1       
 [7] digest_0.6.15      evaluate_0.10.1    tibble_1.4.2      
[10] gtable_0.2.0       lattice_0.20-35    rlang_0.2.0       
[13] Matrix_1.2-12      yaml_2.1.18        expm_0.999-2      
[16] mvtnorm_1.0-7      SparseM_1.77       stringr_1.3.0     
[19] knitr_1.20         MatrixModels_0.4-1 rprojroot_1.3-2   
[22] grid_3.4.1         survival_2.41-3    rmarkdown_1.9     
[25] reshape2_1.4.3     magrittr_1.5       backports_1.1.2   
[28] scales_0.5.0       htmltools_0.3.6    mcmc_0.9-5        
[31] splines_3.4.1      colorspace_1.3-2   quantreg_5.35     
[34] stringi_1.1.6      lazyeval_0.2.1     munsell_0.4.3