### R code from vignette source 'DREAM4.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: init ################################################### printf <- function(...)print(noquote(sprintf(...))) options(width=60) ################################################### ### code chunk number 2: load-01-01 ################################################### library(DREAM4) data(dream4_010_01) show(dream4_010_01) ################################################### ### code chunk number 3: extractFullMatrix ################################################### mtx.all <- assays (dream4_010_01)[[1]] dim(mtx.all) mtx.all[1:10, c(1,2,126)] set.seed(37) print(colnames(mtx.all)[sort(sample(1:ncol(mtx.all), 16))]) ################################################### ### code chunk number 4: retrieveExpressionData ################################################### mtx.goldStandard <- metadata (dream4_010_01)[[1]] ################################################### ### code chunk number 5: colnames ################################################### grep("perturbation.1.", colnames(mtx.all), v=T, fixed=TRUE) ################################################### ### code chunk number 6: extractAndTranspose ################################################### ts1.columns <- grep("perturbation.1.", colnames(mtx.all), fixed=TRUE) mtx.ts1 <- t(mtx.all[, ts1.columns]) print(mtx.ts1[1:3, 1:3]) ################################################### ### code chunk number 7: extractGoldStandard ################################################### print(names(metadata(dream4_010_01))) mtx.goldStandard <- metadata(dream4_010_01)[[1]] ################################################### ### code chunk number 8: transformGS ################################################### idx <- which(mtx.goldStandard == 1) idx.m1 <- idx -1 rows <- idx.m1 %% nrow (mtx.goldStandard) + 1 cols <- idx.m1 %/% nrow (mtx.goldStandard) + 1 tbl.goldStandard <- data.frame(Regulator=rownames(mtx.goldStandard)[rows], Target=colnames(mtx.goldStandard)[cols], Source=rep('goldStandard', length(rows)), stringsAsFactors=FALSE) ################################################### ### code chunk number 9: infer0 ################################################### library(networkBMA) tbl.inferred <- networkBMA(mtx.ts1, nTimePoints=nrow(mtx.ts1), self=FALSE) tbl.inferred <- tbl.inferred[order(tbl.inferred$PostProb, decreasing=TRUE),] ################################################### ### code chunk number 10: aupr0 ################################################### tbl.contingency <- contabs.netwBMA(tbl.inferred, tbl.goldStandard[,-3]) pr <- scores(tbl.contingency, what = c("precision","recall")) colors <- c ("blue", "darkred") plot(pr$recall, pr$precision, type='b', xlab='RECALL', ylab='PRECISION', col=colors[1], xlim=c(0,1), ylim=c(0,1)) ################################################### ### code chunk number 11: printaupr ################################################### print(prc(tbl.contingency, plotit=FALSE)) ################################################### ### code chunk number 12: explicateAUPC ################################################### last.row <- nrow(tbl.contingency) mid.row <- as.integer(round(last.row)/2) print(tbl.contingency[c(1,mid.row,last.row),]) ################################################### ### code chunk number 13: headtblInferred ################################################### print(head(tbl.inferred, n=10)) ################################################### ### code chunk number 14: g10targets ################################################### print(subset(tbl.inferred, Regulator=="G10")) ################################################### ### code chunk number 15: g10inGoldStandard ################################################### mtx.goldStandard["G10",] ################################################### ### code chunk number 16: simulateG10binding ################################################### set.seed(37) tbl.priors <- matrix(data=runif(100, 0, 0.4), nrow=10, ncol=10, dimnames=list(rownames(mtx.goldStandard), colnames(mtx.goldStandard))) tbl.priors["G10", c("G3", "G4")] <- runif(2, 0.8, 1.0) ################################################### ### code chunk number 17: inferWithPriors ################################################### tbl.inferredWithPriors <- networkBMA(mtx.ts1, nTimePoints=nrow(mtx.ts1), prior.prob=tbl.priors, self=FALSE) tbl.inferredWithPriors <- tbl.inferredWithPriors[order(tbl.inferredWithPriors$PostProb, decreasing=TRUE),] print(tbl.inferredWithPriors) ################################################### ### code chunk number 18: aupr2 ################################################### tbl.contingencyWithPriors <- contabs.netwBMA(tbl.inferredWithPriors, tbl.goldStandard[,-3]) plot(pr$recall, pr$precision, type='b', xlab='RECALL', ylab='PRECISION', col=colors[1], xlim=c(0,1), ylim=c(0,1)) prWithPriors <- scores( tbl.contingencyWithPriors, what = c("precision","recall")) lines(prWithPriors$recall, prWithPriors$precision, type='b', xlab='RECALL', ylab='PRECISION', col=colors[2], xlim=c(0,1), ylim=c(0,1)) legend.titles = c("expression only", "with priors") legend (0.6, 0.9, legend.titles, colors) print(prc(tbl.contingencyWithPriors, plotit=FALSE)) ################################################### ### code chunk number 19: printx ################################################### print(tbl.contingencyWithPriors)