inst/doc/qpgraphSimulate.R

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

### R code from vignette source 'qpgraphSimulate.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: qpgraphSimulate.Rnw:118-124
###################################################
library(qpgraph)

args(erGraphParam)
args(erMarkedGraphParam)
args(dRegularGraphParam)
args(dRegularMarkedGraphParam)


###################################################
### code chunk number 4: qpgraphSimulate.Rnw:129-133
###################################################
erGraphParam()
erMarkedGraphParam()
dRegularGraphParam()
dRegularMarkedGraphParam()


###################################################
### code chunk number 5: qpgraphSimulate.Rnw:138-140
###################################################
showClass("graphParam")
showClass("markedGraphParam")


###################################################
### code chunk number 6: qpgraphSimulate.Rnw:151-152
###################################################
args(rgraphBAM)


###################################################
### code chunk number 7: qpgraphSimulate.Rnw:158-160
###################################################
rgraphBAM(erGraphParam())
rgraphBAM(n=2, dRegularGraphParam())


###################################################
### code chunk number 8: 3regulargraph
###################################################
set.seed(1234)
g <- rgraphBAM(dRegularMarkedGraphParam(pI=2, pY=10, d=3))
plot(g)


###################################################
### code chunk number 9: qpgraphSimulate.Rnw:198-199
###################################################
args(rUGgmm)


###################################################
### code chunk number 10: qpgraphSimulate.Rnw:210-219
###################################################
rUGgmm(dRegularGraphParam(p=4, d=2))
rUGgmm(matrix(c(0, 1, 0, 1,
                1, 0, 1, 0,
                0, 1, 0, 1,
                1, 0, 1, 0), nrow=4, byrow=TRUE))
rUGgmm(matrix(c(1, 2,
                2, 3,
                3, 4,
                4, 1), ncol=2, byrow=TRUE))


###################################################
### code chunk number 11: qpgraphSimulate.Rnw:226-229
###################################################
set.seed(12345)
gmm <- rUGgmm(dRegularGraphParam(p=4, d=2))
summary(gmm)


###################################################
### code chunk number 12: qpgraphSimulate.Rnw:233-240
###################################################
class(gmm)
names(gmm)
gmm$X
gmm$p
gmm$g
gmm$mean
gmm$sigma


###################################################
### code chunk number 13: 4cyclegraph
###################################################
plot(gmm)
round(solve(gmm$sigma), digits=1)


###################################################
### code chunk number 14: qpgraphSimulate.Rnw:275-276
###################################################
rmvnorm(10, gmm)


###################################################
### code chunk number 15: qpgraphSimulate.Rnw:281-285
###################################################
set.seed(123)
X <- rmvnorm(10000, gmm)
round(solve(cov(X)), digits=1)
round(solve(gmm$sigma), digits=1)


###################################################
### code chunk number 16: qpgraphSimulate.Rnw:349-350
###################################################
args(rHMgmm)


###################################################
### code chunk number 17: qpgraphSimulate.Rnw:359-368
###################################################
rHMgmm(dRegularMarkedGraphParam(pI=1, pY=3, d=2))
rHMgmm(matrix(c(0, 1, 0, 1,
                1, 0, 1, 0,
                0, 1, 0, 1,
                1, 0, 1, 0), nrow=4, byrow=TRUE), I=1)
rHMgmm(matrix(c(1, 2,
                2, 3,
                3, 4,
                4, 1), ncol=2, byrow=TRUE), I=1)


###################################################
### code chunk number 18: qpgraphSimulate.Rnw:378-381
###################################################
set.seed(12345)
gmm <- rHMgmm(dRegularMarkedGraphParam(pI=1, pY=3, d=2))
summary(gmm)


###################################################
### code chunk number 19: qpgraphSimulate.Rnw:386-399
###################################################
class(gmm)
names(gmm)
gmm$X
gmm$I
gmm$Y
gmm$p
gmm$pI
gmm$pY
gmm$g
gmm$mean()
gmm$sigma
gmm$a
gmm$eta2


###################################################
### code chunk number 20: hmgmm
###################################################
plot(gmm)


###################################################
### code chunk number 21: qpgraphSimulate.Rnw:425-426
###################################################
rcmvnorm(10, gmm)


