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

###################################################
### code chunk number 1: setup
###################################################
pdf.options(useDingbats = FALSE)
options(width = 80)
rm(list=ls())


###################################################
### code chunk number 2: motivatingExample
###################################################
library(Biobase)
library(Rgraphviz)
library(qpgraph)
library(org.EcK12.eg.db)


###################################################
### code chunk number 3: motivatingExample
###################################################
data(EcoliOxygen)
ls()
gds680.eset


###################################################
### code chunk number 4: motivatingExample
###################################################
subset.gds680.eset
dim(subset.filtered.regulon6.1)


###################################################
### code chunk number 5: motivatingExample
###################################################
IQRs <- esApply(subset.gds680.eset, 1, IQR)
top100IQRgenes <- names(sort(IQRs, decreasing=TRUE))[1:100]
eset100 <- subset.gds680.eset[top100IQRgenes, ]
eset100


###################################################
### code chunk number 6: motivatingExample
###################################################
maskTF <- subset.filtered.regulon6.1$EgID_TF %in% top100IQRgenes
maskTG <- subset.filtered.regulon6.1$EgID_TG %in% top100IQRgenes
regulon100 <- subset.filtered.regulon6.1[maskTF & maskTG, ]
dim(regulon100)


###################################################
### code chunk number 7: motivatingExample
###################################################
pcc.est <- cor(t(exprs(eset100)))
dim(pcc.est)
pcc.est[1:5,1:5]


###################################################
### code chunk number 8: motivatingExample
###################################################
pcc.g <- qpAnyGraph(abs(pcc.est), threshold=0.75, return.type="graphNEL")
pcc.g


###################################################
### code chunk number 9: pccnet
###################################################
qpPlotNetwork(pcc.g, annotation="org.EcK12.eg.db")


###################################################
### code chunk number 10: motivatingExample
###################################################
TGgenes <- setdiff(featureNames(eset100), regulon100[, "EgID_TF"])
allTGpairs <- as.matrix(expand.grid(list(TGgenes, TGgenes)))
pcc.est[allTGpairs] <- NA


###################################################
### code chunk number 11: motivatingExample
###################################################
pcc.g <- qpAnyGraph(abs(pcc.est), threshold=0.75, return.type="graphNEL")
pcc.g


###################################################
### code chunk number 12: pccnetonlytftg
###################################################
qpPlotNetwork(pcc.g, annotation="org.EcK12.eg.db")


###################################################
### code chunk number 13: motivatingExample
###################################################
set.seed(123)
rnd.est <- matrix(runif(nrow(pcc.est)^2, min=-1, max=1),
                  nrow=nrow(pcc.est), dimnames=dimnames(pcc.est))
rnd.est <- (rnd.est + t(rnd.est)) / 2
dim(rnd.est)
rnd.est[1:5,1:5]


###################################################
### code chunk number 14: motivatingExample
###################################################
rnd.est[allTGpairs] <- NA
rnd.g <- qpAnyGraph(abs(rnd.est), threshold=0.75, return.type="graphNEL")
rnd.g


###################################################
### code chunk number 15: rndnetonlytftg
###################################################
qpPlotNetwork(rnd.g, annotation="org.EcK12.eg.db")


###################################################
### code chunk number 16: motivatingExample
###################################################
pcc.pr <- qpPrecisionRecall(abs(pcc.est), refGraph=regulon100[, 3:4])
rnd.pr <- qpPrecisionRecall(abs(rnd.est), refGraph=regulon100[, 3:4])


###################################################
### code chunk number 17: prerecpccrnd
###################################################
plot(pcc.pr, type="b", lwd=2, pch=25, xlab="Recall", ylab="Precision", bg="black")
lines(rnd.pr, type="l", lwd=2, lty=2, bg="black")
legend("topright", c("PCC", "Random"), lwd=2, pch=c(25, -1),
       lty=c(1, 2), pt.bg="black", inset=0.01)
TFgenes <- setdiff(featureNames(eset100), TGgenes)
nrow(regulon100)
totalPossibleEdges <- length(TFgenes) * length(TGgenes) + choose(length(TFgenes), 2)
totalPossibleEdges
abline(h=nrow(regulon100) / totalPossibleEdges, lwd=2, lty=3)


