## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----install, eval=FALSE------------------------------------------------------
# ## You can install the released version of barbieQ like so:
# # if (!require("BiocManager", quietly = TRUE))
# #     install.packages("BiocManager")
# #
# # BiocManager::install("barbieQ")
# 
# ## Alternatively, you can install the development version of barbieQ from GitHub
# devtools::install_github("Oshlack/barbieQ")

## ----library------------------------------------------------------------------
suppressPackageStartupMessages({
  library(barbieQ)
  library(magrittr)
  library(tidyr)
  library(dplyr)
  library(ggplot2)
  library(circlize)
  library(logistf)
  library(igraph)
  library(data.table)
  library(ComplexHeatmap)
  library(limma)
  library(SummarizedExperiment)
  library(S4Vectors)
})
set.seed(2025)

## ----data---------------------------------------------------------------------
data(monkeyHSPC, package = "barbieQ")

## ----example, fig.width=8, fig.height=6---------------------------------------
## Passing `object`, `sampleMetadata` and `factorColors` for optional
exampleBBQ <- createBarbieQ(
  object = SummarizedExperiment::assay(monkeyHSPC), 
  sampleMetadata = SummarizedExperiment::colData(monkeyHSPC)$sampleMetadata
)

## ----update metadata----------------------------------------------------------
updateSampleMetadata <- exampleBBQ$sampleMetadata %>% 
  as.data.frame() %>% 
  select(Celltype, Months) %>% 
  mutate(Phase = ifelse(Months < 6, "early", ifelse(Months >=55, "late", "mid"))) %>% 
  mutate(Celltype = gsub("(Gr).*", "\\1", Celltype))

SummarizedExperiment::colData(exampleBBQ)$sampleMetadata <- S4Vectors::DataFrame(updateSampleMetadata)

exampleBBQ$sampleMetadata

## ----subset samples-----------------------------------------------------------
flag_sample <- exampleBBQ$sampleMetadata$Phase == "mid"
exampleBBQ <- exampleBBQ[, flag_sample]
exampleBBQ$sampleMetadata

## ----tag top------------------------------------------------------------------
## Check out minimum group size.
table(exampleBBQ$sampleMetadata$Celltype)
## Tag top Barcodes.
exampleBBQ <- tagTopBarcodes(barbieQ = exampleBBQ, nSampleThreshold = 6)

## ----tag top v1, fig.width=8, fig.height=6------------------------------------
## visualize contribution of top vs. bottom barcodes
plotBarcodePareto(barbieQ = exampleBBQ) |> plot()

## ----tag top v2, fig.width=5, fig.height=3------------------------------------
## visualize collective contribution of top vs. bottom barcodes
plotBarcodeSankey(barbieQ = exampleBBQ) |> plot()

## ----subset top barcodes------------------------------------------------------
flag_barcode <- SummarizedExperiment::rowData(exampleBBQ)$isTopBarcode$isTop
exampleBBQ <- exampleBBQ[flag_barcode,]

## ----sample pair cor, fig.width=6, fig.height=4-------------------------------
## visualize sample pair wise correlation
plotSamplePairCorrelation(barbieQ = exampleBBQ) |> plot()

## ----bartools, fig.width=20, fig.height=10, eval=FALSE------------------------
# # devtools::install_github("DaneVass/bartools", dependencies = TRUE, force = TRUE)
# #
# # dge <- DGEList(
# #   counts = assay(exampleBBQ),
# #   group = exampleBBQ$sampleMetadata$Celltype)
# #
# # bartools::plotBarcodeHistogram(dge)

## ----speckle, eval=FALSE------------------------------------------------------
# ## to inspect variance
# # speckle::plotCellTypeMeanVar(assay(exampleBBQ))
# # speckle::plotCellTypePropsMeanVar(assay(exampleBBQ))

## ----test diffProp, fig.width=6, fig.height=5---------------------------------
## test Barcode differential proportion between sample groups
## Defult transformation: asin-sqrt
asinTrans <- testBarcodeSignif(
  barbieQ = exampleBBQ,
  contrastFormula = "(CelltypeNK_CD56n_CD16p) - (CelltypeB+CelltypeGr+CelltypeT+CelltypeNK_CD56p_CD16n)/4",
  method = "diffProp", transformation = "asin-sqrt"
)

## Alternatively: using logit transformation
logitTrans <- testBarcodeSignif(
  barbieQ = exampleBBQ,
  contrastFormula = "(CelltypeNK_CD56n_CD16p) - (CelltypeB+CelltypeGr+CelltypeT+CelltypeNK_CD56p_CD16n)/4",
  method = "diffProp", transformation = "logit"
)

## Alternatively: no transformation
noTrans <- testBarcodeSignif(
  barbieQ = exampleBBQ,
  contrastFormula = "(CelltypeNK_CD56n_CD16p) - (CelltypeB+CelltypeGr+CelltypeT+CelltypeNK_CD56p_CD16n)/4",
  method = "diffProp", transformation = "none"
)

## ----vis MA plot, fig.width=8, fig.height=4-----------------------------------
(plotBarcodeMA(asinTrans) + coord_trans(x = "log10"))|> plot()
plotBarcodeMA(logitTrans)|> plot()
(plotBarcodeMA(noTrans) + coord_trans(x = "log10"))|> plot()


## ----vis testing HP, fig.width=13, fig.height=5, eval=FALSE-------------------
# plotSignifBarcodeHeatmap(asinTrans) |> plot()
# plotSignifBarcodeHeatmap(logitTrans) |> plot()

## ----vis testing prop, fig.width=12, fig.height=10----------------------------
plotSignifBarcodeProportion(asinTrans) |> plot()
plotSignifBarcodeProportion(logitTrans) |> plot()

## ----test diffOcc, fig.width=6, fig.height=5----------------------------------
## test Barcode differential occurrence between sample groups
## set up the targets (sample conditions)
targets <- exampleBBQ$sampleMetadata %>% 
  as.data.frame() %>% 
  mutate(Group = ifelse(
    Celltype == "NK_CD56n_CD16p",
    "NK_CD56n_CD16p",
    "B.Gr.T.NK_CD56p_CD16n"))

exampleBBQ <- testBarcodeSignif(
  barbieQ = exampleBBQ, 
  sampleMetadata = targets[,"Group", drop=FALSE], 
  method = "diffOcc"
)

## ----vis MA plot occ, fig.width=8, fig.height=4-------------------------------
plotBarcodeMA(exampleBBQ) |> plot()

## ----vis testing HP occ, fig.width=13, fig.height=5, eval=FALSE---------------
# plotSignifBarcodeHeatmap(exampleBBQ) |> plot()

## ----vis testing prop occ, fig.width=12, fig.height=10------------------------
plotSignifBarcodeProportion(exampleBBQ) |> plot()

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