######################################################################
# R commands         				   	Genomic Data Science 
# Lab 14: Advance DE                     University of South Carolina
######################################################################
# This file contains the R commands for the lab.
#
# Lines beginning with the symbol '#' are comments in R.  All other
# lines contain code.
#
# In R for Windows, you may wish to open this file from the menu bar
# (File:Display file); you can then copy commands into the command
# window.  (Use the mouse to highlight one or more lines; then
# right-click and select "Paste to console".)
######################################################################


###################################################
### chunk: eBayes
###################################################
library("limma")
library("hgu95av2.db")
library("genefilter")
library("ALL")
data("ALL")

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: rowSds
###################################################
library("genefilter")
sds = rowSds(exprs(ALL_bcrneg))
sh = shorth(sds)
sh


###################################################
### chunk: histsds
###################################################
hist(sds, breaks=50, col="mistyrose", xlab="standard deviation")
abline(v=sh, col="blue", lwd=3, lty=2)




###################################################
### chunk: unspecific
###################################################
ALLsfilt = ALL_bcrneg[sds>=sh, ]



###################################################
### chunk: ttest
###################################################
table(ALLsfilt$mol.biol)
tt = rowttests(ALLsfilt, "mol.biol")
names(tt)



###################################################
### chunk: Design matrix
###################################################

mutation = as.numeric(ALLsfilt$mol.biol=="BCR/ABL") 
design  = cbind(mean = 1, diff = mutation)

###################################################
### chunk: limma
###################################################
fit = lmFit(exprs(ALLsfilt), design)
fit = eBayes(fit)
topTable(fit, coef=2)
Syms = unlist(mget(featureNames(ALLsfilt), 
    env = hgu95av2SYMBOL)) 
fit$gene<-Syms
topTable(fit, coef = "diff", adjust.method = "fdr", 
    sort.by = "p", genelist = fit$gene)
    
top2<-topTable(fit, coef = "diff", adjust.method = "fdr", 
    sort.by = "p", genelist = fit$gene, number=20)
    


###################################################
mtyp = ALLsfilt$mol.biol
sel = rep(1:2, rev(table(mtyp)))
j<-1
top<-which(rownames(ALLsfilt)==rownames(top2)[j])
plot(exprs(ALLsfilt)[top, order(mtyp)], pch=c(1,15)[sel], 
    col=c("black", "red")[sel], 
    main=featureNames(ALLsfilt)[j], 
    ylab=expression(log[2]~expression~level))
legend("bottomleft", col=c("black", "red"), 
    pch=c(1,15), levels(mtyp), bty="n")



###################################################
### chunk: volcano plots
###################################################

pdf("/Users/yenyiho/Documents/BioinfoClass/Notes/FY13/Lecture27/T-volcano.pdf")
plot(tt$dm, -log10(tt$p.value), xlab=expression(paste(log[2], " fold change", sep="")), ylab=expression(paste(-log[10](p))), ylim=c(0, 14))
points(fit$coefficients[,2], -log10(fit$p.value[,2]),  col=2)
legend(-1.5, 14, c("T test", "eBayes"), text.col=c(1,2), pch=c(1,1), col=c(1,2),cex=2, bty="n")
dev.off()