### R code from vignette source 'DEXSeq'

###################################################
### 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: pData
###################################################
pData(pasillaExons)


###################################################
### code chunk number 5: fData
###################################################
head(fData(pasillaExons)[,c(1,2,7:10)])


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


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


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


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


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


###################################################
### code chunk number 11: estDisp2
###################################################
head(fData(pasillaExons)$dispersion_CR_est)
pasillaExons@dispFitCoefs
head(fData(pasillaExons)$dispersion)


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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