################################################### ### chunk number 1: set width ################################################### options(width=60) options(continue=" ") ################################################### ### chunk number 2: load the packages ################################################### library(eisa) ################################################### ### chunk number 3: load the data ################################################### library(ALL) library(hgu95av2.db) library(affy) data(ALL) ################################################### ### chunk number 4: simple ISA ################################################### thr.gene <- 2.7 thr.cond <- 1.4 set.seed(1) # to get the same results, always # modules <- ISA(ALL, thr.gene=thr.gene, thr.cond=thr.cond) data(ALLModulesSmall) modules <- ALLModulesSmall ################################################### ### chunk number 5: type in name of ISAModules object ################################################### modules ################################################### ### chunk number 6: accessors 1 ################################################### length(modules) dim(modules) ################################################### ### chunk number 7: accessors 2 ################################################### featureNames(modules)[1:5] sampleNames(modules)[1:5] ################################################### ### chunk number 8: number of features and samples ################################################### getNoFeatures(modules) getNoSamples(modules) colnames(pData(modules)) getOrganism(modules) annotation(modules) ################################################### ### chunk number 9: indexing ################################################### modules[[1:5]] ################################################### ### chunk number 10: indexing 2 ################################################### chr <- get(paste(annotation(modules), sep="", "CHR")) chr1features <- sapply(mget(featureNames(modules), chr), function(x) "1" %in% x) modules[chr1features,] ################################################### ### chunk number 11: indexing 3 ################################################### modules[ ,grep("^B", pData(modules)$BT)] ################################################### ### chunk number 12: list genes in modules ################################################### getFeatureNames(modules)[[1]] ################################################### ### chunk number 13: list conditions in modules ################################################### getSampleNames(modules)[[1]] ################################################### ### chunk number 14: query scores ################################################### getFeatureScores(modules, 3) getSampleScores(modules, 3) ################################################### ### chunk number 15: query all scores ################################################### dim(getFeatureMatrix(modules)) dim(getSampleMatrix(modules)) ################################################### ### chunk number 16: seed data ################################################### seedData(modules) ################################################### ### chunk number 17: run data ################################################### runData(modules) ################################################### ### chunk number 18: GO enrichment ################################################### GO <- ISAGO(modules) ################################################### ### chunk number 19: list GO result ################################################### GO ################################################### ### chunk number 20: GO summary ################################################### summary(GO$BP, p=0.001)[[1]][,-6] ################################################### ### chunk number 21: GO gene ids by category ################################################### geneIdsByCategory(GO$BP)[[1]][[3]] ################################################### ### chunk number 22: GO info ################################################### sigCategories(GO$BP)[[1]] library(GO.db) mget(sigCategories(GO$BP)[[1]][1:3], GOTERM) ################################################### ### chunk number 23: KEGG enrichment ################################################### KEGG <- ISAKEGG(modules) KEGG summary(KEGG, p=1e-3)[[1]] library(KEGG.db) mget(sigCategories(KEGG, p=0.001)[[1]], KEGGPATHID2NAME) ################################################### ### chunk number 24: CHR enrichment ################################################### CHR <- ISACHR(modules) summary(CHR,p=0.05)[[2]][,-6] ################################################### ### chunk number 25: CHR enrichment list ################################################### unlist(mget(geneIdsByCategory(CHR)[[2]][[1]], org.Hs.egSYMBOL)) ################################################### ### chunk number 26: miRNA enrichment ################################################### if (require(targetscan.Hs.eg.db)) { miRNA <- ISAmiRNA(modules) summary(miRNA, p=0.1)[[7]] } ################################################### ### chunk number 27: biclust ################################################### library(biclust) Bc <- as(modules, "Biclust") Bc ################################################### ### chunk number 28: heatmap eval=FALSE ################################################### ## ISA2heatmap(modules, 1, ALL, norm="sample", scale="none", ## col=heat.colors(12)) ################################################### ### chunk number 29: heatmap-cb eval=FALSE ################################################### ## sc <- function(x, ab) { ## x <- x-min(x) ## r <- range(x) ## x <- x/(r[2]-r[1]) ## x * (ab[2]-ab[1]) + ab[1] ## } ## cb <- sc(1:12, range(exprs(ALL)[getFeatureNames(modules, 1)[[1]], ## getSampleNames(modules, 2)[[1]] ])) ## par(mar=c(3,4,5,2)+0.1) ## image(rbind(cb), axes=FALSE) ## axis(2, at=sc(1:12, 0:1), cb, label=format(cb, digits=2)) ################################################### ### chunk number 30: save heatmap files ################################################### png("plot-heatmap.png", 400, 400) ISA2heatmap(modules, 1, ALL, norm="sample", scale="none", col=heat.colors(12)) dev.off() png("plot-heatmap-colorbar.png", 100, 400) sc <- function(x, ab) { x <- x-min(x) r <- range(x) x <- x/(r[2]-r[1]) x * (ab[2]-ab[1]) + ab[1] } cb <- sc(1:12, range(exprs(ALL)[getFeatureNames(modules, 1)[[1]], getSampleNames(modules, 2)[[1]] ])) par(mar=c(3,4,5,2)+0.1) image(rbind(cb), axes=FALSE) axis(2, at=sc(1:12, 0:1), cb, label=format(cb, digits=2)) dev.off() ################################################### ### chunk number 31: parallel.pre ################################################### if (require(Cairo)) { CairoPNG("plot-parallel.png", 400, 400 ) } else { png("plot-parallel.png", 400, 400 ) } ################################################### ### chunk number 32: parallel ################################################### profilePlot(modules, 2, ALL, plot="both") ################################################### ### chunk number 33: parallel.post ################################################### dev.off() ################################################### ### chunk number 34: goplot1 ################################################### goplot.2 <- gograph(summary(GO$BP, p=0.05)[[1]]) ################################################### ### chunk number 35: gograph object ################################################### summary(goplot.2) ################################################### ### chunk number 36: gograph object 2 ################################################### goplot.2$width goplot.2$height ################################################### ### chunk number 37: goplot2 eval=FALSE ################################################### ## x11(width=10, height=10*goplot.2$height/goplot.2$width) ## gographPlot(goplot.2) ################################################### ### chunk number 38: goplot2-real eval=FALSE ################################################### ## ps.options(fonts=c("serif", "mono")) ## gographPlot(goplot.2) ################################################### ### chunk number 39: goplot2-real-real ################################################### ps.options(fonts=c("serif", "mono")) gographPlot(goplot.2) ################################################### ### chunk number 40: goplot names ################################################### V(goplot.2)$abbrv[1:5] lapply(V(goplot.2)$desc[1:5], strwrap) ################################################### ### chunk number 41: condplot eval=FALSE ################################################### ## col <- ifelse( grepl("^B", pData(modules)$BT), "darkolivegreen", "orange") ## condPlot(modules, 2, ALL, col=col) ################################################### ### chunk number 42: condplot-real ################################################### col <- ifelse( grepl("^B", pData(modules)$BT), "darkolivegreen", "orange") condPlot(modules, 2, ALL, col=col) ################################################### ### chunk number 43: HTML1 ################################################### htmldir <- tempdir() ISAHTMLTable(modules=modules, target.dir=htmldir, GO=GO, KEGG=KEGG, CHR=CHR) ################################################### ### chunk number 44: HTML2 eval=FALSE ################################################### ## if (interactive()) { ## browseURL(URLencode(paste("file://", htmldir, "/maintable.html", sep=""))) ## } ################################################### ### chunk number 45: HTML modules ################################################### ISAHTMLModules(eset=ALL, modules=modules, target.dir=htmldir, GO=GO, KEGG=KEGG, CHR=CHR) ################################################### ### chunk number 46: coherence ################################################### constantVariance(exprs(ALL), Bc, number=1) additiveVariance(exprs(ALL), Bc, number=1) multiplicativeVariance(exprs(ALL), Bc, number=1) signVariance(exprs(ALL), Bc, number=1) ################################################### ### chunk number 47: coherence all ################################################### cv <- sapply(seq_len(Bc@Number), function(x) constantVariance(exprs(ALL), Bc, number=x)) av <- sapply(seq_len(Bc@Number), function(x) additiveVariance(exprs(ALL), Bc, number=x)) cor(av, cv) ################################################### ### chunk number 48: robustness ################################################### seedData(modules)$rob ################################################### ### chunk number 49: filtering ################################################### varLimit <- 0.5 kLimit <- 4 ALimit <- 5 flist <- filterfun(function(x) var(x)>varLimit, kOverA(kLimit,ALimit)) ALL.filt <- ALL[genefilter(ALL, flist), ] ################################################### ### chunk number 50: Entrez ################################################### ann <- annotation(ALL.filt) library(paste(ann, sep=".", "db"), character.only=TRUE) ENTREZ <- get( paste(ann, sep="", "ENTREZID") ) EntrezIds <- mget(featureNames(ALL.filt), ENTREZ) keep <- sapply(EntrezIds, function(x) length(x) >= 1 && !is.na(x)) ALL.filt.2 <- ALL.filt[keep,] ################################################### ### chunk number 51: Entrez unique ################################################### vari <- apply(exprs(ALL.filt.2), 1, var) larg <- findLargest(featureNames(ALL.filt.2), vari, data=annotation(ALL.filt.2)) ALL.filt.3 <- ALL.filt.2[larg,] ################################################### ### chunk number 52: normed ################################################### ALL.normed <- ISANormalize(ALL.filt.3) ls(assayData(ALL.normed)) dim(featExprs(ALL.normed)) dim(sampExprs(ALL.normed)) ################################################### ### chunk number 53: seeds ################################################### set.seed(3) random.seeds <- generate.seeds(length=nrow(ALL.normed), count=100) ################################################### ### chunk number 54: smart seeds ################################################### type <- as.character(pData(ALL.normed)$BT) ss1 <- ifelse(grepl("^B", type), -1, 1) ss2 <- ifelse(grepl("^B1", type), 1, 0) ss3 <- ifelse(grepl("^B2", type), 1, 0) ss4 <- ifelse(grepl("^B3", type), 1, 0) ss5 <- ifelse(grepl("^B4", type), 1, 0) ss6 <- ifelse(grepl("^T1", type), 1, 0) ss7 <- ifelse(grepl("^T2", type), 1, 0) ss8 <- ifelse(grepl("^T3", type), 1, 0) ss9 <- ifelse(grepl("^T4", type), 1, 0) smart.seeds <- cbind(ss1, ss2, ss3, ss4, ss5, ss6, ss7, ss8, ss9) ################################################### ### chunk number 55: iteration ################################################### modules1 <- ISAIterate(ALL.normed, feature.seeds=random.seeds, thr.feat=1.5, thr.samp=1.8, convergence="cor") modules2 <- ISAIterate(ALL.normed, sample.seeds=smart.seeds, thr.feat=1.5, thr.samp=1.8, convergence="cor") ################################################### ### chunk number 56: unique ################################################### modules1.unique <- ISAUnique(ALL.normed, modules1) modules2.unique <- ISAUnique(ALL.normed, modules2) length(modules1.unique) length(modules2.unique) ################################################### ### chunk number 57: robust ################################################### modules1.robust <- ISAFilterRobust(ALL.normed, modules1.unique) modules2.robust <- ISAFilterRobust(ALL.normed, modules2.unique) length(modules1.robust) length(modules2.robust) ################################################### ### chunk number 58: B vs T ################################################### scores1 <- getSampleMatrix(modules1.robust) tt1 <- colttests(scores1, as.factor(substr(type, 1, 1))) scores2 <- getSampleMatrix(modules2.robust) tt2 <- colttests(scores2, as.factor(substr(type, 1, 1))) sign1 <- which(p.adjust(tt1$p.value, "holm") < 0.05) sign2 <- which(p.adjust(tt2$p.value, "holm") < 0.05) sign1 sign2 ################################################### ### chunk number 59: conditions-plots-dysreg eval=FALSE ################################################### ## color <- ifelse(grepl("T", type), "orange", "darkolivegreen") ## layout(cbind(1:2)) ## condPlot(modules1.robust, which.min(tt1$p.value), ALL.normed, col=color, ## main="Best separator, random seeds") ## condPlot(modules2.robust, which.min(tt2$p.value), ALL.normed, col=color, ## main="Best separator, smart seeds") ################################################### ### chunk number 60: conditions-plots-dysreg-real ################################################### color <- ifelse(grepl("T", type), "orange", "darkolivegreen") layout(cbind(1:2)) condPlot(modules1.robust, which.min(tt1$p.value), ALL.normed, col=color, main="Best separator, random seeds") condPlot(modules2.robust, which.min(tt2$p.value), ALL.normed, col=color, main="Best separator, smart seeds") ################################################### ### chunk number 61: extract separators ################################################### modules1.TB <- modules1.robust[[sign1]] modules2.TB <- modules2.robust[[sign2]] ################################################### ### chunk number 62: dysreg correlation ################################################### cor(getSampleMatrix(modules1.TB), getSampleMatrix(modules2.TB)) cor(getFeatureMatrix(modules1.TB), getFeatureMatrix(modules2.TB)) ################################################### ### chunk number 63: dysreg enrichment ################################################### GO.dysreg1 <- ISAGO(modules1.TB) GO.dysreg2 <- ISAGO(modules2.TB) KEGG.dysreg1 <- ISAKEGG(modules1.TB) KEGG.dysreg2 <- ISAKEGG(modules2.TB) ################################################### ### chunk number 64: dysreg enriched categories ################################################### gocats <- unique(unlist(c(sigCategories(GO.dysreg1$BP), sigCategories(GO.dysreg1$CC), sigCategories(GO.dysreg1$MF), sigCategories(GO.dysreg2$BP), sigCategories(GO.dysreg2$CC), sigCategories(GO.dysreg2$MF)))) keggp <- unique(unlist(c(sigCategories(KEGG.dysreg1), sigCategories(KEGG.dysreg2)))) sapply(mget(gocats, GOTERM), Term) mget(keggp, KEGGPATHID2NAME) ################################################### ### chunk number 65: B1 vs others ################################################### keep <- grepl("^B[1234]", type) type.B <- ifelse(type[keep]=="B1", "B1", "Bx") scores.B1.1 <- scores1[keep,] scores.B1.2 <- scores2[keep,] tt.B1.1 <- colttests(scores.B1.1, as.factor(type.B)) tt.B1.2 <- colttests(scores.B1.2, as.factor(type.B)) min(p.adjust(na.omit(tt.B1.1$p.value))) min(p.adjust(na.omit(tt.B1.2$p.value))) ################################################### ### chunk number 66: B1-cond-plots eval=FALSE ################################################### ## color <- ifelse(type == "B1", "green", ## ifelse(grepl("^T", type), "orange", "darkolivegreen")) ## layout(cbind(1:2)) ## condPlot(modules1.robust, which.min(tt.B1.1$p.value), ALL.normed, col=color, ## main="Best B1 separator, random seeds") ## condPlot(modules2.robust, which.min(tt.B1.2$p.value), ALL.normed, col=color, ## main="Best B1 separator, smart seeds") ################################################### ### chunk number 67: B1-cond-plots-real ################################################### color <- ifelse(type == "B1", "green", ifelse(grepl("^T", type), "orange", "darkolivegreen")) layout(cbind(1:2)) condPlot(modules1.robust, which.min(tt.B1.1$p.value), ALL.normed, col=color, main="Best B1 separator, random seeds") condPlot(modules2.robust, which.min(tt.B1.2$p.value), ALL.normed, col=color, main="Best B1 separator, smart seeds") ################################################### ### chunk number 68: check B1 same module ################################################### B1.cor <- c(cor(scores1[,which.min(tt.B1.1$p.value)], scores2[,which.min(tt.B1.2$p.value)]), cor(getFeatureMatrix(modules1.robust, mods=which.min(tt.B1.1$p.value)), getFeatureMatrix(modules2.robust, mods=which.min(tt.B1.2$p.value)))) B1.cor ################################################### ### chunk number 69: sessioninfo ################################################### toLatex(sessionInfo())