### 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-87
###################################################
setwd(system.file("extdata",package="BitSeq"));


###################################################
### code chunk number 5: BitSeq.Rnw:98-104
###################################################
res1 <- getExpression("data-c0b0.sam",
   "ensSelect1.fasta",
   log = TRUE,
   MCMC_burnIn=200,
   MCMC_samplesN=200,
   MCMC_scaleReduction=1.1)


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


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


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


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


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


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


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


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


###################################################
### code chunk number 14: BitSeq.Rnw:184-193
###################################################
estimateExpression( "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 15: BitSeq.Rnw:205-206 (eval = FALSE)
###################################################
## loadSamples("data-c0b0.rpkm")


###################################################
### code chunk number 16: BitSeq.Rnw:216-218
###################################################
cond1Files = c("data-c0b0.rpkm","data-c0b1.rpkm")
cond2Files = c("data-c1b1.rpkm","data-c1b1.rpkm")


###################################################
### code chunk number 17: BitSeq.Rnw:222-225
###################################################
getMeanVariance(c(cond1Files,cond2Files), 
   outFile = "data.means",
   log = TRUE )


###################################################
### code chunk number 18: BitSeq.Rnw:232-237
###################################################
estimateHyperPar( outFile = "data.par",
   cond1 = cond1Files, 
   cond2 = cond2Files,
   meanFile = "data.means",
   verbose = TRUE )


###################################################
### code chunk number 19: BitSeq.Rnw:245-248
###################################################
estimateDE(cond1Files, cond2Files,
   outFile = "data",
   parFile = "data.par" )


###################################################
### code chunk number 20: BitSeq.Rnw:252-256
###################################################
estimateDE(cond1Files, cond2Files,
   outFile = "data",
   parFile = "data.par",
   samples = TRUE )


###################################################
### code chunk number 21: BitSeq.Rnw:269-277
###################################################
res1 <- getExpression("data-c0b0.sam",
   "ensSelect1.fasta",
   outPrefix="localDir/data-c0b0",
   log = TRUE,
   MCMC_burnIn=200,
   MCMC_samplesN=200,
   MCMC_scaleReduction=1.1,
   pretend=TRUE)