## ----setup, message=FALSE, warning=FALSE--------------------------------------
library(BatchVaria)
library(SummarizedExperiment)
library(BiocStyle)

## ----installation-release, eval=FALSE-----------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# 
# BiocManager::install("BatchVaria")

## ----installation-devel, eval=FALSE-------------------------------------------
# BiocManager::install(version = "devel")
# BiocManager::install("BatchVaria")

## ----runCorrection------------------------------------------------------------
## exampleBatchVaria() draws random data; seed for a reproducible vignette
set.seed(1)
bv <- exampleBatchVaria(confounding = 0.6)

class(bv)

table(
    batch = bv$batch,
    group = bv$group
)

## ----methods------------------------------------------------------------------
availableCorrectionMethods()

## ----corrections--------------------------------------------------------------
bv <- runCorrection(
    bv,
    method = "combat", batch = "batch",
    newAssayName = "naive"
)

bv <- runCorrection(
    bv,
    method = "combat", batch = "batch", preserve = "group",
    newAssayName = "protected"
)

assayNames(bv)

## ----profile------------------------------------------------------------------
assaysUsed <- c("raw", "naive", "protected", "raw_noise")

bv <- profileVariance(
    bv,
    formula = ~ batch + group,
    assays = assaysUsed,
    methods = c("anova", "variancePartition")
)

## ----ambiguous, error=TRUE----------------------------------------------------
try({
varianceTable(bv, assays = assaysUsed)
})

## ----provenance---------------------------------------------------------------
p <- provenance(bv)

## ----provenance-return--------------------------------------------------------
p$corrections[, c("method", "output", "preserve", "no_op")]

## ----table--------------------------------------------------------------------
vt <- varianceTable(bv, assays = assaysUsed, method = "anova")

## ----percent------------------------------------------------------------------
vt$percent

## ----delta--------------------------------------------------------------------
varianceTable(
    bv,
    assays = assaysUsed,
    method = "anova",
    formatDelta = TRUE
)$delta

## ----totalvar-----------------------------------------------------------------
av <- assayVariance(bv, assays = assaysUsed)
av

## ----absolute-----------------------------------------------------------------
vr <- varianceResults(bv, method = "anova")

absVar <- merge(
    vr[, c("assay", "term", "variance_fraction")],
    av[, c("assay", "total_variance")],
    by = "assay"
)
absVar$absolute <- round(absVar$variance_fraction * absVar$total_variance, 1)

## as.data.frame(), because pivot_wider() returns a tibble and tibble
## prints three significant figure.
as.data.frame(tidyr::pivot_wider(
    absVar[, c("assay", "term", "absolute")],
    names_from = "term",
    values_from = "absolute"
))

## ----concordance--------------------------------------------------------------
cmp <- varianceResults(bv, assayName = "raw")[
    , c("method", "term", "variance_fraction")
]

as.data.frame(
    tidyr::pivot_wider(
        cmp,
        names_from = "method",
        values_from = "variance_fraction"
    )
)

## ----orthogonal---------------------------------------------------------------
orth <- profileVariance(
    exampleBatchVaria(confounding = 0),
    ~ batch + group,
    assays = "raw",
    methods = c("anova", "variancePartition")
)

## Printed at full precision deliberately: at three significant figures
## both columns read 0.163 and the agreement below cannot be checked.
as.data.frame(
    tidyr::pivot_wider(
        varianceResults(orth, assayName = "raw")[
            , c("method", "term", "variance_fraction")
        ],
        names_from = "method",
        values_from = "variance_fraction"
    )
)

## ----retention----------------------------------------------------------------
basisRetention(bv, assays = assaysUsed)

## ----summaryfig, fig.width=10, fig.height=6, out.width="100%"-----------------
ev <- evaluateCorrections(
    bv,
    baseline = "raw",
    assays = assaysUsed,
    diagnostics = "variance"
)

plotVarianceDelta(ev)

## ----evalall------------------------------------------------------------------
ev <- evaluateCorrections(
    bv,
    baseline = "raw",
    assays = assaysUsed,
    varianceMethod = "anova"
)

names(ev)

## ----composition, fig.width=8, fig.height=5, out.width="100%"-----------------
plotVarianceComposition(bv, assays = assaysUsed, method = "anova")

## ----radar, fig.width=7, fig.height=6, out.width="80%"------------------------
plotVarianceRadar(bv, assays = assaysUsed, method = "anova")

## ----basisprojection----------------------------------------------------------
coords <- basisProjection(
    bv,
    assays = c("raw", "naive"),
    reference = "raw"
)

head(coords)

## ----basisplot, fig.width=6, fig.height=4.5, out.width="70%", fig.show="hold"----
for (p in plotBasisProjection(
    bv,
    assays = c("raw", "naive", "protected"),
    reference = "raw"
)) {
    print(p)
}

## ----basisplot-group, fig.width=6, fig.height=4.5, out.width="70%"------------
plotBasisProjection(
    bv,
    assays = "naive",
    reference = "raw",
    colourBy = "group"
)[[1]]

## ----trajectory, fig.width=7, fig.height=5, out.width="80%"-------------------
plotPCATrajectory(bv, assayBefore = "raw", assayAfter = "naive")

## ----distance, fig.width=6.5, fig.height=5.5, out.width="70%"-----------------
plotSampleDistance(bv, assayName = "raw")

## ----corrchange, fig.width=7, fig.height=4.5, out.width="80%"-----------------
plotCorrelationChange(ev)

## ----entropy, fig.width=7, fig.height=4, out.width="80%", fig.show="hold"-----
plotBatchEntropy(bv, assayName = "raw", batchVar = "batch", k = 10)
plotBatchEntropy(bv, assayName = "naive", batchVar = "batch", k = 10)

## ----session------------------------------------------------------------------
sessionInfo()