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

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


###################################################
### code chunk number 2: start
###################################################
library("DEXSeq")
library("pasilla")
data("pasillaExons", package="pasilla")


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


###################################################
### code chunk number 4: counts
###################################################
head( counts(pasillaExons), 20 )


###################################################
### code chunk number 5: options
###################################################
options(width=96)


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


###################################################
### code chunk number 7: tabtab1
###################################################
tabGeneIDsPasillaExons = table(geneIDs(pasillaExons))
ngenesPasillaExons = sum(tabGeneIDsPasillaExons>0)
tabtabGeneIDsPasillaExons = table(tabGeneIDsPasillaExons)
stopifnot(tabtabGeneIDsPasillaExons["36"]==1,
          tabtabGeneIDsPasillaExons["16"]==3)


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


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


###################################################
### code chunk number 10: helperfunctions
###################################################
plotDispEsts = function( cds, ymin, linecol="#ff000080",
  xlab = "mean of normalized counts", ylab = "dispersion",
  log = "xy", cex = 0.45, ... )
{
  px = rowMeans( counts( cds, normalized=TRUE ) )
  sel = (px>0)
  px = px[sel]

  py = fData(cds)$dispBeforeSharing[sel]
  if(missing(ymin))
      ymin = 10^floor(log10(min(py[py>0], na.rm=TRUE))-0.1)

  plot(px, pmax(py, ymin), xlab=xlab, ylab=ylab,
    log=log, pch=ifelse(py<ymin, 6, 16), cex=cex, ... )
  xg = 10^seq( -.5, 5, length.out=100 )
  fun = function(x) { cds@dispFitCoefs[1] + cds@dispFitCoefs[2] / x }
  lines( xg, fun(xg), col=linecol, lwd=4)
}

plotMA = function(x, ylim,
  col = ifelse(x$padj>=0.1, "gray32", "red3"),
  linecol = "#ff000080",
  xlab = "mean of normalized counts", ylab = expression(log[2]~fold~change),
  log = "x", cex=0.45, ...)
{
  if (!(is.data.frame(x) && all(c("baseMean", "log2FoldChange") %in% colnames(x))))
    stop("'x' must be a data frame with columns named 'baseMean', 'log2FoldChange'.")

  x = subset(x, baseMean!=0)
  py = x$log2FoldChange
  if(missing(ylim))
      ylim = c(-1,1) * quantile(abs(py[is.finite(py)]), probs=0.99) * 1.1
  plot(x$baseMean, pmax(ylim[1], pmin(ylim[2], py)),
       log=log, pch=ifelse(py<ylim[1], 6, ifelse(py>ylim[2], 2, 16)),
       cex=cex, col=col, xlab=xlab, ylab=ylab, ylim=ylim, ...)
  abline(h=0, lwd=4, col=linecol)
}


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


###################################################
### code chunk number 12: fitDisp1
###################################################
pasillaExons <- fitDispersionFunction(pasillaExons)


###################################################
### code chunk number 13: fitDisp2
###################################################
head(fData(pasillaExons)$dispBeforeSharing)
pasillaExons@dispFitCoefs
head(fData(pasillaExons)$dispFitted)


###################################################
### code chunk number 14: fitdiagnostics
###################################################
plotDispEsts( pasillaExons )


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


###################################################
### code chunk number 16: testForDEU2
###################################################
head( fData(pasillaExons)[, c( "pvalue", "padjust" ) ] )


###################################################
### code chunk number 17: estFC
###################################################
pasillaExons <- estimatelog2FoldChanges( pasillaExons )


###################################################
### code chunk number 18: DEUresults
###################################################
res1 <- DEUresultTable(pasillaExons)
head( res1 )


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


###################################################
### code chunk number 20: tallyGenes
###################################################
table ( tapply( res1$padjust < 0.1, geneIDs(pasillaExons), any ) )


###################################################
### code chunk number 21: MvsA
###################################################
plotMA(with(res1, data.frame(baseMean = meanBase,
                             log2FoldChange = `log2fold(untreated/treated)`,
                             padj = padjust)),
     ylim=c(-4,4), cex=0.8)


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


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


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


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


###################################################
### code chunk number 26: formula3
###################################################
table( res2$padjust < 0.1 )
table(res1$padjust < 0.1, res2$padjust < 0.1)


###################################################
### code chunk number 27: figcomparep
###################################################
bottom = function(x) pmax(x, 1e-6)
plot(bottom(res1$padjust), bottom(res2$padjust), log="xy", pch=20)
abline(a=0,b=1, col="red3")
abline(v=0.1, h=0.1, col="blue")


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


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


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


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


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


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


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


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


###################################################
### code chunk number 36: pkgRoot
###################################################
pkgDir = system.file(package="DEXSeq")
pkgDir
list.files(pkgDir)
list.files(file.path(pkgDir, "python_scripts"))


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


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


###################################################
### code chunk number 39: buildExonCountSet (eval = FALSE)
###################################################
## 
##   library(GenomicFeatures)
##   hse <- makeTranscriptDbFromBiomart(biomart="ensembl", dataset="hsapiens_gene_ensembl")
##   exonicParts <- prepareAnnotationForDEXSeq( hse, aggregateGenes=TRUE )


###################################################
### code chunk number 40: buildExonCountSet2 (eval = FALSE)
###################################################
##   
##   bamDir <- system.file("extdata",package="parathyroidSE",mustWork=TRUE)
##   fls <- list.files(bamDir, pattern="bam$",full=TRUE)
##   bamlst <- BamFileList(fls)
##   SE <- countReadsForDEXSeq( exonicParts, bamlst )


###################################################
### code chunk number 41: buildExonCountSet3 (eval = FALSE)
###################################################
## 
##   ecs <- buildExonCountSet( SE, c("A", "A", "B"), exonicParts )


###################################################
### code chunk number 42: countTableForGene
###################################################
head(countTableForGene(pasillaExons,"FBgn0010909"))


###################################################
### code chunk number 43: countTableForGeneNorm
###################################################
head(countTableForGene(pasillaExons,"FBgn0010909", normalized=TRUE))


###################################################
### code chunk number 44: modelFrame
###################################################
mf <- modelFrameForGene( pasillaExons, "FBgn0010909" )
head( mf )


###################################################
### code chunk number 45: mffg1
###################################################
mf <- estimateExonDispersionsForModelFrame( modelFrameForGene( pasillaExons, "FBgn0010909" ) )


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


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


###################################################
### code chunk number 48: acc
###################################################
head(geneIDs(pasillaExons))
head(exonIDs(pasillaExons))


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