###################################################
### code chunk number 22: qpgraphSimulate.Rnw:431-437
###################################################
set.seed(123)
X <- rcmvnorm(10000, gmm)
csigma <- (1/10000)*sum(X[, gmm$I] == 1)*cov(X[X[, gmm$I]==1, gmm$Y]) +
          (1/10000)*sum(X[, gmm$I] == 2)*cov(X[X[, gmm$I]==2, gmm$Y])
round(solve(csigma), digits=1)
round(solve(gmm$sigma), digits=1)


###################################################
### code chunk number 23: qpgraphSimulate.Rnw:442-447
###################################################
smean <- apply(X[, gmm$Y], 2, function(x, i) tapply(x, i, mean), X[, gmm$I])
smean
gmm$mean()
abs(smean[1, ] - smean[2, ])
gmm$a


###################################################
### code chunk number 24: qpgraphSimulate.Rnw:471-473
###################################################
eQTLcrossParam()
args(eQTLcrossParam)


###################################################
### code chunk number 25: qpgraphSimulate.Rnw:479-480
###################################################
reQTLcross(n=2, eQTLcrossParam())


###################################################
### code chunk number 26: qpgraphSimulate.Rnw:497-506
###################################################
detach("package:qpgraph")
library(qtl)
library(qpgraph)

map <- sim.map(len=rep(100, times=20),
               n.mar=rep(10, times=20),
               anchor.tel=FALSE,
               eq.spacing=FALSE,
               include.x=TRUE)


###################################################
### code chunk number 27: qpgraphSimulate.Rnw:511-514
###################################################
eqtl <- reQTLcross(eQTLcrossParam(map=map, genes=50))
class(eqtl)
eqtl


###################################################
### code chunk number 28: qpgraphSimulate.Rnw:521-522
###################################################
eqtl$model


###################################################
### code chunk number 29: geneticmap
###################################################
par(mfrow=c(1,2))
plot(map)
plot(eqtl, main="eQTL model with cis-associations only")


###################################################
### code chunk number 30: qpgraphSimulate.Rnw:546-549
###################################################
set.seed(123)
eqtl <- reQTLcross(eQTLcrossParam(map=map, genes=50,
                                  cis=0.5, trans=c(5, 5)), a=5)


###################################################
### code chunk number 31: qpgraphSimulate.Rnw:559-561
###################################################
head(ciseQTL(eqtl))
transeQTL(eqtl)


###################################################
### code chunk number 32: transeqtl
###################################################
plot(eqtl, main="eQTL model with trans-eQTL")


###################################################
### code chunk number 33: qpgraphSimulate.Rnw:581-584
###################################################
set.seed(123)
cross <- sim.cross(map, eqtl)
cross


###################################################
### code chunk number 34: qpgraphSimulate.Rnw:598-603
###################################################
allcis <- ciseQTL(eqtl)
allcis[allcis$chrom==1, ]
gene <- allcis[2, "gene"]
chrom <- allcis[2, "chrom"]
location <- allcis[2, "location"]


###################################################
### code chunk number 35: qpgraphSimulate.Rnw:606-609
###################################################
connectedGenes <- graph::inEdges(gene, eqtl$model$g)[[1]]
connectedGenes <- connectedGenes[connectedGenes %in% eqtl$model$Y]
connectedGenes


###################################################
### code chunk number 36: qpgraphSimulate.Rnw:616-617
###################################################
out.mr <- scanone(cross, method="mr", pheno.col=c(gene, connectedGenes))


###################################################
### code chunk number 37: lodprofiles
###################################################
plot(out.mr, chr=chrom, ylab="LOD score", lodcolumn=1:3, col=c("black", "blue", "red"), lwd=2)
abline(v=allcis[allcis$gene == gene, "location"])
legend("topleft", names(out.mr)[3:5], col=c("black", "blue", "red"), lwd=2, inset=0.05)


###################################################
### code chunk number 38: qpgraphSimulate.Rnw:639-642
###################################################
out.perm <- scanone(cross, method="mr", pheno.col=c(gene, connectedGenes),
                    n.perm=100, verbose=FALSE)
summary(out.perm, alpha=c(0.05, 0.10))


###################################################
### code chunk number 39: qpgraphSimulate.Rnw:647-648
###################################################
summary(out.mr, perms=out.perm, alpha=0.05)


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