######################################################################
# R commands         				   	Genomic Data Science
# Lab 13                                 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".)
######################################################################

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)
emat<-exprs(ALL_bcrneg)


featureNames(ALL_bcrneg)[1:10]
featureNames(ALL_bcrneg)

library("genefilter")
sds = rowSds(exprs(ALL_bcrneg))
sh = shorth(sds)
sh
ALLsfilt = ALL_bcrneg[sds>=sh, ]
dim(exprs(ALLsfilt))
ALLsfilt$mol.biol<-factor(ALLsfilt$mol.biol)
table(ALLsfilt$mol.biol)

######## Exercise: Perform T-test for every probe set 

#################################
# Version 1: t.test 
#################################
emat<-exprs(ALLsfilt)
dim(emat)
tout<-matrix(NA, nrow=nrow(emat), ncol=3)
rownames(tout)<-rownames(emat)
colnames(tout)<-c("statisitcs", "dm", "p.value")
for(i in 1:nrow(emat)){
	tmp<-t.test(emat[i,] ~ ALLsfilt$mol.biol, var.equal=T)
	tout[i,1]<-tmp$statistic
	tout[i,3]<-tmp$p.value
	tout[i,2]<-tmp$estimate[1] - tmp$estimat[2]
}
tout[1:10,]
#################################
# Version 2: rowttest (much faster)
#################################
tt = rowttests(ALLsfilt, "mol.biol")
names(tt)

###
pdf("FDR.pdf")
hist(tt$p.value, breaks=50, col="mistyrose", xlab="p-value", 
     main="FDR")
abline(h=135, lty=2)
dev.off()

ALLsrest = ALL_bcrneg[sds<sh, ]
ttrest = rowttests(ALLsrest, "mol.biol")
hist(ttrest$p.value, breaks=50, col="lightblue", xlab="p-value", 
     main="Removed")

###################################################
### chunk: adjustp
###################################################
#BiocManager::install("qvalue")
library("qvalue")
mt = p.adjust(tt$p.value, method="BH")
fdr<-qvalue(tt$p.value)$qvalues
Table<-cbind(tt, fdr)
Table[1:10,]
o<-order(Table$fdr)
Table<-Table[o,]
Table[1:10,]


###################################################
### chunk: Make a Volcano plot 
###################################################
plot(tt$dm, -log10(tt$p.value), pch=1, 
    xlab = expression(mean~log[2]~fold~change),
    ylab = expression(-log[10](p)))
abline(v=1, lty=2, col="grey")
abline(v=-1, lty=2, col="grey")
abline(h=5, lty=2, col="grey")



###################################################
### chunk: Annotation 
###################################################
g<-rownames(Table)
###biocLite("hgu95av2.db")
library("hgu95av2.db")
syms <- unlist(mget(g, hgu95av2SYMBOL))
Table<-cbind(g, syms, Table)
colnames(Table)<-c("probe", "gene symbol", "t", "log2FC", "p.value", "fdr")
Table[1:10,]