###################################################
### code chunk number 18: XYZnet
###################################################
gXYZ <- new("graphNEL", nodes=c("X", "Y", "Z"))
edgemode(gXYZ) <- "directed"
gXYZ <- addEdge(from=c("X", "Y"), to=c("Y", "Z"), gXYZ)
gXYZ <- layoutGraph(gXYZ, layoutType="circo")
graph.par(list(edges=list(lwd=c("X"=4, "Y"=4, "Z"=4)), nodes=list(fontsize=6, lwd=4)))
renderGraph(gXYZ)


###################################################
### code chunk number 19: undirectedGMMs
###################################################
set.seed(123)
X <- rnorm(100)
Y <- X * 2 + rnorm(100)
Z <- Y * 2 + rnorm(100)


###################################################
### code chunk number 20: pccYxZ
###################################################
cov2cor(cov(cbind(X, Y, Z)))
cor(X, Z)
plot(X, Z)


###################################################
### code chunk number 21: undirectedGMMs
###################################################
-cov2cor(solve(cov(cbind(X,Y,Z))))


###################################################
### code chunk number 22: pacXxZ
###################################################
fitX <- lm(X ~ Y)
fitZ <- lm(Z ~ Y)
plot(resid(fitX), resid(fitZ),
     xlab="Residuals X ~ Y",
     ylab="Residuals Z ~ Y")
cor.test(resid(fitX), resid(fitZ))


###################################################
### code chunk number 23: undirectedGMMs
###################################################
G <- matrix(c(0,1,0,0,0,1,0,1,1,0,0,1,0,0,1,0,1,0,0,1,0,0,1,1,0),
            nrow=5, dimnames=list(1:5, 1:5))
G
set.seed(123)
Sigma <- qpG2Sigma(G, rho=0.5)
round(solve(Sigma), digits=2)


###################################################
### code chunk number 24: undirectedGGMs
###################################################
mean(cov2cor(as.matrix(Sigma))[upper.tri(as.matrix(Sigma)) & G])


###################################################
### code chunk number 25: undirectedGMMs
###################################################
library(mvtnorm)
set.seed(123)
X <- rmvnorm(n=50000, sigma=as.matrix(Sigma))
dim(X)
S <- cov(X)
round(solve(S), digits=2)
mean(cov2cor(S)[upper.tri(S) & G])


###################################################
### code chunk number 26: undirectedGMMs
###################################################
set.seed(123)
X <- rmvnorm(n=50, sigma=as.matrix(Sigma))
dim(X)
S <- cov(X)
round(solve(S), digits=2)


###################################################
### code chunk number 27: undirectedGMMs
###################################################
qpCItest(X, i=1, j=5, Q=2:4, long.dim.are.variables=FALSE)
coef(lm(X[,1] ~ X[,2]+X[,3]+X[,4]+X[,5]))


###################################################
### code chunk number 28: undirectedGMMs
###################################################
set.seed(123)
X <- rmvnorm(n=4, sigma=as.matrix(Sigma))
dim(X)
S <- cov(X)


###################################################
### code chunk number 29: undirectedGMMs (eval = FALSE)
###################################################
## round(solve(S), digits=2)


###################################################
### code chunk number 30: undirectedGMMs
###################################################
x <- try(round(solve(S), digits=2), TRUE)
write(x, "")


###################################################
### code chunk number 31: undirectedGMMs
###################################################
qr(S)$rank


###################################################
### code chunk number 32: undirectedGMMs
###################################################
round(qpPAC(X, G, verbose=FALSE)$R, digits=2)


###################################################
### code chunk number 33: thenrr
###################################################
nrr.q5 <- qpNrr(eset100, q=5, pairup.i=TFgenes, pairup.j=featureNames(eset100),
                verbose=FALSE)
nrr.q15 <- qpNrr(eset100, q=15, pairup.i=TFgenes, pairup.j=featureNames(eset100),
                 verbose=FALSE)


###################################################
### code chunk number 34: nrrboxplot
###################################################
mask <- matrix(FALSE, nrow=100, ncol=100, dimnames=dimnames(nrr.q15))
mask[as.matrix(regulon100[, 3:4])] <- TRUE
par(mfrow=c(1, 2))
boxplot(list(Present=as.matrix(nrr.q5)[mask], Absent=as.matrix(nrr.q5)[!mask]),
        main="q=5", col="grey", ylab="NRR")
