# Simulated data set to mimic a small expression quantitative trait loci (eQTL) example Simulated data set to mimic a small expression quantitative trait loci (eQTL) example, with p=150 single nucleotide polymorphisms (SNPs) as explanatory variables, s=10 gene expression features as response variables and data for n=100 observations. Loading the data will load the associated blockList object needed to fit the model with BayesSUR(). The R code for generating the simulated data is given in the Examples paragraph. #importFrom BDgraph rgwish #importFrom gRbase mcsMAT #importFrom scrime simulateSNPs data ## Format An object of class `list` of length 4. ```r exampleEQTL ``` ## Examples ```r # Load the eQTL sample dataset data("exampleEQTL", package = "BayesSUR") str(exampleEQTL) ## Not run: # =============== # The code below is to show how to generate the dataset "exampleEQTL.rda" above # =============== requireNamespace("BDgraph", quietly = TRUE) requireNamespace("gRbase", quietly = TRUE) requireNamespace("scrime", quietly = TRUE) ########################### Problem Dimensions n <- 100 p <- 150 s <- 10 ############################ Select a set of n x p (SNPs) covariates ## The synthetic data in the paper use a subset of the real SNPs as covariates, # but as the NFBC66 dataset is confidential we'll use scrime to sample similar data x <- scrime::simulateSNPs(c(n, 10), p, c(3, 2), prop.explain = c(0.9, 0.95))$data[1:n, ] x <- cbind(rep(1, n), x) #################################################################### graph_pattern <- 2 snr <- 25 corr_param <- 0.9 ### Create the underlying graph if (graph_pattern == 1) { ### 1) Random but full G <- matrix(1, s, s) Prime <- list(c(1:s)) Res <- Prime Sep <- list() } else if (graph_pattern == 2) { ### 2) Block Diagonal structure Prime <- list( c(1:floor(s * 2 / 3)), c((floor(s * 2 / 3) + 1):(ceiling(s * 4 / 5) - 1)), c(ceiling(s * 4 / 5):s) ) Res <- Prime Sep <- lapply(Res, function(x) which(x == -99)) G <- matrix(0, s, s) for (i in Prime) { G[i, i] <- 1 } } else if (graph_pattern == 3) { ### 3) Decomposable model Prime <- list( c(1:floor(s * 5 / 12), ceiling(s * 9 / 10):s), c(floor(s * 2 / 9):(ceiling(s * 2 / 3) - 1)), c(ceiling(s * 2 / 3):(ceiling(s * 4 / 5) - 1)), c(ceiling(s * 4 / 5):s) ) Sep <- list() H <- list() for (i in 2:length(Prime)) { H <- union(H, Prime[[i - 1]]) Sep[[i - 1]] <- intersect(H, Prime[[i]]) } Res <- list() Res[[1]] <- Prime[[1]] for (i in 2:length(Prime)) { Res[[i]] <- setdiff(Prime[[i]], Sep[[i - 1]]) } G <- matrix(0, s, s) for (i in Prime) { G[i, i] <- 1 } ## decomp check dimnames(G) <- list(1:s, 1:s) length(gRbase::mcsMAT(G - diag(s))) > 0 } else if (graph_pattern == 4) { ### 4) Non-decomposable model nblocks <- 5 nElemPerBlock <- c( floor(s / 4), floor(s / 2) - 1 - floor(s / 4), ceiling(s * 2 / 3) - 1 - floor(s / 2), 7 ) nElemPerBlock <- c(nElemPerBlock, s - sum(nElemPerBlock)) res <- 1:s blockIdx <- list() for (i in 1:nblocks) { # blockIdx[[i]] = sample(res,nElemPerBlock[i]) blockIdx[[i]] <- res[1:nElemPerBlock[i]] res <- setdiff(res, blockIdx[[i]]) } G <- matrix(0, s, s) ## add diagonal for (i in 1:nblocks) { G[blockIdx[[i]], blockIdx[[i]]] <- 1 } ## add cycle G[blockIdx[[1]], blockIdx[[2]]] <- 1 G[blockIdx[[2]], blockIdx[[1]]] <- 1 G[blockIdx[[1]], blockIdx[[5]]] <- 1 G[blockIdx[[5]], blockIdx[[1]]] <- 1 G[blockIdx[[2]], blockIdx[[3]]] <- 1 G[blockIdx[[3]], blockIdx[[2]]] <- 1 G[blockIdx[[3]], blockIdx[[5]]] <- 1 G[blockIdx[[5]], blockIdx[[3]]] <- 1 ## decomp check dimnames(G) <- list(1:s, 1:s) length(gRbase::mcsMAT(G - diag(s))) > 0 # Prime = blockIdx Res <- blockIdx ## this is not correct but not used in the non-decomp case } ### Gamma Pattern gamma <- matrix(0, p + 1, s) gamma[1, ] <- 1 ### 2) Extra Patterns ## outcomes (correlated in the decomp model) have some predictors in common gamma[6:10, 6:9] <- 1 ## outcomes (correlated in the decomp model) have some predictors in common # gamma[16:20,14:15] = 1 ## outcomes (sort-of correlated [pair-wise] in the decomp model) # have predictors in common 6:15 gamma[26:30, 4:8] <- 1 ## outcomes (NOT correlated in the decomp model) have predictors in common 16:17 gamma[36:40, c(3:5, 9:10)] <- 1 ## these predictors are associated with ALL the outcomes gamma[46:50, ] <- 1 combn11 <- combn(rep((6:9 - 1) * p, each = length(6:10 - 1)) + rep(6:10 - 1, times = length(6:9)), 2) combn31 <- combn(rep((4:8 - 1) * p, each = length(26:30 - 1)) + rep(26:30 - 1, times = length(4:8)), 2) combn32 <- combn(rep((4:8 - 1) * p, each = length(46:50 - 1)) + rep(46:50 - 1, times = length(4:8)), 2) combn41 <- combn(rep((3:5 - 1) * p, each = length(36:40 - 1)) + rep(36:40 - 1, times = length(3:5)), 2) combn42 <- combn(rep((3:5 - 1) * p, each = length(46:50 - 1)) + rep(46:50 - 1, times = length(3:5)), 2) combn51 <- combn(rep((9:10 - 1) * p, each = length(36:40 - 1)) + rep(36:40 - 1, times = length(9:10)), 2) combn52 <- combn(rep((9:10 - 1) * p, each = length(46:50 - 1)) + rep(46:50 - 1, times = length(9:10)), 2) Gmrf <- rbind(t(combn11), t(combn31), t(combn32), t(combn41), t(combn42), t(combn51), t(combn52)) ## get for every correlated bunch in the decomposable model, if (graph_pattern < 4) { # a different set of predictors for (i in 1:length(Prime)) { gamma[6:10 + (i + 6) * 10, Prime[[i]]] <- 1 } ## for each Prime component ## for every Residual instead for (i in 1:length(Res)) { gamma[6:10 + (i + 10) * 10, Res[[i]]] <- 1 } } else { for (i in 1:length(Prime)) { gamma[6:10 + (i + 4) * 10, Prime[[i]]] <- 1 } ## for each Prime component ## for every Residual instead for (i in 1:length(Res)) { gamma[6:10 + (i + 9) * 10, Res[[i]]] <- 1 } } #### Sample the betas sd_b <- 1 b <- matrix(rnorm((p + 1) * s, 0, sd_b), p + 1, s) xb <- matrix(NA, n, s) for (i in 1:s) { if (sum(gamma[, i]) > 1) { xb[, i] <- x[, gamma[, i] == 1] %*% b[gamma[, i] == 1, i] } else { xb[, i] <- rep(1, n) * b[1, i] } } ## Sample the variance v_r <- mean(diag(var(xb))) / snr nu <- s + 1 M <- matrix(corr_param, s, s) diag(M) <- rep(1, s) P <- BDgraph::rgwish(n = 1, adj = G, b = 3, D = v_r * M) var <- solve(P) factor <- 10 factor_min <- 0.01 factor_max <- 1000 count <- 0 maxit <- 10000 factor_prev <- 1 repeat{ var <- var / factor * factor_prev ### Sample the errors and the Ys cVar <- chol(as.matrix(var)) # err = matrix(rnorm(n*s),n,s) %*% cVar err <- matrix(rnorm(n * s, sd = 0.5), n, s) %*% cVar y <- xb + err ## Reparametrisation ( assuming PEO is 1:s ) cVar <- t(cVar) # make it lower-tri S <- diag(diag(cVar)) sigma <- S * S L <- cVar %*% solve(S) rho <- diag(s) - solve(L) ### S/N Ratio emp_snr <- mean(diag(var(xb) %*% solve(sigma))) emp_g_snr <- mean(diag(var((err) %*% t(rho)) %*% solve(sigma))) ############## if (abs(emp_snr - snr) < (snr / 10) | count > maxit) { break } else { if (emp_snr < snr) { # increase factor factor_min <- factor } else { # decrease factor factor_max <- factor } factor_prev <- factor factor <- (factor_min + factor_max) / 2 } count <- count + 1 } ################# colnames(y) <- paste("GEX", 1:ncol(y), sep = "") colnames(G) <- colnames(y) Gy <- G gamma <- gamma[-1, ] mrfG <- Gmrf[!duplicated(Gmrf), ] data <- cbind(y, x[, -1]) # leave out the intercept because is coded inside already exampleEQTL <- list(data = data, blockList = list(1:s, s + 1:p)) ## Write data file to the user's directory by save() ## End(Not run) ```

exampleEQTL dataset