---
title: "barbieQ: An R package for analysing barcode count data from clonal tracking experiments"
output: 
  BiocStyle::html_document:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Quick start to barbieQ}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
date: "`r Sys.Date()`"
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# Introduction

The barbieQ package provides a series of robust statistical tools for analysing barcode count data generated from cell *clonal tracking* (*lineage tracing*) experiments. 

In these experiments, an initial cell and its offspring collectively form a *clone* (or *lineage*). A unique DNA barcode, incorporated into the genome of an initial cell, is inherited by all its progeny within the clone. This one-to-one mapping of barcodes to clones enables  tracking of clonal behaviours. By quantifying barcode counts, researchers can measure the abundance of individual clones under various experimental conditions or perturbations. 

While existing tools for barcode count data analysis primarily rely on qualitative interpretation through visualizations, they often lack robust methods to model the sources of barcode variability. ^[[barcodetrackR](https://bioconductor.org/packages/barcodetrackR/), [CellDestiny](https://github.com/TeamPerie/CellDestiny), [genBaRcode](https://CRAN.R-project.org/package=genBaRcode)]^

To address this gap, this R software package, **barbieQ**, provides advanced statistical methods to model barcode variability. The package supports preprocessing, visualization, and statistical testing to identify barcodes with significant differences in proportions or occurrences across experimental conditions. Key functionalities include initializing data structures, filtering barcodes, and applying regression models to test for significant clonal changes.

The main functions include:

 * `createBarbieQ()`
 
 * `tagTopBarcodes()`
 
 * `plotBarcodePairCorrelation()`
 
 * `clusterCorrelatingBarcodes()`
 
 * `plotSamplePairCorrelation()`
 
 * `plotBarcodeProportion()`
 
 * `testBarcodeSignif()`
 
 * `plotSignifBarcodeProportion()`

# Intallation

```{r install, eval=FALSE}
# You can install the stable 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")
```

# Load Dependecy

```{r 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)
```

# Load Package Data
 
- Monkey HSPC Cell Barcoding Data (monkeyHSPC)

A subset of data from a study on monkey HSPC cell expansion using barcoding technique.^[NK clonal expansion](http://dx.doi.org/10.1126/sciimmunol.aat9781)),  [barcodetrackRData](https://github.com/dunbarlabNIH/barcodetrackRData)^
Barcode counts within different samples were used to interpret the patterns of HSPC differentiation.

It is a `SummarizedExperiment` object created using function `barbieQ::createBarbie`, containing a barcode count matrix with 16,603 rows and 62 columns, and a data frame of sample metadata.

```{r data}
data(monkeyHSPC, package = "barbieQ")
```

# Example

## Create `barbieQ` Object

Please start with creating a `barbieQ` structure by taking the barcode count matrix as input, using function `createBarbieQ()`.

By creating a `barbieQ` object, you've automatically processed a series of data transformation and saved transformed data in the `barbieQ` object, for the convenience of the following analyses.

```{r 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
)
```

## Subset Dataset

Here we subset the object by selecting samples of certain stage of collection time.

Define "early", "mid", and "late" stages based on `Months`, and clean up `Celltype`.

```{r 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 dataset and keep samples of "mid" stage only.

```{r subset samples}
flag_sample <- exampleBBQ$sampleMetadata$Phase == "mid"
exampleBBQ <- exampleBBQ[, flag_sample]
exampleBBQ$sampleMetadata
```


## Tag Top Contributing Barcodes

A filtering step is recommended to filter out barcodes that consistently show low counts across the dataset. The retained barcodes are considered with essential contribution, called *"top barcodes"*. 

 - By applying the function `tagTopBarcodes()` to the `barbieQ` object, you're determining which barcodes are *"top barcodes"* while tagging them in the object.
 
 - In this example dataset, we are interested in the differences on barcode outcomes between cell types, so we will group samples by cell types. We set up the `nSampleThreshold` as the minimum group size, which is `6`.

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

 - Once *"top barcodes"* are determined and tagged, it's good to check their contributions before cutting off the *"bottom barcodes"* that are considered as non-essential contributors. 
 
     - By applying the function `plotBarcodePareto()` to the `barbieQ` object, you are visualizing the contribution of each barcode colour coded by *"top"* or *"bottom"*. (Here, the contribution refers to the average proportion of individual barcodes across samples in the dataset.)
     
     - By applying the function `plotBarcodeSankey()` to the `barbieQ` object, you are visualizing the collective contribution of the *"top"* or *"bottom"* barcode groups. 

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

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

 - Once you are happy with the classification of *"top"* and *"bottom"* barcodes, you can filter out the *"bottom"* barcodes by subsetting the `SE` object based on the tagged array.
 
```{r subset top barcodes}
flag_barcode <- SummarizedExperiment::rowData(exampleBBQ)$isTopBarcode$isTop
exampleBBQ <- exampleBBQ[flag_barcode,]
```

## Visualize Sample Correlation

To generally understand the sample similarity, you can visualize sample pair-wise correlations in a check board style by applying the function `plotSamplePairCorrelation()` to the `barbieQ()` object.

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

## Test Barcodes and Identify Significant Changes

Based on the understanding of sample conditions that likely to be the source of variability in barcode outcomes, you can robustly test the significance of the barcode changes between the sample conditions, by applying the function `testBarcodeSignif()` to the `barbieQ` object. The testing results will be saved in the object, and can be further visualized using functions: `plotBarcodeMA()`, `plotSignifBarcodeHeatmap()`, `plotSignifBarcodeProportion()`, and etc.
 
 - By specifying the parameter `method` by `diffProp` (defaulted), you are testing individual barcodes' **differential proportion** between conditions.

 - By specifying the parameter `method` by `diffOcc`, you are testing individual barcodes' **differential occurrence** between conditions.

```{r test diffProp, fig.width=6, fig.height=5}
## test Barcode differential proportion between sample groups
exampleBBQ <- testBarcodeSignif(
  barbieQ = exampleBBQ,
  contrastFormula = "(CelltypeNK_CD56n_CD16p) - (CelltypeB+CelltypeGr+CelltypeT+CelltypeNK_CD56p_CD16n)/4",
  method = "diffProp"
)
```


```{r vis MA plot, fig.width=8, fig.height=4}
plotBarcodeMA(exampleBBQ) |> plot()
```

```{r vis testing HP, fig.width=14, fig.height=6}
plotSignifBarcodeHeatmap(exampleBBQ) |> plot()
```

```{r vis testing prop, fig.width=8, fig.height=5}
plotSignifBarcodeProportion(exampleBBQ) |> plot()
```


# References

We are currently writing a paper to introduce the methods and approaches implemented in this `barbieQ` package and will update with a citation once available.

# SessionInfo

```{r session info}
sessionInfo()
```