## ----style-knitr, eval=TRUE, echo=FALSE, results="asis"---------------------------------
BiocStyle::latex()

## ----load_library, message=FALSE, eval=TRUE, echo=FALSE---------------------------------
library(INSPEcT)

## ----features_quantification, message=FALSE, eval=TRUE----------------------------------
require(TxDb.Mmusculus.UCSC.mm9.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm9.knownGene
paths_4su <- system.file('extdata', '4sURNA_0h.bam', package="INSPEcT")
paths_total <- system.file('extdata', 'totalRNA_0h.bam', package="INSPEcT")
makeRPKMsOut <- makeRPKMs(txdb, paths_4su, paths_total)
rpkms <- makeRPKMsOut$rpkms

## ----rate_estimation, message=TRUE------------------------------------------------------
data('rpkms', package='INSPEcT')
tpts <- c(0, 1/6, 1/3, 1/2, 1, 2, 4, 8, 16)
tL <- 1/6
mycerIds <- newINSPEcT(tpts, tL, rpkms$foursu_exons, rpkms$total_exons, 
	rpkms$foursu_introns, rpkms$total_introns, BPPARAM=SerialParam())
head(ratesFirstGuess(mycerIds, 'synthesis'))

## ----rate_estimation_ddp, message=FALSE-------------------------------------------------
mycerIdsDdp <- newINSPEcT(tpts, tL, rpkms$foursu_exons, rpkms$total_exons, 
	rpkms$foursu_introns, rpkms$total_introns, BPPARAM=SerialParam(), degDuringPulse=TRUE)
head(ratesFirstGuess(mycerIdsDdp, 'synthesis'))

## ----pre_model_heatmap, message=FALSE, fig.width=6, fig.height=4, fig.cap="Heatmap of the pre-model rates representing the concentrations of total mRNA and pre-mRNA and the rates of the mRNA life cycle", fig.align="center"----
mycerIds10 <- mycerIds[1:10]
inHeatmap(mycerIds10, clustering=FALSE)

## ----model_rates, message=TRUE, eval=FALSE----------------------------------------------
#  mycerIds10 <- modelRates(mycerIds10, seed=1)

## ----post_model_heatmap, message=TRUE, fig.width=6, fig.height=4, fig.cap="Heatmap of the modeled rates representing the concentrations of total mRNA and pre-mRNA and the rates of the mRNA life cycle", fig.align="center", eval=TRUE----
data('mycerIds10', package='INSPEcT')
head(viewModelRates(mycerIds10, 'synthesis'))
inHeatmap(mycerIds10, type='model', clustering=FALSE)

## ----gene_class, message=TRUE-----------------------------------------------------------
geneClass(mycerIds10)

## ----gene_plot, message=TRUE, fig.width=10, fig.height=3, fig.cap="Pre-model and modeled concentrations and rates for a selected gene with standard deviation, where available", fig.align="center"----
plotGene(mycerIds10, 2, fix.yaxis=TRUE)

## ----goodness_of_fit, message=TRUE, fig.width=3.5, fig.height=3.5, hold=TRUE, fig.cap="Histogram of the p-values from the goodness of fit test for selected models", fig.align="center"----
chisq <- chisqmodel(mycerIds10)
hist(log10(chisq), main='', xlab='log10 chi-squared p-value')
discard <- which(chisq>.05)
featureNames(mycerIds10)[discard]
mycerIds10new <- mycerIds10[-discard]

## ----simulated_data_notrun, message=FALSE, eval=FALSE-----------------------------------
#  simRates <- makeSimModel(mycerIds, 1000, newTpts=NULL, seed=1)

## ----simulated_data_notrun2, message=FALSE, eval=FALSE----------------------------------
#  simData1rep <- makeSimDataset(simRates, tpts, 1, seed=1)
#  simData1rep <- modelRates(simData1rep, seed=1)
#  newTpts <- c(0, 1/6, 1/3, 1/2, 1, 1.5, 2, 4, 8, 12, 16, 24)
#  simData3rep <- makeSimDataset(simRates, newTpts, 3, seed=1)
#  simData3rep <- modelRates(simData3rep, seed=1)

## ----simulated_data_run, message=FALSE, eval=TRUE, fig.width=10, fig.height=5, fig.cap="For each rate, INSPEcT classification performance is measured in terms of sensitivity and specificity using a ROC curve analysis (rocCurve method). False negatives (FN) represent cases where the rate is identified as constant while it was simulated as varying. False positives (FP) represent cases where INSPEcT identified a rate as varying while it was simulated as constant. On the contrary, true positives (TP) and negatives (TN) are cases of correct classification of varying and constant rates, respectively. Consequently, sensitivity and specificity are computed using increasing thresholds for the Brown's method used to combine multiple p-values derived from the log-likelihood ratio tests", fig.align="center"----
data('simRates', package='INSPEcT')
data('simData1rep', package='INSPEcT')
data('simData3rep', package='INSPEcT')
par(mfrow=c(1,2))
rocCurve(simRates, simData1rep); title("1 replicate - 9 time points", line=3)
rocCurve(simRates, simData3rep); title("3 replicates - 12 time points", line=3)

## ----simulated_data_run2, message=FALSE, eval=TRUE, fig.width=8, fig.height=4, fig.cap="Plot of the sensitivity (black curve) and specificity (red curve) that is achieved after performing the log-likelihood ratio and Brown's method for combining p-values with selected thresholds. Thresholds that can be set for chi-squared test to accept models that will undergo the log-likelihood ratio test and for Brown's p-value to assess variability of rates."----
rocThresholds(simRates, simData3rep, bTsh=c(.01,.01,.05), cTsh=.1)

## ----params1, message=FALSE, eval=FALSE-------------------------------------------------
#  modelingParams(mycerIds10)$useSigmoidFun <- FALSE
#  modelingParams(mycerIds10)$nInit <- 20
#  modelingParams(mycerIds10)$nIter <- 1000

## ----params2, message=FALSE, eval=FALSE-------------------------------------------------
#  thresholds(mycerIds10)$chisquare <- .1
#  thresholds(mycerIds10)$brown <- c(alpha=.01, beta=.01, gamma=.05)
#  llrtests(mycerIds10)$processing <- list(c('0','c'))

## ----params3, message=FALSE, eval=FALSE-------------------------------------------------
#  ## modeling
#  modelingParams(mycerIds10)
#  ## model selection and testing framework
#  modelSelection(mycerIds10)
#  thresholds(mycerIds10)
#  llrtests(mycerIds10)

## ----steady1, message=FALSE, eval=TRUE--------------------------------------------------
newTpts <- c(0, 1/6)
simData3rep_2 <- makeSimDataset(simRates, newTpts, 3, seed=2)

## ----steady2, message=FALSE, eval=TRUE--------------------------------------------------
diffrates <- compareSteady(simData3rep, simData3rep_2)

## ----steady3, message=TRUE, eval=TRUE---------------------------------------------------
diffrates

## ----steady4, message=TRUE, eval=FALSE--------------------------------------------------
#  head(synthesis(diffrates))
#  head(processing(diffrates))
#  head(degradation(diffrates))

## ----steady5, message=TRUE, eval=TRUE, fig.width=5, fig.height=5, hold=TRUE, fig.cap="Example of the plotMA generated image. Orange triangles correspond to genes whose rates are differentially used between the two conditions, blue cloud correnspond to the whole ditribution of rates.", fig.align="center"----
plotMA(diffrates, rate='synthesis', alpha=.5)

## ----r----------------------------------------------------------------------------------
sessionInfo()