library(qtl) library(qpgraph) library(doParallel) ## parameters controling the entire simulation experiment useMPIcluster <- TRUE ## set n.cores to an appropriate value below if this is set to TRUE n.cores <- detectCores() ## use all available cores (you might want to change this!!) if (useMPIcluster) ## number of cores to use in parallel in an MPI cluster n.cores <- 120 ## THIS NUMBER SHOULD MATCH THE SGE SCRIPT LINE #$ -pe ompi 120 n <- 100 ## number of individuals per cross data set rho <- seq(0.25, 0.75, by=0.25) ## mean marginal correlation between genes alpha <- c(0.5, 1.0, 2.5, 5.0) ## additive effects of the eQTL n.simxeffect <- 10 ## number of simulated eQTL models per effect size n.sim <- 1000 ## number of simulated data sets per simulated eQTL model n.genes <- 5 ## number of genes map <- sim.map(len=100, ## genetic map consisting of one chromosome 100 cM long n.mar=10, ## with 10 markers equally spaced along the chromosome anchor.tel=FALSE, ## and no telomeric markers eq.spacing=TRUE, include.x=FALSE) selBiasLODcutoff <- 3 ## selection bias LOD score cutoff cl <- NULL if (useMPIcluster) { cl <- makeCluster(n.cores, type="MPI") clusterEvalQ(cl, library(qtl)) clusterEvalQ(cl, library(qpgraph)) clusterSetRNGStream(cl) clusterExport(cl, list("n", "rho", "alpha", "map", "n.sim")) registerDoParallel(cl=cl) } else registerDoParallel(cores=n.cores) ## eQTL model: one eQTL and a chain of associated ## genes where the first one has the eQTL, using a ## genetic map consisting of one chromosome 100 cM long ## ensure the eQTL is located right over a marker eqtl <- eQTLcross(map) eqtl <- addGenes(eqtl, n.genes) eqtl <- addeQTL(eqtl, "g1", location=map[[1]][floor(nmar(map)/2)]) for (i in 1:(n.genes-1)) eqtl <- addGeneAssociation(eqtl, paste0("g", i), paste0("g", i+1)) ## build data structure to store the estimated additive effects ## and LOD scores from the eQTL to each different gene under model1 LODahat <- vector(mode="list", length=length(rho)) names(LODahat) <- paste0("r", rho) LODahat <- lapply(LODahat, function(x) { x <- vector(mode="list", length=length(alpha)) names(x) <- paste0("a", alpha) lapply(x, function(x) matrix(NA, nrow=n.simxeffect*n.sim, ncol=2*eqtl$model$pY, dimnames=list(NULL, c(paste0(eqtl$model$Y, "lod"), paste0(eqtl$model$Y, "ahat"))))) }) for (r in rho) { ## for each mean marginal correlation value cat("rho=", r, "\n") for (a in alpha) { ## for each additive effect value for (i in 1:n.simxeffect) { ## simulate n.simxeffect eQTL models sim.eqtl <- reQTLcross(eqtl, rho=r, a=a) if (useMPIcluster) clusterExport(cl, list("sim.eqtl")) LODahat[[paste0("r", r)]][[paste0("a", a)]][((i-1)*n.sim+1):(i*n.sim), ] <- foreach (icount(n.sim), .combine=rbind) %dopar% { cross <- sim.cross(map, sim.eqtl, n.ind=n) ## simulate eQTL data from the simulated eQTL model out <- scanone(cross, pheno.col=paste0("g", 1:n.genes), ## calculate LOD scores for each gene expression method="mr") ## at the eQTL ## get the LOD scores for every gene at the marker where the eQTL is located lod <- as.numeric(out[which(out$pos - alleQTL(sim.eqtl)$location == 0), paste0("g", 1:n.genes)]) names(lod) <- paste0("g", 1:n.genes, "lod") ahat <- apply(cross$pheno, 2, ## estimate the additive effect of the eQTL on each gene function(y, i) abs(sum(c(1, -1)*tapply(y, i, mean))), cross$qtlgeno) names(ahat) <- paste0("g", 1:n.genes, "ahat") c(lod, ahat) } } } } if (useMPIcluster) ## stop the MPI cluster, if necessary stopCluster(cl) ## plot the figure labelPanel <- function(lab, font=2, cex=2, offsetx=0.05, offsety=0.05) { par(xpd=TRUE) w <- par("usr")[2] - par("usr")[1] h <- par("usr")[4] - par("usr")[3] text(par("usr")[1]-w*offsetx, par("usr")[4]+h*offsety, lab, font=font, cex=cex) par(xpd=FALSE) } pdf(file="Figure1.pdf", height=7, width=7, useDingbats=FALSE) par(mfrow=c(2, 2), mar=c(4, 5, 2, 2)) plot(eqtl$model, layoutType="twopi", lwd=3) labelPanel(LETTERS[1], offsetx=0.2, offsety=0.08) for (i in seq(along=rho)) { rng <- 0 for (j in seq(along=alpha)) rng <- max(c(rng, colMeans(LODahat[[paste0("r", rho[i])]][[paste0("a", alpha[j])]][, (n.genes+1):(2*n.genes)]))) rng <- c(0, rng) plot(colMeans(LODahat[[paste0("r", rho[i])]][[1]][, (n.genes+1):(2*n.genes)]), lty=1, xaxt="n", panel.first=grid(lwd=2), las=1, type="n", ylim=rng, xlab="Genes", ylab="eQTL additive effect", cex.axis=1.2, cex.lab=1.5) axis(1, at=1:n.genes, labels=paste0("g", 1:n.genes), cex.axis=1.2, cex.lab=1.5) axis(4, las=1, cex.axis=1.2) for (k in seq(along=alpha)) ## plot additive effects from all simulations lines(colMeans(LODahat[[paste0("r", rho[i])]][[paste0("a", alpha[k])]][, (n.genes+1):(2*n.genes)]), type="l", lwd=3, lty=k, col=grey(0.70)) for (k in seq(along=alpha)) { ## plot additive effects under selection bias, that is, meanSelBias <- ## only those from eQTLs whose LOD score meets the minimum cutoff in 'selBiasLODcutoff' mapply(function(selBiasMask, ahat) mean(ahat[selBiasMask]), as.data.frame(LODahat[[paste0("r", rho[i])]][[paste0("a", alpha[k])]][, 1:n.genes] > selBiasLODcutoff), as.data.frame(LODahat[[paste0("r", rho[i])]][[paste0("a", alpha[k])]][, (n.genes+1):(n.genes*2)])) lines(meanSelBias, type="l", lwd=2, lty=k, col="black") } legend("topright", sprintf("a = %.1f", rev(alpha)), lty=length(alpha):1, lwd=2, seg.len=3, inset=0.025, bg="white", title=eval(substitute(expression(rho==r), list(r=rho[i])))) labelPanel(LETTERS[i+1], offsetx=0.2, offsety=0.08) } dev.off()