|
1 | 1 |
|
2 | | -#' Simulate cis gene expression |
3 | | -#' |
4 | | -#' This function simulates cis gene expression based on genotype data. |
5 | | -#' It generates gene expression matrix E using genotype matrix G and adjacency matrix A. |
6 | | -#' |
7 | | -#' @param G A matrix of genotypes with dimensions n x p. |
8 | | -#' @param A A binary adjacency matrix of dimensions p x m indicating direct effects. |
9 | | -#' @param phi_v The per-SNP heritability value. |
10 | | -#' |
11 | | -#' @return A matrix E of dimensions n x m representing the simulated gene expression. |
12 | | -#' |
13 | | -#' @examples |
14 | | -#' n <- 1000 |
15 | | -#' p <- 40 |
16 | | -#' m <- 8 |
17 | | -#' G <- matrix(rbinom(n * p, size = 2, prob = 0.3), ncol = p) |
18 | | -#' A <- matrix(sample(0:1, p * m, replace = TRUE), nrow = p) |
19 | | -#' phi_v <- 0.05 |
20 | | -#' E <- simulate_cis_expression(G, A, phi_v) |
21 | | -#' |
22 | | -#' @export |
23 | | -simulate_cis_expression <- function(G, A, phi_v) { |
24 | | - n <- nrow(G) |
25 | | - m <- ncol(A) |
26 | | - E <- matrix(NA, n, m) |
27 | | - |
28 | | - for (j in 1:m) { |
29 | | - connected_snps <- which(A[, j] == 1) |
30 | | - num_snps <- length(connected_snps) |
31 | | - |
32 | | - beta <- rep(0, ncol(G)) |
33 | | - for (i in connected_snps) { |
34 | | - if (i == connected_snps[1]) { |
35 | | - beta[i] <- 1 |
36 | | - } else { |
37 | | - beta[i] <- sqrt(beta[connected_snps[1]]^2 * var(G[, connected_snps[1]]) / var(G[, i])) |
38 | | - } |
39 | | - } |
40 | | - |
41 | | - # Corrected calculation of variance_sum and sigma2_cis |
42 | | - variance_sum <- var(G[, connected_snps, drop = FALSE] %*% beta[connected_snps]) |
43 | | - sigma2_cis <- var(G[, connected_snps[1]] * beta[connected_snps[1]]) / phi_v - variance_sum |
44 | | - while (sigma2_cis <= 0) { |
45 | | - phi_v <- phi_v - 0.01 |
46 | | - sigma2_cis <- var(G[, connected_snps[1]] * beta[connected_snps[1]]) / phi_v - variance_sum |
47 | | - } |
48 | | - |
49 | | - # Simulate gene expression |
50 | | - E_tmp <- G %*% beta + rnorm(n, mean = 0, sd = sqrt(sigma2_cis)) |
51 | | - # E[, j] <- scale(E_tmp) |
52 | | - E[, j] <- E_tmp |
53 | | - } |
54 | | - |
55 | | - return(E) |
56 | | -} |
57 | | - |
58 | | - |
59 | | - |
60 | 2 | #' Simulate trans gene expression for different simulation scenarios |
61 | 3 | #' |
62 | 4 | #' This function can perform Type I error simulations or Accuracy and False Discovery Rate (ACC/FDR) simulations |
|
0 commit comments