Skip to content.

bioconductor.org

Bioconductor is an open source and open development software project
for the analysis and comprehension of genomic data.

Sections

GSEA.R

################################################### ### 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()

News
2009-10-26

BioC 2.5, consisting of 352 packages and designed to work with R 2.10.z, was released today.

2009-01-07

R, the open source platform used by Bioconductor, featured in a series of articles in the New York Times.