################################################### ### chunk number 1: Setup ################################################### library("Biobase") library("annotate") library("Category") library("hgu95av2") library("genefilter") ################################################### ### chunk number 2: Example-subset ################################################### library("ALL") data(ALL) Bcell <- grep("^B", as.character(ALL$BT)) bcrAblOrNegIdx <- which(as.character(ALL$mol.biol) %in% c("NEG", "BCR/ABL")) eset <- ALL[, intersect(Bcell, bcrAblOrNegIdx)] eset$mol.biol <- factor(eset$mol.biol) numBN <- length(eset$mol.biol) ################################################### ### chunk number 3: noEntrez ################################################### entrezIds <- mget(geneNames(eset), envir=hgu95av2LOCUSID) haveEntrezId <- names(entrezIds)[sapply(entrezIds, function(x) !is.na(x))] numNoEntrezId <- length(geneNames(eset)) - length(haveEntrezId) eset <- eset[haveEntrezId, ] ################################################### ### chunk number 4: simplefiltering ################################################### ## Non-specific filtering based on IQR lowQ = rowQ(eset, floor(0.25 * numBN)) upQ = rowQ(eset, ceiling(0.75 * numBN)) iqrs = upQ - lowQ selected <- iqrs > 0.5 nsFiltered <- eset[selected, ] ################################################### ### chunk number 5: reduceto1to1 ################################################### ## Reduce to unique probe <--> gene mapping by keeping largest IQR ## This gives us "unique genes" in the non-specific filtered gene ## set which simplifies further calculations. nsFilteredIqr <- iqrs[selected] uniqGenes <- findLargest(geneNames(nsFiltered), nsFilteredIqr, "hgu95av2") nsFiltered <- nsFiltered[uniqGenes, ] ## basic stats on our non-specific filter result numSelected <- length(geneNames(nsFiltered)) numBcrAbl <- sum(nsFiltered$mol.biol == "BCR/ABL") numNeg <- sum(nsFiltered$mol.biol == "NEG") ################################################### ### chunk number 6: noKEGG ################################################### ## Remove genes with no PATH mapping havePATH <- sapply(mget(geneNames(nsFiltered), hgu95av2PATH), function(x) if (length(x) == 1 && is.na(x)) FALSE else TRUE) numNoPATH<- sum(!havePATH) nsF <- nsFiltered[havePATH, ] ################################################### ### chunk number 7: compAmat ################################################### Am = PWAmat("hgu95av2") egN = unlist(mget(geneNames(nsF), hgu95av2LOCUSID)) sub1 = match(egN, row.names(Am)) Am = Am[sub1,] dim(Am) table(colSums(Am)) ################################################### ### chunk number 8: compttests ################################################### rtt = rowttests(nsF, "mol.biol") rttStat = rtt$statistic ################################################### ### chunk number 9: reducetoInt ################################################### Amat = t(Am) rs = rowSums(Amat) Amat2 = Amat[rs>10,] rs2 = rs[rs>10] nCats = length(rs2) ################################################### ### chunk number 10: pctests ################################################### tA = as.vector(Amat2 %*% rttStat) tAadj = tA/sqrt(rs2) names(tA) = names(tAadj) = row.names(Amat2) ################################################### ### chunk number 11: qqplot ################################################### qqnorm(tAadj) ################################################### ### chunk number 12: findSmPW ################################################### smPW = tAadj[tAadj < -5] pwName = KEGGPATHID2NAME[[names(smPW)]] pwName ################################################### ### chunk number 13: mnplot ################################################### KEGGmnplot(names(smPW), nsF, "hgu95av2", nsF$"mol.biol") ################################################### ### chunk number 14: heatmap ################################################### KEGG2heatmap(names(smPW), nsF, "hgu95av2") ################################################### ### chunk number 15: ttperm ################################################### NPERM = 100 ttp = ttperm(exprs(nsF), nsF$mol.biol, NPERM) permDm = do.call("cbind", lapply(ttp$perms, function(x) x$statistic)) permD = Amat2 %*% permDm pvals = matrix(NA, nr=nCats, ncol=2) dimnames(pvals) = list(row.names(Amat2), c("Lower", "Upper")) for(i in 1:nCats) { pvals[i,1] = sum(permD[i,] < tA[i])/NPERM pvals[i,2] = sum(permD[i,] > tA[i])/NPERM } ord1 = order(pvals[,1]) lowC = (row.names(pvals)[ord1])[pvals[ord1,1]< 0.05] highC = row.names(pvals)[pvals[,2] < 0.05] getPathNames(lowC) getPathNames(highC) lnhC = length(highC) ################################################### ### chunk number 16: findAmap ################################################### AmChr = MAPAmat("hgu95av2", minCount=5) ################################################### ### chunk number 17: sub2ourData ################################################### subC = row.names(AmChr) %in% egN AmChr = AmChr[subC,] dim(AmChr) table(colSums(AmChr)) ################################################### ### chunk number 18: ################################################### sessionInfo()