################################################### ### chunk number 1: Roptions ################################################### ## Some options for R options(width=60) options(continue=" ") if (require(Cairo)) { png <- CairoPNG } dir.create("graphics") ################################################### ### chunk number 2: toymatrix ################################################### library(isa2) library(Matrix) data <- data2 <- isa.in.silico(20,20,2,10,10)[[1]] perm1 <- sample(seq_len(nrow(data))) data2 <- data2[perm1,] perm2 <- sample(seq_len(ncol(data))) data2 <- data2[,perm2] alldata <- structure(c(data2, data), dim=c(nrow(data),ncol(data),2), dimnames=list(NULL,NULL,c("Original matrix", "Reordered matrix"))) rx <- range(alldata, finite = TRUE) nn <- 100 n0 <- min(nn, max(0, round((0 - rx[1])/(rx[2] - rx[1]) * nn))) col.regions <- c(colorRampPalette(c("blue3", "gray80"))(n0), colorRampPalette(c("gray75", "red3"))(nn - n0)) lp <- levelplot(alldata, xlab="", ylab="", scale=list(y=list(draw=FALSE), x=list(draw=FALSE)), col.regions=col.regions, between=list(x=0.5)) print(lp) ################################################### ### chunk number 3: insilico ################################################### set.seed(10) library(isa2) data <- isa.in.silico(200,200,2,50,50)[[1]] ################################################### ### chunk number 4: blockplot.pre ################################################### png("graphics/blockplot.png", 600, 400) ################################################### ### chunk number 5: blockplot ################################################### images(list(data), xlab="", ylab="", strip=FALSE) ################################################### ### chunk number 6: blockplot.post ################################################### dev.off() ################################################### ### chunk number 7: isa1 ################################################### modules <- isa(data) names(modules) ################################################### ### chunk number 8: no.mods. ################################################### ncol(modules$rows) ################################################### ### chunk number 9: simpleplot eval=FALSE ################################################### ## layout(cbind(seq(ncol(modules$rows)))) ## sapply(seq(ncol(modules$rows)), function(x) { ## par(mar=c(3,5,1,1)) ## plot(modules$rows[,x], ylim=c(-1,1), ## ylab="scores", xlab=NA, ## cex.lab=2, cex.axis=1.5) ## }) ################################################### ### chunk number 10: simpleplot2 ################################################### layout(cbind(seq(ncol(modules$rows)))) sapply(seq(ncol(modules$rows)), function(x) { par(mar=c(3,5,1,1)) plot(modules$rows[,x], ylim=c(-1,1), ylab="scores", xlab=NA, cex.lab=2, cex.axis=1.5) }) ################################################### ### chunk number 11: isaresult.rundata ################################################### names(modules$rundata) ################################################### ### chunk number 12: seedcols ################################################### colnames(modules$seeddata) ################################################### ### chunk number 13: seeddata ################################################### modules$seeddata ################################################### ### chunk number 14: mymodules ################################################### mymodules <- lapply(seq(ncol(modules$rows)), function(x) { list(rows=which(modules$rows[,x] != 0), columns=which(modules$columns[,x] != 0)) }) length(mymodules) mymodules[[1]]$rows mymodules[[1]]$columns ################################################### ### chunk number 15: insilico10 ################################################### data <- isa.in.silico(200, 100, 10) ################################################### ### chunk number 16: tenmods.pre ################################################### png("graphics/tenmods.png", 600, 400) ################################################### ### chunk number 17: tenmods ################################################### images(list(data[[1]]), strip=FALSE, xlab="", ylab="") ################################################### ### chunk number 18: tenmods.post ################################################### dev.off() ################################################### ### chunk number 19: normalize ################################################### normed.data <- isa.normalize(data[[1]]) names(normed.data) dim(normed.data$Er) dim(normed.data$Ec) ################################################### ### chunk number 20: checkscalecenter ################################################### max(abs(rowSums(normed.data$Er))) max(abs(rowSums(normed.data$Ec))) max(abs(apply(normed.data$Er, 1, sd)-1)) max(abs(apply(normed.data$Ec, 1, sd)-1)) ################################################### ### chunk number 21: generate.seeds ################################################### row.seeds <- generate.seeds(length=nrow(data[[1]]), count=1000, sparsity=c(2,5,10,100)) ################################################### ### chunk number 22: isa.iterate ################################################### modules <- isa.iterate(normed.data, row.seeds=row.seeds, thr.row=2, thr.col=2, direction="up") ################################################### ### chunk number 23: isa.iterate.res ################################################### ncol(modules$rows) ################################################### ### chunk number 24: all.zero.modules ################################################### sum(apply(modules$rows==0, 2, all)) ################################################### ### chunk number 25: unique ################################################### modules2 <- isa.unique(normed.data, modules, cor.limit=0.9) ################################################### ### chunk number 26: isa.unique.res ################################################### ncol(modules2$rows) ################################################### ### chunk number 27: correctness ################################################### found.rows <- cor(data[[2]], modules2$rows!=0) found.cols <- cor(data[[3]], modules2$columns!=0) found <- pmin(found.rows,found.cols) apply(found, 1, max) ################################################### ### chunk number 28: robustness ################################################### rob <- robustness(normed.data, modules2$rows, modules2$columns) summary(rob) ################################################### ### chunk number 29: robbox ################################################### par(cex.lab=1.5, cex.axis=1.5) boxplot(rob, ylab="Robustness") ################################################### ### chunk number 30: robbox2 ################################################### par(cex.lab=1.5, cex.axis=1.5) boxplot(rob, ylab="Robustness") ################################################### ### chunk number 31: filter.robust ################################################### modules3 <- isa.filter.robust(data[[1]], normed.data, modules2, perms=1, row.seeds=row.seeds) ################################################### ### chunk number 32: robbox3 ################################################### plot(modules3$seeddata$rob, ylab="Robustness", cex.lab=1.5) ################################################### ### chunk number 33: robbox4 ################################################### plot(modules3$seeddata$rob, ylab="Robustness", cex.lab=1.5) ################################################### ### chunk number 34: bestmods ################################################### bestmods <- order(modules3$seeddata$rob, decreasing=TRUE)[1:10] mod.cor <- pmin(cor(modules3$rows[,bestmods]!=0, data[[2]]), cor(modules3$columns[,bestmods]!=0, data[[3]])) apply(mod.cor, 1, max) apply(mod.cor, 2, max) ################################################### ### chunk number 35: othermods.pre ################################################### png("graphics/othermods.png", 400, 400) ################################################### ### chunk number 36: othermods ################################################### othermods <- order(modules3$seeddata$rob)[1:5] print(plotModules(modules3, othermods, data=data[[1]])) ################################################### ### chunk number 37: othermods.post ################################################### dev.off() ################################################### ### chunk number 38: biclust ################################################### library(biclust) five.mods <- isa.in.silico(200,200,5,20,20) modules <- isa(five.mods[[1]],thr.row=2, thr.col=2) Bc <- isa.biclust(modules) Bc ################################################### ### chunk number 39: drawHeatmap.pre ################################################### png("graphics/drawheatmap.png", 400,400) ################################################### ### chunk number 40: drawHeatmap ################################################### drawHeatmap(five.mods[[1]], Bc, number=1,local=FALSE) ################################################### ### chunk number 41: drawHeatmap.post ################################################### dev.off() ################################################### ### chunk number 42: parallel.pre ################################################### parallelCoordinates(five.mods[[1]], Bc, number=1, plotBoth=TRUE) ################################################### ### chunk number 43: parallel ################################################### parallelCoordinates(five.mods[[1]], Bc, number=1, plotBoth=TRUE) ################################################### ### chunk number 44: plotclust.pre ################################################### plotclust(Bc, five.mods[[1]], Titel="") ################################################### ### chunk number 45: plotclust ################################################### plotclust(Bc, five.mods[[1]], Titel="") ################################################### ### chunk number 46: noisegen ################################################### noise <- seq(0.1,1,by=0.1) data <- lapply(noise, function(x) isa.in.silico(500,200,10,noise=x)) ################################################### ### chunk number 47: noiseplot.pre ################################################### png("graphics/noise.png", width=600, height=400) ################################################### ### chunk number 48: noiseplot ################################################### images(lapply(data, "[[", 1), names=as.character(noise), xlab="", ylab="") ################################################### ### chunk number 49: noiseplot.post ################################################### dev.off() ################################################### ### chunk number 50: noise.isa ################################################### modules <- lapply(data, function(x) isa(x[[1]])) ################################################### ### chunk number 51: sensitivity.noise ################################################### best <- lapply(seq_along(modules), function(i) { cc <- pmin(cor(modules[[i]]$rows!=0, data[[i]][[2]]), cor(modules[[i]]$columns!=0, data[[i]][[3]])) apply(cc, 2, max) }) best.mean <- sapply(best, mean) ################################################### ### chunk number 52: sensbox ################################################### boxplot(best, names=noise, xlab="Noise level", ylab="Sensitivity", cex.lab=1.5) ################################################### ### chunk number 53: sensbox2 ################################################### boxplot(best, names=noise, xlab="Noise level", ylab="Sensitivity", cex.lab=1.5) ################################################### ### chunk number 54: spec ################################################### spec <- lapply(seq_along(modules), function(i) { cc <- pmin(cor(modules[[i]]$rows != 0, data[[i]][[2]]), cor(modules[[i]]$columns != 0, data[[i]][[3]])) apply(cc, 1, max) }) ################################################### ### chunk number 55: specbox ################################################### boxplot(spec, names=noise, xlab="Noise level", ylab="Specificity", cex.lab=1.5) ################################################### ### chunk number 56: specbox2 ################################################### boxplot(spec, names=noise, xlab="Noise level", ylab="Specificity", cex.lab=1.5) ################################################### ### chunk number 57: overlap ################################################### set.seed(1) two.over <- isa.in.silico(100,100,2,40,40,overlap.row=10) ################################################### ### chunk number 58: overlap.isa ################################################### ov.normed <- isa.normalize(two.over[[1]]) ov.seeds <- generate.seeds(count=100, length=100) ov.modules.1 <- isa.iterate(ov.normed, ov.seeds, convergence="cor", thr.row=0, thr.col=0, direction="up") ov.modules.1 <- isa.unique(ov.normed, ov.modules.1) ov.modules.1 <- isa.filter.robust(two.over[[1]], ov.normed, ov.modules.1) ov.modules.2 <- isa.iterate(ov.normed, ov.seeds, convergence="cor", thr.row=1, thr.col=1, direction="up") ov.modules.2 <- isa.unique(ov.normed, ov.modules.2) ov.modules.2 <- isa.filter.robust(two.over[[1]], ov.normed, ov.modules.2) ncol(ov.modules.1$rows) ncol(ov.modules.2$rows) ################################################### ### chunk number 59: overlapplot.pre ################################################### png("graphics/overlapplot.png", 600, 400) ################################################### ### chunk number 60: overlapplot ################################################### no.modules <- ncol(ov.modules.1$rows) + ncol(ov.modules.2$rows) plotModules(data=two.over[[1]], list(rows=cbind(ov.modules.1$rows, ov.modules.2$rows), columns=cbind(ov.modules.1$columns, ov.modules.2$column)), names=c("Input matrix", paste("Module", seq_len(no.modules)))) ################################################### ### chunk number 61: overlapplot.post ################################################### dev.off() ################################################### ### chunk number 62: sessioninfo ################################################### toLatex(sessionInfo())