###################################################
### chunk number 1: rawdata
###################################################
library(edgeR)
WD <- getwd()
dataPath <- paste(.find.package("edgeR"),"data",sep="/")
fn <- dir(dataPath, "txt")
fn
setwd(dataPath)
head(read.table(fn[1],header=TRUE))
dataList <- readDGE(fn,header=TRUE)
setwd(WD)


###################################################
### chunk number 2: createDGEList
###################################################
cls <- gsub("[0-9]","",colnames(dataList$data))
cls
keep <- rowSums(dataList$data) > 8  # to keep the compute time lower
sum(keep)
d<-DGEList(data=dataList$data[keep,],group=cls,lib.size=dataList$lib.size)  # call constructor
d


###################################################
### chunk number 3: moderation
###################################################
alpha <- alpha.approxeb(d)
alpha
ms <- deDGE(d,alpha=alpha$alpha)
ms


###################################################
### chunk number 4: MAplot
###################################################
exactP <- exactTestNB(ms$pseudo, ms$group, pair=c("Tu","NC"), ms$M * ms$ps$p.common,  ms$r, verbose=TRUE)
exactPadj <- p.adjust(exactP,"fdr")
k <- (exactPadj<.05)
plotMA( ms, xlim=c(-16,-5), ylim=c(-6,6), xlab="Relative Abundance", ylab="log Ratio", col=c("black","blue")[k+1])


###################################################
### chunk number 5: topTags
###################################################
topTags(ms)


###################################################
### chunk number 6: poisson
###################################################
set.seed(101)
n <- 1000
lib.sizes<-c(40000,50000,38000,40000)
p<-runif(n,min=.0001,.001)
mu<-outer(p,lib.sizes)
mu[1:5,3:4]<-mu[1:5,3:4]*8
y<-matrix(rpois(4*n,lambda=mu),nrow=n)
dP<-DGEList(data=y,group=rep(1:2,each=2),lib.size=lib.sizes)
msP<-deDGE(dP,doPoisson=TRUE)
msP


###################################################
### chunk number 7: poissonstats
###################################################
topTags(msP)
plotMA(msP,col=c( rep("blue",5), rep("black",n-5) ))


###################################################
### chunk number 8: setup
###################################################
sessionInfo()