###################################################
### 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"])