################################################### ### chunk number 1: options ################################################### options(width=50, digits=2) ################################################### ### chunk number 2: ALL ################################################### library("ALL") data("ALL") ################################################### ### chunk number 3: bcrabl ################################################### bcell = grep("^B", as.character(ALL$BT)) # B-cell ALL moltyp = which(as.character(ALL$mol.biol) %in% c("NEG", "BCR/ABL")) bcrneg = ALL[, intersect(bcell, moltyp)] bcrneg$mol.biol = factor(bcrneg$mol.biol) # drop unused levels ################################################### ### chunk number 4: bcrabl-result ################################################### dim(bcrneg) ################################################### ### chunk number 5: nsFilter ################################################### library("genefilter") bcrneg_filt1 = nsFilter(bcrneg, var.cutoff=0.5)$eset dim(bcrneg_filt1) ################################################### ### chunk number 6: gsc-filter ################################################### library(GSEABase) gsc <- GeneSetCollection(bcrneg_filt1, setType=KEGGCollection()) ################################################### ### chunk number 7: gsc-egs ################################################### gsc gsc[[2]] ################################################### ### chunk number 8: gsc-filter-2 ################################################### ok <- sapply(geneIds(gsc), length)>10 gsc <- gsc[ok] length(gsc) uids <- unique(unlist(geneIds(gsc))) bcrneg_filt2 <- bcrneg_filt1[uids,] dim(bcrneg_filt2) ################################################### ### chunk number 9: rowtttest ################################################### rtt <- rowttests(bcrneg_filt2, "mol.biol") rttStat <- rtt$statistic names(rttStat) <- featureNames(bcrneg_filt2) head(rttStat) ################################################### ### chunk number 10: zstats ################################################### zCalc <- function(ids, tStat) { sum(tStat[ids]) / sqrt(length(ids)) } z <- sapply(geneIds(gsc), zCalc, tStat=rttStat) names(z) <- names(gsc) head(z) ################################################### ### chunk number 11: qqnorm ################################################### qqnorm(z) qqline(z) ################################################### ### chunk number 12: outlier ################################################### z[which.min(z)] ################################################### ### chunk number 13: kegg-id ################################################### keggId <- names(z[which.min(z)]) keggGS <- gsc[[keggId]] keggES <- bcrneg_filt2[keggGS,] library(Category) getPathNames(keggId) ################################################### ### chunk number 14: keggInfo ################################################### KEGGmnplot(keggId, bcrneg_filt2, annotation(bcrneg_filt2), bcrneg_filt2$mol.biol, pch=16, col="darkblue") ################################################### ### chunk number 15: olap ################################################### overlap <- gsc[["04512"]] & gsc[["04510"]] length(GSEABase::geneIds(overlap))