################################################### ### chunk number 1: load0 ################################################### library("BiocCaseStudies") library("Biobase") library("Category") options(width=56) ################################################### ### chunk number 2: ALL ################################################### library("ALL") data("ALL") ################################################### ### chunk number 3: bcrabl ################################################### bcell = grep("^B", as.character(ALL$BT)) moltyp = which(as.character(ALL$mol.biol) %in% c("NEG", "BCR/ABL")) ALL_bcrneg = ALL[, intersect(bcell, moltyp)] ALL_bcrneg$mol.biol = factor(ALL_bcrneg$mol.biol) ################################################### ### chunk number 4: nsFilter ################################################### library("genefilter") ALLfilt_bcrneg = nsFilter(ALL_bcrneg, var.cutoff=0.5)$eset ################################################### ### chunk number 5: ans1 ################################################### table(ALLfilt_bcrneg$mol.biol) ################################################### ### chunk number 6: KEGGimat ################################################### library("GSEABase") gsc = GeneSetCollection(ALLfilt_bcrneg, setType=KEGGCollection()) Am = incidence(gsc) dim(Am) ################################################### ### chunk number 7: EsetfromKEGG ################################################### nsF = ALLfilt_bcrneg[colnames(Am),] ################################################### ### chunk number 8: compttests ################################################### rtt = rowttests(nsF, "mol.biol") rttStat = rtt$statistic ################################################### ### chunk number 9: reducetoInt ################################################### selectedRows = (rowSums(Am)>10) Am2 = Am[selectedRows, ] ################################################### ### chunk number 10: pctests ################################################### tA = as.vector(Am2 %*% rttStat) tAadj = tA/sqrt(rowSums(Am2)) names(tA) = names(tAadj) = rownames(Am2) ################################################### ### chunk number 11: Aufpassen ################################################### # stopifnot(sum(tAadj<(-5))==1L) ################################################### ### chunk number 12: qqplot ################################################### qqnorm(tAadj) ################################################### ### chunk number 13: findSmPW ################################################### library("KEGG.db") smPWs = tAadj[tAadj < (-5)] pwNames = KEGGPATHID2NAME[names(smPWs)] cbind(smPWs, toTable(pwNames)) smPW = tAadj[tAadj < (-8)] pwName = KEGGPATHID2NAME[[names(smPW)]] ################################################### ### chunk number 14: mnplot ################################################### KEGGmnplot(names(smPW), nsF, "hgu95av2", nsF$"mol.biol", pch=16, col="darkblue") ################################################### ### chunk number 15: heatmap ################################################### sel = as.integer(nsF$mol.biol) KEGG2heatmap(names(smPW), nsF, "hgu95av2", col=colorRampPalette(c("white", "darkblue"))(256), ColSideColors=c("black", "white")[sel]) ################################################### ### chunk number 16: nlevelsnsFmolbiol ################################################### stopifnot(nlevels(nsF$mol.biol)==2L) ################################################### ### chunk number 17: ttperm ################################################### library(Category) set.seed(123) NPERM = 1000 pvals = gseattperm(nsF, nsF$mol.biol, Am2, NPERM) pvalCut = 0.025 lowC = names(which(pvals[, 1]<=pvalCut)) highC = names(which(pvals[, 2]<=pvalCut)) ################################################### ### chunk number 18: fmt ################################################### fmt = function(x) as.vector(mapply(paste, names(x), x, sep=": ")) ################################################### ### chunk number 19: lowC1 eval=FALSE ################################################### ## getPathNames(lowC) ################################################### ### chunk number 20: lowC2 ################################################### fmt(getPathNames(lowC)) ################################################### ### chunk number 21: highC1 eval=FALSE ################################################### ## getPathNames(highC) ################################################### ### chunk number 22: highC2 ################################################### fmt(getPathNames(highC))[1:5] ################################################### ### chunk number 23: applypvals ################################################### apply(pvals, 2, min) rownames(pvals)[apply(pvals, 2, which.min)] ################################################### ### chunk number 24: pvalcomp ################################################### permpvs = pmin(pvals[,1], pvals[,2]) pvsparam = pnorm(tAadj) pvspara = pmin(pvsparam, 1-pvsparam) ################################################### ### chunk number 25: pvplot ################################################### plot(permpvs, pvspara, xlab="Permutation p-values", ylab="Parametric p-values") ################################################### ### chunk number 26: filterMissingMAP ################################################### ## depending on which annotation infrastructure we use ## hgu95av2MAP will either be an environment or an ## AnnDbBimap object fnames = featureNames(ALLfilt_bcrneg) if(is(hgu95av2MAP, "environment")){ chrLocs = mget(fnames, hgu95av2MAP) mapping = names(chrLocs[sapply(chrLocs, function(x) !all(is.na(x)))]) }else{ mapping = toTable(hgu95av2MAP[fnames])$probe_id } psWithMAP = unique(mapping) nsF2 = ALLfilt_bcrneg[psWithMAP, ] ################################################### ### chunk number 27: findAmap ################################################### EGtable = toTable(hgu95av2ENTREZID[featureNames(nsF2)]) entrezUniv = unique(EGtable$gene_id) chrMat = MAPAmat("hgu95av2", univ=entrezUniv) rSchr = rowSums(chrMat) ################################################### ### chunk number 28: MAPAmat ################################################### dim(chrMat) ################################################### ### chunk number 29: reduceCols ################################################### chrMat = chrMat[rowSums(chrMat) >= 5, ] dim(chrMat) ################################################### ### chunk number 30: reorderChrMat ################################################### EGlist = mget(featureNames(nsF2), hgu95av2ENTREZID) EGIDs = sapply(EGlist, "[", 1) idx = match(EGIDs, colnames(chrMat)) chrMat = chrMat[, idx] ################################################### ### chunk number 31: extraParanoia ################################################### stopifnot(!any(is.na(idx)), all(listLen(EGlist)==1L)) ################################################### ### chunk number 32: overlap ################################################### Ams = Am2[union(lowC, highC),] Amx = Ams %*% t(Ams) minS = outer(diag(Amx), diag(Amx), pmin) overlapIndex = Amx/minS ################################################### ### chunk number 33: overlapSetsPlot ################################################### library("lattice") myFun = function(x) { dd.row = as.dendrogram(hclust(dist(x))) row.ord = order.dendrogram(dd.row) dd.col = as.dendrogram(hclust(dist(t(x)))) col.ord = order.dendrogram(dd.col) colnames(x) = sapply(getPathNames(colnames(x)), "[", 1L) levelplot(x[row.ord,col.ord], scales = list(x = list(rot = 90)), xlab="", ylab="", main="overlapIndex", col.regions=colorRampPalette(c("white", "darkblue"))(256)) } print(myFun(overlapIndex)) ################################################### ### chunk number 34: ijmin ################################################### ijord = function(m, i) { m[lower.tri(m, diag=TRUE)] = NA o=order(m, decreasing=TRUE, na.last=TRUE)[i] cbind((o-1L)%/%nrow(m)+1L, (o-1L)%%nrow(m)+1L) } nm = getPathNames(colnames(overlapIndex)) pathnames = as.vector(mapply(paste, names(nm), sapply(nm, "[", 1L), sep=": ")) k=ijord(overlapIndex, 1:2) ################################################### ### chunk number 35: geneoverlap ################################################### rowSums(Ams)[c("04510", "04512", "04514", "04940")] Amx["04512", "04510"] Amx["04940", "04514"] ################################################### ### chunk number 36: lmfits ################################################### P04512 = Ams["04512",] P04510 = Ams["04510",] lm1 = lm(rttStat ~ P04512) summary(lm1)$coefficients lm2 = lm(rttStat ~ P04510) summary(lm2)$coefficients lm3 = lm(rttStat ~ P04510+P04512) summary(lm3)$coefficients ################################################### ### chunk number 37: stopifnot ################################################### stopifnot( coefficients(summary(lm1))["P04512", "Pr(>|t|)"]<0.05, coefficients(summary(lm2))["P04510", "Pr(>|t|)"]<0.05, coefficients(summary(lm3))["P04512", "Pr(>|t|)"]>0.05, coefficients(summary(lm3))["P04510", "Pr(>|t|)"]<0.05) ################################################### ### chunk number 38: threegroups ################################################### P04512.Only = ifelse(P04512 != 0 & P04510 == 0, 1, 0) P04510.Only = ifelse(P04512 == 0 & P04510 != 0, 1, 0) Both = ifelse(P04512 != 0 & P04510 != 0, 1, 0) lm4 = lm(rttStat ~ P04510.Only + P04512.Only + Both) ################################################### ### chunk number 39: tgshow eval=FALSE ################################################### ## summary(lm4) ################################################### ### chunk number 40: tgdo ################################################### fixedWidthCat(summary(lm4)) ################################################### ### chunk number 41: pairpw2 ################################################### P04514 = Ams["04514",] P04940 = Ams["04940",] P04514.Only = ifelse(P04514 != 0 & P04940 == 0, 1, 0) P04940.Only = ifelse(P04514 == 0 & P04940 != 0, 1, 0) Both = ifelse(P04514 != 0 & P04940 != 0, 1, 0) lm5 = lm(rttStat ~ P04514.Only + P04940.Only + Both) summary(lm5) ################################################### ### chunk number 42: sessionInfo ################################################### mySessionInfo()