### R code from vignette source 'vignettes/EDASeq/inst/doc/EDASeq.Rnw'

###################################################
### code chunk number 1: EDASeq.Rnw:37-38
###################################################
options(width=60)


###################################################
### code chunk number 2: data
###################################################
require(EDASeq)
require(yeastRNASeq)
require(leeBamViews)


###################################################
### code chunk number 3: import-raw
###################################################
files <- list.files(file.path(system.file(package = "yeastRNASeq"),
                    "reads"), pattern = "fastq", full.names = TRUE)
names(files) <- gsub("\\.fastq.*", "", basename(files))
met <- DataFrame(conditions=c(rep("mut",2),rep("wt",2)),
                  row.names=names(files))
fastq <- FastqFileList(files)
elementMetadata(fastq) <- met
fastq


###################################################
### code chunk number 4: import-aligned
###################################################
files <- list.files(file.path(system.file(package = "leeBamViews"),
                    "bam"), pattern = "bam$", full.names = TRUE)
names(files) <- gsub("\\.bam", "", basename(files))

gt <- gsub(".*/", "", files)
gt <- gsub("_.*", "", gt)
lane <- gsub(".*(.)$", "\\1", gt)
geno <- gsub(".$", "", gt)

pd <- DataFrame(geno=geno, lane=lane, row.names=paste(geno,lane,sep="."))

bfs <- BamFileList(files)
elementMetadata(bfs) <- pd
bfs


###################################################
### code chunk number 5: plot-total
###################################################
colors <- c(rep(rgb(1,0,0,alpha=0.7),2),rep(rgb(0,0,1,alpha=0.7),2),
            rep(rgb(0,1,0,alpha=0.7),2),rep(rgb(0,1,1,alpha=0.7),2))
barplot(bfs,las=2,col=colors)


###################################################
### code chunk number 6: plot-mean-qual
###################################################
plotQuality(bfs,col=colors,lty=1)
legend("topright",unique(elementMetadata(bfs)[,1]),
       fill=unique(colors))


###################################################
### code chunk number 7: plot-qual-1
###################################################
plotQuality(bfs[[1]],cex.axis=.8)


###################################################
### code chunk number 8: plot-mapped-1
###################################################
barplot(bfs[[1]],las=2)


###################################################
### code chunk number 9: plot-nt-1
###################################################
plotNtFrequency(bfs[[1]])


###################################################
### code chunk number 10: load-gene-level
###################################################
data(geneLevelData)
head(geneLevelData)


###################################################
### code chunk number 11: load-lgc
###################################################
data(yeastGC)
head(yeastGC)
data(yeastLength)
head(yeastLength)


###################################################
### code chunk number 12: filter
###################################################
filter <- apply(geneLevelData,1,function(x) mean(x)>10)
table(filter)
common <- intersect(names(yeastGC),rownames(geneLevelData[filter,])) 
length(common)


###################################################
### code chunk number 13: create-object
###################################################
feature <- data.frame(gc=yeastGC,length=yeastLength)
data <- newSeqExpressionSet(exprs=as.matrix(geneLevelData[common,]),
                            featureData=feature[common,],
                            phenoData=data.frame(
                              conditions=c(rep("mut",2),rep("wt",2)),
                              row.names=colnames(geneLevelData)))
data


###################################################
### code chunk number 14: show-data
###################################################
head(exprs(data))
pData(data)
head(fData(data))


###################################################
### code chunk number 15: show-offset
###################################################
head(offst(data))


###################################################
### code chunk number 16: boxplot-genelevel
###################################################
boxplot(data,col=colors[1:4])


###################################################
### code chunk number 17: md-plot
###################################################
MDPlot(data,c(1,3))


###################################################
### code chunk number 18: plot-mean-var-mut
###################################################
meanVarPlot(data[,1:2],log=T)


###################################################
### code chunk number 19: plot-mean-var-all
###################################################
meanVarPlot(data,log=T)


###################################################
### code chunk number 20: plot-gc
###################################################
biasPlot(data,"gc",log=T,ylim=c(1,5))


###################################################
### code chunk number 21: boxplot-gc
###################################################
lfc <- log(exprs(data)[,3]+0.5)-log(exprs(data)[,1]+0.5)
biasBoxplot(lfc,fData(data)$gc)


###################################################
### code chunk number 22: normalization
###################################################
dataWithin <- withinLaneNormalization(data,"gc",which="full")
dataNorm <- betweenLaneNormalization(dataWithin,which="full")


###################################################
### code chunk number 23: plot-gc-norm
###################################################
biasPlot(dataNorm,"gc",log=T,ylim=c(1,5))


###################################################
### code chunk number 24: boxplot-norm
###################################################
boxplot(dataNorm,col=colors)


###################################################
### code chunk number 25: norm-offset
###################################################
dataOffset <- withinLaneNormalization(data,"gc",
                                      which="full",offset=TRUE)
dataOffset <- betweenLaneNormalization(dataOffset,
                                       which="full",offset=TRUE)


###################################################
### code chunk number 26: edger
###################################################
library(edgeR)
design <- model.matrix(~conditions,data=pData(data))
disp <- estimateGLMCommonDisp(exprs(dataOffset), 
                              design, offst(dataOffset))
fit <- glmFit(exprs(dataOffset),design,disp,offst(dataOffset))
lrt <- glmLRT(exprs(dataOffset),fit,coef=2)
head(lrt$table)


###################################################
### code chunk number 27: deseq (eval = FALSE)
###################################################
## library(DESeq)
## counts <- as(dataNorm,"CountDataSet")
## sizeFactors(counts) <- rep(1,4)
## counts <- estimateDispersions(counts)
## res <- nbinomTest(counts, "wt", "mut" )
## head(res)


###################################################
### code chunk number 28: sessionInfo
###################################################
toLatex(sessionInfo())