boxplot(list(Present=as.matrix(nrr.q15)[mask], Absent=as.matrix(nrr.q15)[!mask]),
        main="q=15", col="grey", ylab="NRR")


###################################################
### code chunk number 35: prerecpccrndnrr
###################################################
nrr.q5.pr <- qpPrecisionRecall(nrr.q5, refGraph=regulon100[, 3:4], decreasing=FALSE)
nrr.q15.pr <- qpPrecisionRecall(nrr.q15, refGraph=regulon100[, 3:4], decreasing=FALSE)
par(mfrow=c(1, 1))
plot(pcc.pr, type="b", lwd=2, pch=25, xlab="Recall", ylab="Precision", bg="black")
lines(rnd.pr, type="l", lwd=2, lty=2, bg="black")
lines(nrr.q5.pr, type="b", lwd=2, pch=24, bg="black")
lines(nrr.q15.pr, type="b", lwd=2, pch=23, bg="black")
legend("topright", c("PCC", "Random", "NRR q=5", "NRR q=15"), lwd=2,
       pch=c(25, -1, 24, 23), lty=c(1, 2, 1), pt.bg="black", inset=0.01)


###################################################
### code chunk number 36: nrrnet
###################################################
nrr.thr.q5 <- qpPRscoreThreshold(nrr.q5.pr, level=0.5, recall.level=FALSE, max.score=0)
nrr.g.q5 <- qpGraph(nrr.q5, threshold=nrr.thr.q5, return.type="graphNEL")
nrr.g.q5
qpTopPairs(nrr.q5[nodes(nrr.g.q5), nodes(nrr.g.q5)],
           nrr.g.q5, annotation="org.EcK12.eg.db",
           n=10)
qpPlotNetwork(nrr.g.q5, annotation="org.EcK12.eg.db")


###################################################
### code chunk number 37: eQTLyeast
###################################################
load("qpgraphWorkshop/YeastGG.RData")
ls(pattern="sacCer3")


###################################################
### code chunk number 38: eQTLyeast
###################################################
dim(BremGGsacCer3chr3n4)
BremGGsacCer3chr3n4[1:5, 265:275]
sum(is.na(BremGGsacCer3chr3n4)) / (269*112)


###################################################
### code chunk number 39: eQTLyeast
###################################################
head(sacCer3markerPos)
head(sacCer3genePos)
sacCer3chrLen


###################################################
### code chunk number 40: eQTLyeast
###################################################
allci <- qpAllCItests(BremGGsacCer3chr3n4, I=rownames(sacCer3markerPos),
                      pairup.i=rownames(sacCer3markerPos),
                      pairup.j=rownames(sacCer3genePos), verbose=FALSE)


###################################################
### code chunk number 41: eQTLyeast
###################################################
allci$p.value[rownames(sacCer3markerPos)[1:5], rownames(sacCer3genePos)[1:5]]
allci$statistic
allci$n


###################################################
### code chunk number 42: yeasteqtlmap
###################################################
par(mar=c(4, 5, 1, 1))
sel.pairs <- qpPlotMap(allci$p.value, sacCer3markerPos, sacCer3genePos,
                       sacCer3chrLen, ylab="", cex=2, p.value=0.01)
mtext("Ordered Genes", 2, line=4)


###################################################
### code chunk number 43: eQTLyeast
###################################################
dim(sel.pairs)
head(sel.pairs, n=2)


###################################################
### code chunk number 44: eQTLyeast
###################################################
markerlinks <- sapply(split(sel.pairs[,3], sel.pairs[,1]), length)
head(sacCer3markerPos[names(sort(markerlinks, decreasing=TRUE)), ], n=10)


###################################################
### code chunk number 45: eQTLyeast
###################################################
markerhotspots <- c("6768_at_x07", "6829_at_x01", "2435_at_x00")


###################################################
### code chunk number 46: eQTLyeast
###################################################
yeastRegInt <- read.table("qpgraphWorkshop/yeastractRegTwoColTable.tsv",
                          colClasses="character")
head(yeastRegInt, n=2)


