### R code to replicate: ### ### Zhao Z, Banterle M, Bottolo L, Richardson S, Lewin A, Zucknick M (2020). ### BayesSUR: An R Package for High-Dimensional Multivariate Bayesian ### Variable and Covariance Selection in Linear Regression ### ### (Corresponding R package: https://CRAN.R-project.org/package=BayesSUR) ### ### zhi.zhao@medisin.uio.no, manuela.zucknick@medisin.uio.no ### Note: ### In this script, the plot() method for S3 class "BayesSUR" ### (BayesSUR:::plot.BayesSUR()) is called with the argument option "fig.tex = TRUE". ### This option requires that LaTeX is installed on the machine and will lead to ### errors if LaTeX is not installed. ### If this is the case, we therefore recommend to set the option to "fig.tex = FALSE". ################################################### ### 4. Quick start with a simple example ################################################### set.seed(82193) n <- 10; s <- 3; p <- 15 X <- matrix(rnorm(n * p, 2, 1), nrow = n) B <- matrix(c(0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, rep(0, s * p - 12)), nrow = p, byrow = TRUE) E <- matrix(rnorm(n * s, 0, 0.2), nrow = n) Y <- X %*% B + E B library("BayesSUR") library("tictoc") tic("Time of model fitting") set.seed(1294) fit <- BayesSUR(Y = Y, X = X, outFilePath = "results", output_CPO = TRUE) toc() fit plot(fit, estimator = c("beta", "gamma"), type = "heatmap", fig.tex = TRUE, output = "exampleEst", xlab = "Predictors", ylab = "Responses") G <- matrix(0, ncol = s * p, nrow = s * p) combn1 <- combn(rep((1:2 - 1) * p, each = length(2:3)) + rep(2:3, times = length(1:2)), 2) combn2 <- combn(rep((3-1) * p, each = length(c(1,4))) + rep(c(1,4), times = length(3)), 2) G[c(combn1[1,], combn2[1]), c(combn1[2,], combn2[2])] <- 1 tic("Time of model fitting") set.seed(5294) fit <- BayesSUR(Y = Y, X = X, outFilePath = "results", gammaPrior = "MRF", mrfG = G) toc() plot(fit, estimator = c("beta", "gamma"), type = "heatmap", fig.tex = TRUE, output = "exampleEst2", xlab = "Predictors", ylab = "Responses") ################################################### ### 5. Two extended examples based on real data ################################################### ### 5.1. Simulated eQTL data ################################################### data("exampleEQTL", package = "BayesSUR") str(exampleEQTL) attach(exampleEQTL) library("tikzDevice") tikz("ParamTrue.tex", width = 6, height = 3, standAlone = TRUE, packages = c("\\usepackage{tikz}", "\\usepackage{bm}", "\\usepackage[active,tightpage]{preview}", "\\PreviewEnvironment{pgfpicture}")) par(mfrow = c(1,2), mar = c(6, 6, 4, 2) + 0.1) image(z = gamma, x = 1:150, y = 1:10, col = grey(1:0), xlab = "SNPs Index", ylab = "Responses", main = paste("True", "$\\bm{\\gamma}$")); box() image(z = t(Gy), x = 1:10, y = 1:10, col = grey(1:0), xlab = "Responses", ylab = "Responses", main = paste("True", "$\\mathcal{G}$")); box() dev.off() tools::texi2pdf("ParamTrue.tex") ### Please note that running this example with the given MCMC parameter values ### (nIter, nChains, burnin) might take up to half an hour. set.seed(28173) tic("Time of model fitting") fit <- BayesSUR(data = data, Y = blockList[[1]], X = blockList[[2]], outFilePath = "results", nIter = 200000, nChains = 3, burnin = 100000, covariancePrior = "HIW", gammaPrior = "hotspot") toc() plot(fit, estimator = c("beta", "gamma", "Gy"), type = "heatmap", fig.tex = TRUE) pdf("ResponseGraph.pdf", height = 6, width = 12) layout(matrix(1:2, ncol = 2)) plot(fit, estimator = "Gy", type = "graph") plotGraph(Gy) dev.off() pdf("Network.pdf", height = 10, width = 10) plot(fit, estimator = c("gamma", "Gy"), type = "network", name.predictors = "SNPs", name.responses = "Gene expression") dev.off() pdf("Manhattan.pdf", height = 7, width = 8) plot(fit, estimator = "gamma", type = "Manhattan") dev.off() pdf("MCMCdiag.pdf", height = 8, width = 10) plot(fit, estimator = "logP", type = "diagnostics") dev.off() detach(exampleEQTL) ################################################### ### 5.2. The genomics of drug sensitivity in ### cancer data ################################################### data("exampleGDSC", package = "BayesSUR") attach(exampleGDSC) ### Please note that running this example with the given MCMC parameter values ### (nIter, nChains, burnin) might take several hours. hyperpar <- list(mrf_d = -3, mrf_e = 0.2) set.seed(6437) tic("Time of model fitting") fit <- BayesSUR(data = data, Y = blockList[[1]], X_0 = blockList[[2]], X = blockList[[3]], outFilePath = "results", nIter = 200000, burnin = 100000, nChains = 6, covariancePrior = "HIW", gammaPrior = "MRF", hyperpar = hyperpar, mrfG = mrfG) toc() plotEstimator(fit, estimator = "Gy", name.responses = c("Methotrexate", "RDEA119", "PD.0325901", "CI.1040", "AZD6244", "Nilotinib", "Axitinib"), fig.tex = TRUE, output = "ResponseGraphGDSC1") pdf("ResponseGraphGDSC2.pdf", height = 8, width = 8) plotGraph(fit, estimator = "Gy") dev.off() pdf("ResponseNetwork.pdf", height = 8, width = 8) plotNetwork(fit, estimator = c("gamma", "Gy"), label.predictor = "", name.predictors = "Genes", name.responses = "Drugs", nodesizePredictor = 2) dev.off() pdf("ResponseNetwork1.pdf", height = 8, width = 8) data("targetGene", package = "BayesSUR") plotNetwork(fit, estimator = c("gamma", "Gy"), includeResponse = c("RDEA119", "PD.0325901", "CI.1040", "AZD6244"), includePredictor = names(targetGene$group1)) dev.off() pdf("ResponseNetwork2.pdf", height = 4, width = 8) layout(matrix(1:2, ncol = 2)) plotNetwork(fit, estimator = c("gamma", "Gy"), edge.weight = TRUE, includeResponse = c("Nilotinib", "Axitinib"), includePredictor = names(targetGene$group2)) plotNetwork(fit, estimator = c("gamma", "Gy"), edge.weight = TRUE, PmaxPredictor = 0.01, includeResponse = c("Nilotinib", "Axitinib"), includePredictor = names(targetGene$group2)) dev.off()