cis-trans mediation simulations

cis_trans_mediation.R

# Null is true
nsim <- 1000
p <- matrix(0, nsim, 3)

n <- 1000
g1 <- rbinom(n, 2, 0.2)
g2 <- rbinom(n, 2, 0.3)
g3 <- rbinom(n, 2, 0.4)
g <- cbind(g1, g2, g3)


for(i in 1:nsim)
{

	b <- rnorm(3)

	cis <- g %*% b + rnorm(n)
	trans <- rnorm(1) * cis + rnorm(n)

	gscore <- rowSums(g)
	p[i, 1] <- coefficients(summary(lm(trans ~ gscore + cis)))[2,4]

	gscore2 <- g %*% b
	p[i, 2] <- coefficients(summary(lm(trans ~ gscore2 + cis)))[2,4]

	gscore3 <- g %*% -b
	p[i, 3] <- coefficients(summary(lm(trans ~ gscore3 + cis)))[2,4]
}

hist(p[,1])
hist(p[,2])
hist(p[,3])

apply(p, 2, mean)



# Null is true for SNPs 1-3 and SNP 4 influences trans directly

nsim <- 1000
p <- matrix(0, nsim, 3)

n <- 1000
g1 <- rbinom(n, 2, 0.2)
g2 <- rbinom(n, 2, 0.3)
g3 <- rbinom(n, 2, 0.4)
g4 <- rbinom(n, 2, 0.4)
g <- cbind(g1, g2, g3, g4)

for(i in 1:nsim)
{
	cat(i, "\n")
	b <- rnorm(4)

	cis <- g[,1:3] %*% b[1:3] + rnorm(n)
	trans <- rnorm(1) * cis + rnorm(n) + b[4] * g[,4]

	gscore <- rowSums(g)
	p[i, 1] <- coefficients(summary(lm(trans ~ gscore + cis)))[2,4]

	gscore2 <- g %*% b
	p[i, 2] <- coefficients(summary(lm(trans ~ gscore2 + cis)))[2,4]

	gscore3 <- g %*% -b
	p[i, 3] <- coefficients(summary(lm(trans ~ gscore3 + cis)))[2,4]
}

hist(p[,1])
hist(p[,2])
hist(p[,3])

apply(p, 2, function(x) mean(-log10(x)))