###################################################
### code chunk number 47: eQTLyeast
###################################################
library("org.Sc.sgd.db")
tfIDs <- toupper(unique(yeastRegInt[, 1]))
tfIDs <- unlist(mget(tfIDs, revmap(org.Sc.sgdGENENAME), ifnotfound=NA))
tfIDs[is.na(tfIDs)] <- names(tfIDs)[is.na(tfIDs)]
tfIDs <- tfIDs[tfIDs %in% rownames(sacCer3genePos)]
tfIDs <- intersect(tfIDs, sel.pairs$j)
length(tfIDs)


###################################################
### code chunk number 48: eQTLyeast
###################################################
nrr <- qpNrr(BremGGsacCer3chr3n4, q=20, I=rownames(sacCer3markerPos),
             pairup.i=c(markerhotspots, tfIDs), pairup.j=unique(sel.pairs$j),
             restrict.Q=rownames(sacCer3genePos), verbose=FALSE)
dim(nrr)


###################################################
### code chunk number 49: eQTLyeast
###################################################
nrr <- nrr[c(markerhotspots, tfIDs, unique(sel.pairs$j)),
           c(markerhotspots, tfIDs, unique(sel.pairs$j))]
dim(nrr)


###################################################
### code chunk number 50: eQTLyeast
###################################################
g <- qpGraph(nrr, threshold=0, return.type="graphNEL")
g


###################################################
### code chunk number 51: yeasteqtlnet
###################################################
qpPlotNetwork(g, pairup.i=c(markerhotspots, tfIDs),
              pairup.j=unique(sel.pairs$j), annotation="org.Sc.sgd.db")
chr3hs <- qpTopPairs(refGraph=g, pairup.i="2435_at_x00", pairup.j=nodes(g),
                     annotation="org.Sc.sgd.db", n=Inf)
chr3hs[, 3:4]


###################################################
### code chunk number 52: eQTLyeast
###################################################
library(rtracklayer)
chr <- c(paste0("chr", as.roman(sacCer3markerPos[unique(chr3hs$i), 1])),
         paste0("chr", as.roman(sacCer3genePos[unique(chr3hs$j), 1])))
rng <- IRanges(c(sacCer3markerPos[unique(chr3hs$i), 2],
                 sacCer3genePos[unique(chr3hs$j), 2]), width=1,
               names=c(unique(chr3hs$iSymbol), unique(chr3hs$jSymbol)))
gr <- GRangesForUCSCGenome("sacCer3", chr, rng, strand="+")
head(gr)
export(gr, "chr3hsint.bed")


###################################################
### code chunk number 53: phennet
###################################################
head(pData(eset100)[, 1:3])


###################################################
### code chunk number 54: phennet
###################################################
nrr.q5.st.gp <- qpNrr(eset100, q=5, pairup.i=c(TFgenes, "Strain"),
                      pairup.j=featureNames(eset100),
                      fix.Q="GrowthProtocol", verbose=FALSE)


###################################################
### code chunk number 55: phennet
###################################################
nrr.q5.st.gp.pr <- qpPrecisionRecall(nrr.q5.st.gp, refGraph=regulon100[, 3:4],
                                     pairup.i=TFgenes, pairup.j=featureNames(eset100),
                                     decreasing=FALSE)
nrr.thr.q5.st.gp <- qpPRscoreThreshold(nrr.q5.st.gp.pr, level=0.5,
                                       recall.level=FALSE, max.score=0)
nrr.g.q5.st.gp <- qpGraph(nrr.q5.st.gp, threshold=nrr.thr.q5.st.gp,
                          return.type="graphNEL")
nrr.g.q5.st.gp


###################################################
### code chunk number 56: phennet
###################################################
library(limma)
design <- model.matrix(~ factor(Strain) + factor(GrowthProtocol), data=eset100)
fit <- lmFit(eset100, design)
fit <- eBayes(fit)
deGenes <- topTable(fit, coef=2, p.value=0.1, n=Inf)$ID
length(deGenes)


###################################################
### code chunk number 57: phennet
###################################################
qpPlotNetwork(nrr.g.q5.st.gp, pairup.i=c(TFgenes, "Strain"),
              pairup.j=nodes(nrr.g.q5.st.gp),
              highlight=deGenes, annotation="org.EcK12.eg.db")


###################################################
### code chunk number 58: sessionInfo
###################################################
toLatex(sessionInfo())