###################################################
### chunk number 1: set width
###################################################
options(width=60)
options(continue=" ")
try(X11.options(type="xlib"), silent=TRUE)


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