inst/doc/BitSeq.R

In BitSeq: Transcript expression inference and differential expression analysis for RNA-seq data

### R code from vignette source 'BitSeq.Rnw'

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: BitSeq.Rnw:30-31
###################################################
options(width = 60)


###################################################
### code chunk number 3: BitSeq.Rnw:51-54 (eval = FALSE)
###################################################
## if (!requireNamespace("BiocManager", quietly=TRUE))
##     install.packages("BiocManager")
## BiocManager::install("BitSeq")


###################################################
### code chunk number 4: BitSeq.Rnw:58-59
###################################################
library(BitSeq)


###################################################
### code chunk number 5: BitSeq.Rnw:88-94
###################################################
# 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 6: BitSeq.Rnw:107-114
###################################################
res1 <- getExpression("data-c0b0.sam",
   "ensSelect1.fasta",
   log = TRUE,
   MCMC_burnIn=200,
   MCMC_samplesN=200,
   MCMC_samplesSave=50,
   seed=47)


###################################################
### code chunk number 7: BitSeq.Rnw:121-122 (eval = FALSE)
###################################################
## hist(res1$exp$mean)


###################################################
### code chunk number 8: BitSeq.Rnw:126-127
###################################################
samples1 <- loadSamples(res1$fn)


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


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


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


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


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


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


###################################################
### code chunk number 15: 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 16: 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 17: BitSeq.Rnw:236-237 (eval = FALSE)
###################################################
## loadSamples("data-c0b0.rpkm")


###################################################
### code chunk number 18: BitSeq.Rnw:250-254
###################################################
estimateVBExpression( "data-c0b0.prob", 
   outFile = "data-c0b0-vb", 
   outputType = "RPKM",
   trInfoFile = "data.tr")


###################################################
### code chunk number 19: BitSeq.Rnw:266-268
###################################################
allConditions = list(c("data-c0b0.rpkm","data-c0b1.rpkm"),
                     c("data-c1b1.rpkm","data-c1b1.rpkm"))


###################################################
### code chunk number 20: BitSeq.Rnw:272-273
###################################################
getMeanVariance(allConditions, outFile = "data.means", log = TRUE )


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


###################################################
### code chunk number 22: BitSeq.Rnw:292-302
###################################################
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 23: BitSeq.Rnw:306-308
###################################################
estimateDE(allConditions, outFile = "data", parFile = "data.par",
   samples = TRUE )


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


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