###################################################
### chunk number 1:
###################################################
set.seed(1234)
x <- rnorm(5)
print(x)

###################################################
### chunk number 2:
###################################################
print(3:10)

###################################################
### chunk number 3:
###################################################
print(x[2:4])

###################################################
### chunk number 4:
###################################################
inds <- c(2,3,4)
print(x[inds])

###################################################
### chunk number 5:
###################################################
inds <- c(low=2,mid=3,hi=4)
print(inds)
print(x[ inds[c("low","hi")] ])

###################################################
### chunk number 6:
###################################################
myList <- list(a=x, b=inds)
print(myList[[1]])
print(myList[["a"]])
print(myList$b)


###################################################
### chunk number 7:
###################################################
library(cluster)
library(Biobase)
library(annotate)
library(genefilter)
library(AggPred)
library(sma)


###################################################
### chunk number 8:
###################################################
library(golubEsets)
data(golubTrain)
data(golubTest)
print(golubTrain)
print(table(golubTrain$ALL.AML))

###################################################
### chunk number 9:
###################################################
LS <- exprs(golubTrain)
cl <- golubTrain$ALL.AML
TS <- exprs(golubTest)
clts <- golubTest$ALL.AML

###################################################
### chunk number 10:
###################################################
LS[LS<100]<-100
TS[TS<100]<-100
LS[LS>16000]<-16000
TS[TS>16000]<-16000


###################################################
### chunk number 11:
###################################################
mmfilt <- function(r=5, d=500, na.rm=TRUE) {
function(x) {
minval <- min(x, na.rm=na.rm)
maxval <- max(x, na.rm=na.rm)
(maxval/minval > 5) && (maxval-minval > 500)
}
}

mmfun <- mmfilt()

ffun <- filterfun(mmfun)

good <- genefilter(cbind(LS, TS), ffun )

sum(good) ##I get 3571 not 3051!

LSsub <- log10(LS[good,])
TSsub <- log10(TS[good,])
print(dim(LSsub))


###################################################
### chunk number 12:
###################################################
cor.all<-cor(LSsub)
plot.cor(cor.all[order(cl), order(cl)], labels=sort(cl),
labcols=unclass(sort(cl)),
title="Correlation matrix, all genes",zlim=c(-1,1))

###################################################
### chunk number 13:
###################################################
library(ellipse)
plotcorr(cor.all[order(cl), order(cl)], labels=sort(cl),
labcols=unclass(sort(cl)),
title="Correlation matrix, all genes")

###################################################
### chunk number 14:
###################################################
gf <- gapFilter(250, 500, .1)
ff <- filterfun(gf)
good <- genefilter(10^LSsub, gf)
sum(good) # should be 1091

cor.gf <- cor(LS[good,])


###################################################
### chunk number 15:
###################################################
plot.cor(cor.gf[order(cl), order(cl)], labels=sort(cl),
labcols=unclass(sort(cl)),
title="Correlation matrix, filtered genes",zlim=c(-1,1))

###################################################
### chunk number 16:
###################################################
plotcorr(cor.gf[order(cl), order(cl)], labels=sort(cl),
labcols=unclass(sort(cl)),
title="Correlation matrix, filtered genes")

###################################################
### chunk number 17:
###################################################
library(edd)
ALLtrDists <- edd.unsupervised( golubTrain[, golubTrain$ALL=="ALL"] )
print(table(ALLtrDists))

###################################################
### chunk number 18:
###################################################
golURL <- url("http://www-genome.wi.mit.edu/mpr/data_set_ALL_AML_train.txt","r")

###################################################
### chunk number 19:
###################################################
golVec <- scan(golURL, sep="\t", what="")


###################################################
### chunk number 20:
###################################################
print(length(golVec))
print(length(golVec)/79)

###################################################
### chunk number 21:
###################################################
golVec <- c(golVec[1:78],"",golVec[-(1:78)])
print(length(golVec)/79)
golMat <- t(matrix(golVec,nr=79))
print(dim(golMat))
print(golMat[1:5,1:5])

###################################################
### chunk number 22:
###################################################
cExprs <- golMat[-1,seq(3,78,2)]
numExprs <- t(apply(cExprs,1,as.numeric))
print(numExprs[1:4,1:4])

###################################################
### chunk number 23:
###################################################
accnos <- golMat[-1,2]
dimnames(numExprs) <- list(accnos,NULL)

###################################################
### chunk number 24:
###################################################
exprSamps <- golMat[1,seq(3,78,2)]

###################################################
### chunk number 25:
###################################################
pheu <- url("http://www-genome.wi.mit.edu/mpr/table_ALL_AML_samples.txt","r")
PH <- readLines(pheu)
print(length(PH))
print(PH[1:5])

###################################################
### chunk number 26:
###################################################
print(PH[10])
print(PH[47])

###################################################
### chunk number 27:
###################################################
print(strsplit(PH[10],"\t"))

###################################################
### chunk number 28:
###################################################
samp <- rep(NA,38)
dx <- rep(NA,38)
cellType <- rep(NA,38)
for (i in 10:47) {
tmp <- strsplit(PH[i],"\t")[[1]]
samp[i-9] <- tmp[1]
dx[i-9] <- tmp[3]
cellType[i-9] <- tmp[6]
}
dx <- substring(dx,1,3)
samp <- c(substring(samp[1:9],1,1),substring(samp[10:38],1,2))
phenoDF <- data.frame(samp=as.numeric(samp), dx=dx,
cellType=cellType)
print(phenoDF[1:4,])

###################################################
### chunk number 29:
###################################################
library(Biobase)
myPh <- new("phenoData", pData=phenoDF[as.numeric(exprSamps),],
varLabels=list(samp="sample code",
dx="diagnosis", cellType="cell type"))
myES <- new("exprSet", exprs=numExprs,
phenoData=myPh, description=names(phenoDF),
annotation="affy", notes="from web")

###################################################
### chunk number 30:
###################################################
print(myES)
print(myES[1:10,])
print(myES$dx[c(1:5,37:38)])

###################################################
### chunk number 31:
###################################################
print(geneNames(myES)[c(1,1000,2000)])
print(exprs(myES)[c(1,1000,2000),1:3])