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

###################################################
### code chunk number 1: BitSeq.Rnw:29-30
###################################################
options(width = 40)


###################################################
### code chunk number 2: BitSeq.Rnw:50-52 (eval = FALSE)
###################################################
## source("http://www.bioconductor.org/biocLite.R")
## biocLite("BitSeq")


###################################################
### code chunk number 3: BitSeq.Rnw:56-57
###################################################
library(BitSeq)


###################################################
### code chunk number 4: BitSeq.Rnw:86-92
###################################################
# save the current directory
# (we move back to old_directory at the end of vignette)
old_directory = getwd();
on.exit(setwd(old_directory))
# move to directory with the data
setwd(system.file("extdata",package="BitSeq"));


###################################################
### code chunk number 5: BitSeq.Rnw:105-112
###################################################
res1 <- getExpression("data-c0b0.sam",
   "ensSelect1.fasta",
   log = TRUE,
   MCMC_burnIn=200,
   MCMC_samplesN=200,
   MCMC_samplesSave=50,
   seed=47)


###################################################
### code chunk number 6: BitSeq.Rnw:119-120 (eval = FALSE)
###################################################
## hist(res1$exp$mean)


###################################################
### code chunk number 7: BitSeq.Rnw:124-125
###################################################
samples1 <- loadSamples(res1$fn)


###################################################
### code chunk number 8: BitSeq.Rnw:128-130 (eval = FALSE)
###################################################
## plot( unlist(s2["ENST00000436661",]), 
##    unlist(s2["ENST00000373501",]))


###################################################
### code chunk number 9: BitSeq.Rnw:141-143
###################################################
cond1Files = c(res1$fn,"data-c0b1.rpkm")
cond2Files = c("data-c1b1.rpkm","data-c1b1.rpkm")


###################################################
### code chunk number 10: BitSeq.Rnw:149-152
###################################################
de1 <- getDE(list(cond1Files, cond2Files), 
   samples=TRUE)
print(de1$fn)


###################################################
### code chunk number 11: BitSeq.Rnw:156-157
###################################################
head( de1$pplr[ order(abs(0.5-de1$pplr$pplr), decreasing=TRUE ), ], 5)


###################################################
### code chunk number 12: BitSeq.Rnw:166-167
###################################################
setwd(system.file("extdata",package="BitSeq"))


###################################################
### code chunk number 13: BitSeq.Rnw:177-182
###################################################
parseAlignment( "data-c0b0.sam", 
   outFile = "data-c0b0.prob", 
   trSeqFile = "ensSelect1.fasta",
   trInfoFile = "data.tr",
   verbose = TRUE )


###################################################
### code chunk number 14: BitSeq.Rnw:196-204
###################################################
estimateExpression( "data-c0b0.prob", 
   outFile = "data-c0b0", 
   outputType = "RPKM",
   trInfoFile = "data.tr",
   MCMC_burnIn = 200,
   MCMC_samplesN = 200,
   MCMC_samplesSave = 100,
   MCMC_chainsN = 2 )


###################################################
### code chunk number 15: BitSeq.Rnw:220-229
###################################################
estimateExpressionLegacy( "data-c0b0.prob", 
   outFile = "data-c0b0", 
   outputType = "RPKM",
   trInfoFile = "data.tr",
   MCMC_burnIn = 200,
   MCMC_samplesN = 200,
   MCMC_samplesSave = 100,
   MCMC_scaleReduction = 1.1,
   MCMC_chainsN = 2 )


###################################################
### code chunk number 16: BitSeq.Rnw:235-236 (eval = FALSE)
###################################################
## loadSamples("data-c0b0.rpkm")


###################################################
### code chunk number 17: BitSeq.Rnw:246-248
###################################################
allConditions = list(c("data-c0b0.rpkm","data-c0b1.rpkm"),
                     c("data-c1b1.rpkm","data-c1b1.rpkm"))


###################################################
### code chunk number 18: BitSeq.Rnw:252-255
###################################################
getMeanVariance(allConditions, 
   outFile = "data.means",
   log = TRUE )


###################################################
### code chunk number 19: BitSeq.Rnw:262-266
###################################################
estimateHyperPar( outFile = "data.par",
   conditions = allConditions,
   meanFile = "data.means",
   verbose = TRUE )


###################################################
### code chunk number 20: BitSeq.Rnw:274-286
###################################################
estimateDE(allConditions,
   outFile = "data",
   parFile = "data.par" )
##
## pretend run with three conditions and normalization constants
##
cond3Files = c("data-c2b0.rpkm","data-c2b1.rpkm", "data-c2b2.rpkm")
estimateDE(list( allConditions[[1]], allConditions[[2]], cond3Files), 
           outFile="mydata", 
           parFile="mydata.par", 
           norm=c(1.0, 0.999, 1.0017, 0.9998, 1.0, 0.99, 0.97), 
           pretend=TRUE)


###################################################
### code chunk number 21: BitSeq.Rnw:290-294
###################################################
estimateDE(allConditions,
   outFile = "data",
   parFile = "data.par",
   samples = TRUE )


###################################################
### code chunk number 22: BitSeq.Rnw:306-316
###################################################
res1 <- getExpression("data-c0b0.sam",
   "ensSelect1.fasta",
   outPrefix="localDir/data-c0b0",
   log = TRUE,
   MCMC_burnIn=200,
   MCMC_samplesN=200,
   pretend=TRUE)

# restore the original working directory
setwd(old_directory);