inst/doc/eQTLnetworks.R

In qpgraph: Estimation of genetic and molecular regulatory networks from high-throughput genomics data

### R code from vignette source 'eQTLnetworks.Rnw'

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex(use.unsrturl=FALSE)


###################################################
### code chunk number 2: setup
###################################################
pdf.options(useDingbats=FALSE)
options(width=80)
rm(list=ls())
try(detach("package:mvtnorm"), silent=TRUE)
try(detach("package:qtl"), silent=TRUE)


###################################################
### code chunk number 3: eQTLnetworks.Rnw:71-74
###################################################
library(GenomeInfoDb)
library(qtl)
library(qpgraph)


###################################################
### code chunk number 4: eQTLnetworks.Rnw:79-84
###################################################
map <- sim.map(len=rep(100, times=9),
               n.mar=rep(10, times=9),
               anchor.tel=FALSE,
               eq.spacing=TRUE,
               include.x=FALSE)


###################################################
### code chunk number 5: eQTLnetworks.Rnw:96-99
###################################################
set.seed(12345)
sim.eqtl <- reQTLcross(eQTLcrossParam(map=map, genes=50, cis=0.5, trans=rep(5, 5)),
                       a=2, rho=0.5)


###################################################
### code chunk number 6: simeqtlnet
###################################################
plot(sim.eqtl, main="Simulated eQTL network")


###################################################
### code chunk number 7: eQTLnetworks.Rnw:106-109
###################################################
set.seed(12345)
cross <- sim.cross(map, sim.eqtl, n.ind=100)
cross


###################################################
### code chunk number 8: eQTLnetworks.Rnw:128-134
###################################################
annot <- data.frame(chr=as.character(sim.eqtl$genes[, "chr"]),
                    start=sim.eqtl$genes[, "location"],
                    end=sim.eqtl$genes[, "location"],
                    strand=rep("+", nrow(sim.eqtl$genes)),
                    row.names=rownames(sim.eqtl$genes),
                    stringsAsFactors=FALSE)


###################################################
### code chunk number 9: eQTLnetworks.Rnw:141-147
###################################################
pMap <- lapply(map, function(x) x * 5)
class(pMap) <- "map"
annot$start <- floor(annot$start * 5)
annot$end <- floor(annot$end * 5)
genome <- Seqinfo(seqnames=names(map), seqlengths=rep(100 * 5, nchr(pMap)),
                  NA, "simulatedGenome")


###################################################
### code chunk number 10: eQTLnetworks.Rnw:153-155
###################################################
param <- eQTLnetworkEstimationParam(cross, physicalMap=pMap,
                                    geneAnnotation=annot, genome=genome)


###################################################
### code chunk number 11: eQTLnetworks.Rnw:159-161
###################################################
eqtlnet.q0 <- eQTLnetworkEstimate(param, ~ marker + gene, verbose=FALSE)
eqtlnet.q0


###################################################
### code chunk number 12: eQTLnetworks.Rnw:166-169
###################################################
eqtlnet.q0.fdr <- eQTLnetworkEstimate(param, estimate=eqtlnet.q0,
                                      p.value=0.05, method="fdr")
eqtlnet.q0.fdr


###################################################
### code chunk number 13: simeqtlnetvsfdr
###################################################
par(mfrow=c(1, 2))
plot(sim.eqtl, main="Simulated eQTL network")
plot(eqtlnet.q0.fdr, main="Esiimated eQTL network")


###################################################
### code chunk number 14: eQTLnetworks.Rnw:191-194
###################################################
eqtlnet.q0.fdr.nrr <- eQTLnetworkEstimate(param, ~ marker + gene | gene(q=3),
                                          estimate=eqtlnet.q0.fdr, verbose=FALSE)
eqtlnet.q0.fdr.nrr


###################################################
### code chunk number 15: eQTLnetworks.Rnw:199-202
###################################################
eqtlnet.q0.fdr.nrr <- eQTLnetworkEstimate(param, estimate=eqtlnet.q0.fdr.nrr,
                                          epsilon=0.1)
eqtlnet.q0.fdr.nrr


###################################################
### code chunk number 16: simeqtlnetvsnrr
###################################################
par(mfrow=c(1, 2))
plot(sim.eqtl, main="Simulated eQTL network")
plot(eqtlnet.q0.fdr.nrr, main="Esiimated eQTL network")


###################################################
### code chunk number 17: eQTLnetworks.Rnw:225-227
###################################################
eqtls <- alleQTL(eqtlnet.q0.fdr.nrr)
median(sapply(split(eqtls$QTL, eqtls$gene), length))


###################################################
### code chunk number 18: eQTLnetworks.Rnw:237-240
###################################################
eqtlnet.q0.fdr.nrr.sel <- eQTLnetworkEstimate(param, estimate=eqtlnet.q0.fdr.nrr,
                                              alpha=0.05)
eqtlnet.q0.fdr.nrr.sel


###################################################
### code chunk number 19: simeqtlnetvsnrrsel
###################################################
par(mfrow=c(1, 2))
plot(sim.eqtl, main="Simulated eQTL network")
plot(eqtlnet.q0.fdr.nrr.sel, main="Esiimated eQTL network")


###################################################
### code chunk number 20: hiveplot
###################################################
library(grid)
library(graph)
grid.newpage()
pushViewport(viewport(layout=grid.layout(3, 3)))
for (i in 1:3) {
  for (j in 1:3) {
    chr <- (i-1) * 3 + j
    pushViewport(viewport(layout.pos.col=j, layout.pos.row=i))
    plot(eqtlnet.q0.fdr.nrr.sel, type="hive", chr=chr)
    grid.text(paste0("chr", as.roman(chr)), x=unit(0.05, "npc"),
              y=unit(0.9, "npc"), just="left")
    grid.text("genes", x=unit(0.08, "npc"), y=unit(0.1, "npc"), just="left", gp=gpar(cex=0.9))
    grid.text("all chr", x=unit(0.92, "npc"), y=unit(0.2, "npc"), just="right", gp=gpar(cex=0.9))
    grid.text("genes", x=unit(0.92, "npc"), y=unit(0.1, "npc"), just="right", gp=gpar(cex=0.9))
    grid.text("markers", x=unit(0.5, "npc"), y=unit(0.95, "npc"), just="centre", gp=gpar(cex=0.9))
    popViewport(2)
  }
}


###################################################
### code chunk number 21: info
###################################################
toLatex(sessionInfo())