################################################### ### chunk number 1: setup ################################################### options(width = 40) ################################################### ### chunk number 2: get dataset ################################################### 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 library(genefilter) filt_bcrneg <- nsFilter(ALL_bcrneg, require.entrez=TRUE, require.GOBP=TRUE, remove.dupEntrez=TRUE, feature.exclude="^AFFX", var.cutoff=0.5) ALLfilt_bcrneg <- filt_bcrneg$eset ################################################### ### chunk number 3: rowttests ################################################### tt <- rowttests(ALLfilt_bcrneg, "mol.biol") plot(tt$dm, -log10(tt$p.value), pch=".", xlab=expression(mean~log[2]~fold~change), ylab=expression(-log[2](p))) ################################################### ### chunk number 4: linear model 2 ################################################### model.matrix(~mol.biol + 0, ALLfilt_bcrneg) ################################################### ### chunk number 5: linear model 1 ################################################### model.matrix(~ mol.biol, ALLfilt_bcrneg) ################################################### ### chunk number 6: design matrix ################################################### library(limma) #cl = as.numeric(ALLfilt_bcrneg$mol.biol=="BCR/ABL") #design <- cbind(mean=1, diff=cl) design <- model.matrix( ~mol.biol + 0, ALLfilt_bcrneg) colnames(design) <- c("BCR_ABL", "NEG") # contr <- makeContrasts(BCR_ABL-NEG, levels=design) contr <- c(1, -1) ################################################### ### chunk number 7: linear models and ebayes ################################################### fit <- lmFit(exprs(ALLfilt_bcrneg), design) fit1 <- contrasts.fit(fit, contr) fit2 <- eBayes(fit1) #syms <- unlist(mget(featureNames(ALLfilt_bcrneg), hgu95av2SYMBOL)) topTable(fit2, adjust.method="BH", number=5) ################################################### ### chunk number 8: html eval=FALSE ################################################### ## library(annotate) ## top20Gene <- topTable(fit2, adjust.method="BH", ## number=20, genelist=syms) ## htmlpage(genelist=as.data.frame(top20Gene$ID), ## othernames=top20Gene, ## filename="top20gene.html", ## table.head=c("probe ID", names(top20Gene))) ## broweURL("top20Gene.html") ################################################### ### chunk number 9: annotation ################################################### library(org.Hs.eg.db) org.Hs.eg() org.Hs.eg_dbInfo() org.Hs.egGENENAME ################################################### ### chunk number 10: Lkeys ################################################### map <- org.Hs.egGENENAME map head(Lkeys(map)) ## probeset id map[["1000"]] revmap(map)[["adenosine deaminase"]] ## reversible ################################################### ### chunk number 11: working with GO ################################################### library(GO.db) ls("package:GO.db") ## find children as.list(GOMFCHILDREN["GO:0008094"]) ## all the descendants (children, grandchildren, and so on) as.list(GOMFOFFSPRING["GO:0008094"])