################################################### ### chunk number 1: setup ################################################### options(width = 40) ################################################### ### chunk number 2: getting ALLfilt_bcrneg ################################################### library(ALL) library(hgu95av2.db) data(ALL) bcell <- grep("^B", as.character(ALL$BT)) types <- c("NEG", "BCR/ABL") moltyp <- which(as.character(ALL$mol.biol) %in% types) # subsetting ALL_bcrneg <- ALL[, intersect(bcell, moltyp)] ALL_bcrneg$BT <- factor(ALL_bcrneg$BT) ALL_bcrneg$mol.biol <- factor(ALL_bcrneg$mol.biol) # nonspecific filter: remove genes that does not ## show much variation across samples library(genefilter) filt_bcrneg <- nsFilter(ALL_bcrneg, var.cutoff=0.5) ALLfilt_bcrneg <- filt_bcrneg$eset ################################################### ### chunk number 3: using KEGG ################################################### library(KEGG.db) library(GSEABase) gsc <- GeneSetCollection(ALLfilt_bcrneg, setType=KEGGCollection()) Am <- incidence(gsc) ################################################### ### chunk number 4: nsF ################################################### nsF <- ALLfilt_bcrneg[colnames(Am), ] ################################################### ### chunk number 5: exercise 1 ################################################### dim(nsF) dim(Am) nGene <- rowSums(Am) rownames(Am)[nGene < 10] sort(nGene, decreasing=TRUE)[1] KEGGPATHID2NAME[["05200"]] ################################################### ### chunk number 6: compute per-gene test ################################################### rtt <- rowttests(nsF, "mol.biol") names(rtt) rttStats <- rtt$statistic ################################################### ### chunk number 7: reduce am ################################################### selectedRows <- (rowSums(Am) > 10) Am2 <- Am[selectedRows, ] ################################################### ### chunk number 8: compute z score ################################################### tA <- as.vector(Am2 %*% rttStats) tAadj <- tA /sqrt(rowSums(Am2)) names(tAadj) <- rownames(Am2) ################################################### ### chunk number 9: exercise 2 ################################################### smPW <- tAadj[tAadj < -5] mget(names(smPW),KEGGPATHID2NAME) lgPW <- tAadj[tAadj > 5] mget(names(lgPW), KEGGPATHID2NAME) ################################################### ### chunk number 10: KEGG2heatmap ################################################### KEGG2heatmap("03010", nsF, "hgu95av2") ################################################### ### chunk number 11: permutation ################################################### library(Category) set.seed(123) pvals <- gseattperm(nsF, nsF$mol.biol, Am2, 1000) pvalCut <- 0.05 lowC <- rownames(pvals)[pvals[, 1] <= pvalCut] unlist(getPathNames(lowC), use.names=FALSE) ################################################### ### chunk number 12: remove probes ################################################### ALLfilt <- nsFilter(ALL_bcrneg, require.entez=TRUE, remove.dupEntrez=TRUE, require.CytoBand=TRUE, var.func=IQR, var.cutoff=0.5)$eset ################################################### ### chunk number 13: t-stats for ALLfilt ################################################### library(limma) design <- model.matrix(~0 + ALLfilt$mol.biol) colnames(design) <- c("BCR/ABL", "NEG") contr <- c(1, -1) fit1 <- lmFit(ALLfilt, design) fit2 <- contrasts.fit(fit1, contr) fit3 <- eBayes(fit2) tlimma <- topTable(fit3, number=nrow(fit3), adjust.method="none") ## annotation entrezUniverse <- unlist(mget(tlimma$ID, hgu95av2ENTREZID)) tstats <- tlimma$t names(tstats) <- entrezUniverse ################################################### ### chunk number 14: Chrmaplinearmparams ################################################### library(Category) params <- new("ChrMapLinearMParams", conditional=FALSE, testDirection="up", universeGeneIds=entrezUniverse, geneStats=tstats, annotation="hgu95av2", pvalueCutoff=0.01, minSize=4L) ################################################### ### chunk number 15: linearMTest ################################################### lman <- linearMTest(params) lman summary(lman) ################################################### ### chunk number 16: conditional test ################################################### slotNames(params) paramsCond <- params paramsCond@conditional <- TRUE lmanCond <- linearMTest(paramsCond) summary(lmanCond) ################################################### ### chunk number 17: gene selection by t-test ################################################### ttests <- rowttests(ALLfilt_bcrneg, "mol.biol") smPV <- ttests[ttests$p.value < 0.005, ] selectedEntrezIds <- unlist(mget(rownames(smPV), hgu95av2ENTREZID)) entrezUniverse=unlist(mget(featureNames(ALLfilt_bcrneg), hgu95av2ENTREZID)) #library(limma) #design <- model.matrix(~0 + ALLfilt_bcrneg$mol.biol) #colnames(design) <- c("BCR/ABL", "NEG") #contr <- c(1, -1) #fit1 <- lmFit(ALLfilt_bcrneg, design) #fit2 <- contrasts.fit(fit1, contr) #fit3 <- eBayes(fit2) #tlimma <- topTable(fit3, number=nrow(fit3), # adjust.method="none") #smPV <- tlimma[tlimma$P.Value < 0.01, ] #dim(smPV) #selectedEntrezIds <- unlist(mget(smPV$ID, hgu95av2ENTREZID)) #entrezUniverse=unlist(mget(featureNames(ALLfilt_bcrneg), hgu95av2ENTREZID)) ################################################### ### chunk number 18: GOHyperGParams ################################################### library(GOstats) hgCutoff <- 0.001 GOparams <- new("GOHyperGParams", geneIds=selectedEntrezIds, universeGeneIds=entrezUniverse, annotation="hgu95av2.db", ontology="BP", pvalueCutoff=0.001, conditional=TRUE, testDirection="over") ################################################### ### chunk number 19: hyperGTest eval=FALSE ################################################### ## hgOver <- hyperGTest(GOparams) ## class(hgOver) ## summary(hgOver) ################################################### ### chunk number 20: use htmlReport to generate an html report eval=FALSE ################################################### ## showMethods("htmlReport") ## htmlReport(hgOver, file="hgResult.html")