## ---- echo = FALSE, results = "hide", message = FALSE-------------------------
require(knitr)
opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE, dpi = 30)
knitr::opts_chunk$set(dev="png")

## ----library, echo = FALSE----------------------------------------------------
library("BASiCS")
library("BiocStyle")
library("SingleCellExperiment")
library("cowplot")
library("ggplot2")
theme_set(theme_bw())

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

## ----ExampleDataTest----------------------------------------------------------
set.seed(1)
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)) 

## ----ExampleDataNoSpikes------------------------------------------------------
set.seed(1)
CountsNoSpikes <- matrix(rpois(50*40, 2), ncol = 40)
rownames(CountsNoSpikes) <- paste0("Gene", 1:50)

# With batch structure
DataExampleNoSpikes <- SingleCellExperiment(assays = list(counts = CountsNoSpikes), 
                      colData = data.frame(BatchInfo = rep(c(1,2), each = 20))) 

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

## ----quick-start-HVGdetection, fig.height = 8, fig.width = 20-----------------
HVG <- BASiCS_DetectHVG(ChainSC, VarThreshold = 0.6, Plot = FALSE)
LVG <- BASiCS_DetectLVG(ChainSC, VarThreshold = 0.2, Plot = FALSE)

plot_grid(
  BASiCS_PlotVG(HVG, "Grid"),
  BASiCS_PlotVG(HVG, "VG")
)
plot_grid(
  BASiCS_PlotVG(LVG, "Grid"),
  BASiCS_PlotVG(LVG, "VG")
)

## ----quick-start-HVGdetectionTable--------------------------------------------
format(HVG)
format(LVG)

## ----quick-start-HVGdetectionPlot, fig.width = 8, fig.height = 8, eval = FALSE----
#  SummarySC <- Summary(ChainSC)
#  plot(SummarySC, Param = "mu", Param2 = "delta", log = "xy")
#  HTable <- format(HVG)
#  LTable <- format(LVG)
#  with(HTable, points(Mu, Delta))
#  with(LTable, points(Mu, Delta))

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

## ----quick-start-TestDE, fig.width = 16, fig.height = 8-----------------------
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, PlotOffset = TRUE, Plot = TRUE)

## ----testing, message = TRUE--------------------------------------------------
Test

## ----format-------------------------------------------------------------------
head(format(Test, Which = "Mean"))
head(format(Test, Which = "Disp"))
## This object doesn't contain residual overdispersion tests as the chains
## were run using the non-regression version of BASiCS
# head(format(DE, Which = "Disp"))

## ----rowData------------------------------------------------------------------
rowData(Test)
rowData(Test) <- cbind(rowData(Test), Index = 1:nrow(rowData(Test)))
format(Test, "Mean")

## ----plots, fig.height=16, fig.width=20---------------------------------------
BASiCS_PlotDE(Test)
BASiCS_PlotDE(Test, Plots = c("MA", "Volcano"))
BASiCS_PlotDE(Test, Plots = "MA", Parameters = "Mean")

## ----quick-start-TestDE-2, fig.width = 16, fig.height = 8---------------------
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, PlotOffset = FALSE, Plot = FALSE)

## ---- fig.height = 8, fig.width = 8-------------------------------------------
DataNoSpikes <- newBASiCS_Data(Counts, Tech, SpikeInfo = NULL, 
                              BatchInfo = rep(c(1,2), each = 20))

# Alternatively
DataNoSpikes <- SingleCellExperiment(assays = list(counts = Counts), 
                      colData = data.frame(BatchInfo = rep(c(1,2), each = 20))) 

ChainNoSpikes <- BASiCS_MCMC(Data = DataNoSpikes, N = 1000, 
                             Thin = 10, Burn = 500, 
                             WithSpikes = FALSE,  Regression = TRUE,
                             PrintProgress = FALSE)

## ---- fig.height = 8, fig.width = 8-------------------------------------------
DataRegression <- newBASiCS_Data(Counts, Tech, SpikeInfo, 
                              BatchInfo = rep(c(1,2), each = 20))
ChainRegression <- BASiCS_MCMC(Data = DataRegression, N = 1000, 
                               Thin = 10, Burn = 500, 
                               Regression = TRUE, 
                               PrintProgress = FALSE)

## ---- fig.height = 12, fig.width = 12-----------------------------------------
data("ChainRNAReg")
BASiCS_ShowFit(ChainRNAReg)

## ---- fig.height = 8, fig.width = 16------------------------------------------
data("ChainSCReg")
Test <- BASiCS_TestDE(Chain1 = ChainSCReg, Chain2 = ChainRNAReg,
                      GroupLabel1 = "SC", GroupLabel2 = "PaS",
                      EpsilonM = log2(1.5), EpsilonD = log2(1.5),
                      EpsilonR = log2(1.5)/log2(exp(1)),
                      EFDR_M = 0.10, EFDR_D = 0.10,
                      Offset = TRUE, PlotOffset = FALSE, Plot = FALSE)

## ----testde-resdisp, fig.width=24, fig.height=16------------------------------
head(format(Test, Which = "ResDisp"))
BASiCS_PlotDE(Test, Parameters = "ResDisp")

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

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

## ----Traceplots, fig.height = 8, fig.width = 16-------------------------------
plot(Chain, Param = "mu", Gene = 1, log = "y")

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

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

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

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

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

## ----DispVsExp, fig.width = 16, fig.height = 8--------------------------------
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:2, 1:2]

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

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

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