##load up required libraries

library(cluster) library(Biobase) library(annotate) library(genefilter) library(golubEsets)

############################ ## data processing ###########################

data(golubTrain) golubTrain data(golubTest) golubTest

LS<-exprs(golubTrain) cl<-golubTrain$ALL.AML TS<-exprs(golubTest) clts<-golubTest$ALL.AML

##mmfilt is in the library

LS[LS<100]<-100 TS[TS<100]<-100 LS[LS>16000]<-16000 TS[TS>16000]<-16000

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

LSsub <- log10(LS[good,]) TSsub <- log10(TS[good,]) LTsub <- cbind(LSsub, TSsub)

############################## ## Annotate : optional ##############################

##probably want to do a ttest filter here ##but we'll leave that to you

ll<- read.annotation("hgu68ll") #locus link gname <- row.names(LSsub)[50] llname<-get(gname, env=ll) locuslink(llname) ##open your browser to the right thing

chroms <- read.annotation("hgu68Chrom")

whichC <- multiget(row.names(LSsub[1:100,]), env=chroms, iffail=NA) cc <-as.numeric(unlist(whichC)) hist(cc)

################################### ## Clustering #################################

gf <- gapFilter(550, 800, .1) ff <- filterfun(gf) xx <- 10^LSsub ##to get back to original units good <- genefilter(xx, gf) sum(good) # should be 1093

LS2 <- LSsub[good,]

library(mva) clust.res<-hclust(as.dist(1-cor(LS2)), method = "complete") plot.hclust(clust.res, labels=cl,hang = 0.1)

clust.res<-hclust(dist(t(LS2)), method = "average") plot.hclust(clust.res, labels=cl,hang = 0.1)

pamLS <- pam(as.dist(1-cor(LS2)), 3, diss=TRUE) plot(pamLS)

####################################### ## classification ##################################

library(class)

ALL <- c(as.character(golubTrain$ALL), as.character(golubTest$ALL)) BT <- c(as.character(golubTrain$T.B), as.character(golubTest$T.B)) BTs <- ifelse(BT=="B-cell", "B", "T") BTs <- ifelse(BT=="T-cell", "T", BTs) BTs <- ifelse(BT=="NA", "", BTs) newALL <- factor(paste(ALL, BTs, sep="")) nALL.train <- newALL[1:38] nALL.test <-newALL[39:72]

##now do an Anova filter af <- Anova(nALL.train, p=0.0000001) ffA <- filterfun(af) wh.ALL <- genefilter(LSsub, af) #use only training data sum(wh.ALL) #how many

##how does it work on the test set pred.test <- knn(t(LSsub[wh.ALL,]), t(TSsub[wh.ALL,]), nALL.train, k=5) table(pred.test, nALL.test)

##how does it work on the training set pred.tr <- knn(t(LSsub[wh.ALL,]), t(LSsub[wh.ALL,]), nALL.train, k=5) table(pred.tr, nALL.train)

## ##use knn and cross-validation to choose the best k!

knn.cvk<-function(LS, cl, k=1:5, l=0, prob=FALSE, use.all=TRUE) { classes<-unique(cl) cv.err<-rep(0,length(k)) cl.pred<-matrix(NA, nrow(LS), length(k)) for(j in (1:length(k))) { cl.pred[,j]<-knn.cv(LS, cl, k[j], l, prob, use.all) cv.err[j]<-sum(cl!=cl.pred[,j]) } k0<-k[which.min(cv.err)] return(k=k0,pred=classes[cl.pred[,which.min(cv.err)]]) }

whichone <- knn.cvk(t(LSsub[wh.ALL,]), nALL.train)

##Tree based classifiers

library(rpart)

dorpart <- function(LS, TS, cl) { LS <- as.data.frame(LS) TS <- as.data.frame(TS) dimnames(TS)[[2]] <- dimnames(LS)[[2]] <- paste("V", 1:ncol(LS), sep = "") fit <- rpart(factor(cl) ~ ., data = LS) predict.rpart(fit, TS, type = "class") }

all.rp <- dorpart(t(LSsub[wh.ALL,]), t(TSsub[wh.ALL,]), nALL.train) table(all.rp, nALL.test)

##Discriminant Analysis

library(MASS)

dolda <- function(LS, TS, cl, prior, pmethod="plug-in") { ndata <- data.frame(LS) ndata$cl <- cl if( missing(prior) ) mn <- lda(cl~., data=ndata) else mn <- lda(cl~., data=ndata, prior=prior) pdata<-data.frame(TS) preds <- predict(mn, pdata, method=pmethod) rval <- list(lda = mn, preds=preds) class(rval) <- "mylda" rval }

lda.test <- dolda(t(LSsub[wh.ALL,]), t(TSsub[wh.ALL,]), nALL.train) table(lda.test$preds$class, nALL.test)

##how does it work on the training set ##not very meaningful lda.pred.tr <- dolda(t(LSsub[wh.ALL,]), t(LSsub[wh.ALL,]), nALL.train) table(lda.pred.tr$preds$class, nALL.train)