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

###################################################
### code chunk number 1: options
###################################################
options(digits=3)


###################################################
### code chunk number 2: start
###################################################
library("DEXSeq")


###################################################
### code chunk number 3: datapasilla
###################################################
data("pasillaExons", package="pasilla")


###################################################
### code chunk number 4: installpasilla (eval = FALSE)
###################################################
## source("http://www.bioconductor.org/biocLite.R")
## biocLite("pasilla")


###################################################
### code chunk number 5: pData
###################################################
design(pasillaExons)


###################################################
### code chunk number 6: fData
###################################################
head(fData(pasillaExons)[,c(1,2,9:12)])


###################################################
### code chunk number 7: tabtab1
###################################################
tg = table(geneIDs(pasillaExons)) 
tt = table(tg)
stopifnot(tt["36"]==1, tt["16"]==3)


###################################################
### code chunk number 8: tabtab2
###################################################
table(table(geneIDs(pasillaExons))) 


###################################################
### code chunk number 9: sizeFactors1
###################################################
pasillaExons <- estimateSizeFactors(pasillaExons)


###################################################
### code chunk number 10: sizeFactors2
###################################################
sizeFactors(pasillaExons)


###################################################
### code chunk number 11: mffg1
###################################################
head( modelFrameForGene( pasillaExons, "FBgn0010909" ) )
estimateExonDispersionsForModelFrame( modelFrameForGene( pasillaExons, "FBgn0010909" ) )


###################################################
### code chunk number 12: estDisp1
###################################################
pasillaExons <- estimateDispersions(pasillaExons)


###################################################
### code chunk number 13: fitDisp
###################################################
pasillaExons <- fitDispersionFunction(pasillaExons)
head(fData(pasillaExons)$dispBeforeSharing)
pasillaExons@dispFitCoefs
head(fData(pasillaExons)$dispFitted)
head(fData(pasillaExons)$dispersion)


###################################################
### code chunk number 14: fitdiagnostics
###################################################
meanvalues <- rowMeans(counts(pasillaExons))
plot(meanvalues, fData(pasillaExons)$dispBeforeSharing, log="xy", main="mean vs CR dispersion")
x <- 0.01:max(meanvalues)
y <- pasillaExons@dispFitCoefs[1] + pasillaExons@dispFitCoefs[2] / x
lines(x, y, col="red")


###################################################
### code chunk number 15: mffg2
###################################################
testGeneForDEU( pasillaExons, "FBgn0010909" )


###################################################
### code chunk number 16: mffg3
###################################################
tgdeu = testGeneForDEU( pasillaExons, "FBgn0010909" )
specialExon = "E010"
stopifnot(tgdeu[specialExon,"pvalue"]<1e-7)


###################################################
### code chunk number 17: testForDEU1
###################################################
pasillaExons <- testForDEU( pasillaExons )


###################################################
### code chunk number 18: testForDEU2
###################################################
pasillaExons <- estimatelog2FoldChanges( pasillaExons )
res1 <- DEUresultTable(pasillaExons)


###################################################
### code chunk number 19: testForDEU3
###################################################
table ( res1$padjust < 0.1 )


###################################################
### code chunk number 20: MvsA
###################################################
plot(res1$meanBase, res1[,"log2fold(untreated/treated)"], log="x", col=ifelse(res1$padjust < 0.1, "red", "black"), ylim=c(-4,4), main="pasilla MvsA")


###################################################
### code chunk number 21: design
###################################################
design(pasillaExons)


###################################################
### code chunk number 22: formuladispersion
###################################################
formuladispersion <- count ~ sample + ( exon + type ) * condition
pasillaExons <- estimateDispersions( pasillaExons, formula = formuladispersion )
pasillaExons <- fitDispersionFunction(pasillaExons)


###################################################
### code chunk number 23: formula1
###################################################
formula0 <- count ~ sample + type * exon + condition
formula1 <- count ~ sample + type * exon + condition * I(exon == exonID)


###################################################
### code chunk number 24: formula2
###################################################
pasillaExons <- testForDEU( pasillaExons, formula0=formula0, formula1=formula1 )
res2 <- DEUresultTable( pasillaExons )


###################################################
### code chunk number 25: formula3
###################################################
table( res2$padjust < 0.1 )


###################################################
### code chunk number 26: figcomparep
###################################################
bottom = function(x) pmax(x, 1e-6)
plot(bottom(res1$padjust), bottom(res2$padjust), log="xy")


###################################################
### code chunk number 27: plot1
###################################################
plotDEXSeq(pasillaExons, "FBgn0010909", cex.axis=1.2, cex=1.3, lwd=2, legend=TRUE)


###################################################
### code chunk number 28: checkClaim
###################################################
wh = (fData(pasillaExons)$geneID=="FBgn0010909")
stopifnot(sum(fData(pasillaExons)$padjust[wh] < formals(plotDEXSeq)$FDR)==1)


###################################################
### code chunk number 29: plot2
###################################################
plotDEXSeq(pasillaExons, "FBgn0010909", displayTranscripts=TRUE, cex.axis=1.2, cex=1.3, lwd=2, legend=TRUE)


###################################################
### code chunk number 30: plot3
###################################################
plotDEXSeq(pasillaExons, "FBgn0010909", expression=FALSE, norCounts=TRUE, cex.axis=1.2, cex=1.3, lwd=2, legend=TRUE)


###################################################
### code chunk number 31: plot4
###################################################
plotDEXSeq(pasillaExons, "FBgn0010909", cex.axis=1.2, cex=1.3, lwd=2, legend=TRUE, expression=FALSE, splicing=TRUE)


###################################################
### code chunk number 32: DEXSeqHTML (eval = FALSE)
###################################################
## DEXSeqHTML( pasillaExons, FDR=0.1, color=c("#FF000080", "#0000FF80") )


###################################################
### code chunk number 33: para1 (eval = FALSE)
###################################################
## data("pasillaExons", package="pasilla")
## library(multicore)
## pasillaExons <- estimateSizeFactors( pasillaExons )
## pasillaExons <- estimateDispersions( pasillaExons, nCores=3, quiet=TRUE)
## pasillaExons <- fitDispersionFunction( pasillaExons )
## pasillaExons <- testForDEU( pasillaExons, nCores=3)


###################################################
### code chunk number 34: alldeu
###################################################
data("pasillaExons", package="pasilla")
pasillaExons <- makeCompleteDEUAnalysis(pasillaExons, formulaDispersion=formuladispersion, formula0=formula0, formula1=formula1, nCores=1)


###################################################
### code chunk number 35: dirpasilla
###################################################
dir(system.file("extdata",package="pasilla"))


###################################################
### code chunk number 36: ecswithout
###################################################
bare <- newExonCountSet(
   countData = counts(pasillaExons), 
   design=design(pasillaExons), 
   geneIDs=geneIDs(pasillaExons), 
   exonIDs=exonIDs(pasillaExons))


###################################################
### code chunk number 37: gct
###################################################
head(geneCountTable(pasillaExons))


###################################################
### code chunk number 38: sessionInfo
###################################################
sessionInfo()