### 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())