## ---- echo=FALSE, results="hide", message=TRUE---------------------------
require(knitr)
opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)

## ----style, echo=FALSE, results='asis'-----------------------------------
BiocStyle::markdown()

## ----library, echo=FALSE-------------------------------------------------
library(BASiCS)

## ----ExampleDataTest-----------------------------------------------------
set.seed(1)
Counts = Counts = matrix(rpois(50*40, 2), ncol = 40)
rownames(Counts) <- c(paste0("Gene", 1:40), paste0("Spike", 1:10))
Tech = c(rep(FALSE,40),rep(TRUE,10))
set.seed(2)
SpikeInput = rgamma(10,1,1)
SpikeInfo <- data.frame("SpikeID" = paste0("Spike", 1:10), 
                        "SpikeInput" = SpikeInput)

# No batch structure
DataExample = newBASiCS_Data(Counts, Tech, SpikeInfo)

# With batch structure
DataExample = newBASiCS_Data(Counts, Tech, SpikeInfo, 
                             BatchInfo = rep(c(1,2), each = 20)) 

## ----quick-start-MCMC----------------------------------------------------
Data <- makeExampleBASiCS_Data()
Chain <- BASiCS_MCMC(Data = Data, N = 1000, Thin = 10, Burn = 500, 
                     PrintProgress = FALSE)

## ----LoadSingleData------------------------------------------------------
data(ChainSC)

## ----quick-start-HVGdetection, fig.height = 6, fig.width = 6-------------
par(mfrow = c(2,2))
HVG <- BASiCS_DetectHVG(ChainSC, VarThreshold = 0.6, Plot = TRUE)
LVG <- BASiCS_DetectLVG(ChainSC, VarThreshold = 0.2, Plot = TRUE)

## ----quick-start-HVGdetectionTable---------------------------------------
head(HVG$Table)
head(LVG$Table)

## ----quick-start-HVGdetectionPlot----------------------------------------
SummarySC <- Summary(ChainSC)
plot(SummarySC, Param = "mu", Param2 = "delta", log = "xy")
with(HVG$Table[HVG$Table$HVG == TRUE,], points(Mu, Delta))
with(LVG$Table[LVG$Table$LVG == TRUE,], points(Mu, Delta))

## ----quick-start-LoadBothData--------------------------------------------
data(ChainSC)
data(ChainRNA)

## ----quick-start-TestDE, fig.width=10, fig.height=5----------------------
Test <- BASiCS_TestDE(Chain1 = ChainSC, Chain2 = ChainRNA,
                      GroupLabel1 = "SC", GroupLabel2 = "PaS",
                      EpsilonM = log2(1.5), EpsilonD = log2(1.5),
                      EFDR_M = 0.10, EFDR_D = 0.10,
                      Offset = TRUE, OffsetPlot = TRUE, Plot = TRUE)

## ----quick-start-DE------------------------------------------------------
head(Test$TableMean)
head(Test$TableDisp)

## ----quick-start-TestDE-2, fig.width=10, fig.height=5--------------------
Test <- BASiCS_TestDE(Chain1 = ChainSC, Chain2 = ChainRNA,
                      GroupLabel1 = "SC", GroupLabel2 = "PaS",
                      EpsilonM = 0, EpsilonD = log2(1.5),
                      EFDR_M = 0.10, EFDR_D = 0.10,
                      Offset = TRUE, OffsetPlot = TRUE, Plot = TRUE)

## ----MCMCNotRun----------------------------------------------------------
Data <- makeExampleBASiCS_Data()
Chain <- BASiCS_MCMC(Data, N = 1000, Thin = 10, Burn = 500, 
                     PrintProgress = FALSE, StoreChains = TRUE, 
                     StoreDir = tempdir(), RunName = "Example")

## ----LoadChainNotRun-----------------------------------------------------
Chain <- BASiCS_LoadChain("Example", StoreDir = tempdir()) 

## ----Traceplots, fig.height = 3, fig.width = 6---------------------------
plot(Chain, Param = "mu", Gene = 1, log = "y")
plot(Chain, Param = "phi", Cell = 1)

## ----AccessChains--------------------------------------------------------
displayChainBASiCS(Chain, Param = "mu")[1:5,1:5]

## ----Summary-------------------------------------------------------------
ChainSummary <- Summary(Chain)

## ----SummaryExample------------------------------------------------------
head(displaySummaryBASiCS(ChainSummary, Param = "mu"))

## ----OtherHPD, fig.width = 7, fig.height = 7-----------------------------
par(mfrow = c(2,2))
plot(ChainSummary, Param = "mu", main = "All genes", log = "y")
plot(ChainSummary, Param = "mu", Genes = 1:10, main = "First 10 genes")
plot(ChainSummary, Param = "delta", main = "All genes")
plot(ChainSummary, Param = "delta", Genes = c(2,5,10,50), main = "5 customized genes")

## ----Normalisation, fig.width = 7, fig.height = 3.5----------------------
par(mfrow = c(1,2))
plot(ChainSummary, Param = "phi")
plot(ChainSummary, Param = "s", Cells = 1:5)

## ----DispVsExp, fig.width = 7, fig.height = 3.5--------------------------
par(mfrow = c(1,2))
plot(ChainSummary, Param = "mu", Param2 = "delta", log = "x", SmoothPlot = FALSE)
plot(ChainSummary, Param = "mu", Param2 = "delta", log = "x", SmoothPlot = TRUE)

## ----DenoisedCounts------------------------------------------------------
DenoisedCounts = BASiCS_DenoisedCounts(Data = Data, Chain = Chain)
DenoisedCounts[1:5, 1:5]

## ----DenoisedProp--------------------------------------------------------
DenoisedRates <- BASiCS_DenoisedRates(Data = Data, Chain = Chain, 
                                      Propensities = FALSE)
DenoisedRates[1:5, 1:5]

## ----DenoisedRates-------------------------------------------------------
DenoisedProp = BASiCS_DenoisedRates(Data = Data, Chain = Chain, 
                                    Propensities = TRUE)
DenoisedProp[1:5, 1:5]

## ----SessionInfo---------------------------------------------------------
sessionInfo()