## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    crop = NULL,
    fig.width = 6,
    fig.height = 4.5,
    dpi = 96,
    fig.retina = 1
)

## ----install, eval = FALSE----------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("scProkaR")

## ----load-packages------------------------------------------------------------
library(scProkaR)
library(SingleCellExperiment)
library(ggplot2)

## ----load-data----------------------------------------------------------------
data("sce1", package = "scProkaR")
sce1

## ----describe-design----------------------------------------------------------
table(sce1$treatments, sce1$timepoints)

## ----describe-coldata---------------------------------------------------------
head(as.data.frame(colData(sce1)), 4)

## ----describe-features--------------------------------------------------------
head(rownames(sce1), 8)

## ----describe-depth-----------------------------------------------------------
summary(Matrix::colSums(assay(sce1, "counts")))

## ----subset-cells-------------------------------------------------------------
set.seed(101)
keep_cells <- unlist(lapply(
    split(seq_len(ncol(sce1)), sce1$sample),
    function(idx) sample(idx, 250)
))
sce_small <- sce1[, sort(keep_cells)]
dim(sce_small)

## ----create-object------------------------------------------------------------
counts_matrix <- assay(sce_small, "counts")
cell_meta <- as.data.frame(colData(sce_small))

set.seed(101)
sce <- CreateBacObject(
    counts_matrix,
    cell_metadata = cell_meta,
    sample_col = "sample",
    batch_col = "treatments",
    condition_col = "treatments",
    time_col = "timepoints",
    organism = "bacteria",
    run_unintegrated = TRUE,
    unintegrated_dims = seq_len(20),
    unintegrated_k = 20,
    unintegrated_resolution = 1
)
sce

## ----inspect-object-----------------------------------------------------------
reducedDimNames(sce)
names(S4Vectors::metadata(sce)$scProkaR)

## ----tenx-entry-point---------------------------------------------------------
tenx_call <- quote(CreateBacObject(
    "path/to/filtered_feature_bc_matrix",
    sample_id_value = "PMB2_7",
    sample_col = "sample"
))
tenx_call

## ----seurat-entry-point-------------------------------------------------------
if (requireNamespace("SeuratObject", quietly = TRUE)) {
    demo_cells <- seq_len(200)
    seu <- SeuratObject::CreateSeuratObject(
        counts = counts_matrix[, demo_cells],
        meta.data = cell_meta[demo_cells, , drop = FALSE]
    )
    sce_from_seurat <- CreateBacObject(
        seu,
        seurat_assay = "RNA",
        seurat_layer = "counts",
        sample_col = "sample",
        batch_col = "treatments"
    )
    print(dim(sce_from_seurat))
    print(head(colnames(colData(sce_from_seurat))))
} else {
    message("SeuratObject is not installed; skipping the Seurat entry point.")
}

## ----run-qc-------------------------------------------------------------------
sce <- RunBacQC(sce)
head(as.data.frame(colData(sce))[, c(
    "total_counts",
    "detected_features",
    "rrna_fraction",
    "ribo_fraction"
)], 4)

## ----check-rrna---------------------------------------------------------------
sum(grepl("^(rrs|rrl|rrf)", rownames(sce), ignore.case = TRUE))
summary(sce$rrna_fraction)

## ----summarise-ribo-----------------------------------------------------------
round(tapply(sce$ribo_fraction, sce$treatments, median), 4)

## ----filter-cells-------------------------------------------------------------
sce <- FilterBacCells(
    sce,
    min_counts = 50,
    min_features = 30,
    max_rrna_fraction = 0.2
)
dim(sce)
table(sce$treatments)

## ----integrate-mnn------------------------------------------------------------
if (requireNamespace("batchelor", quietly = TRUE)) {
    set.seed(101)
    sce <- IntegrateBacData(
        sce,
        batch_col = "treatments",
        method = "mnn",
        dims = seq_len(20)
    )
    print(reducedDimNames(sce))
} else {
    message("batchelor is not installed; skipping MNN integration.")
}

## ----integrate-harmony--------------------------------------------------------
if (requireNamespace("harmony", quietly = TRUE)) {
    set.seed(101)
    sce <- IntegrateBacData(
        sce,
        batch_col = "treatments",
        method = "harmony",
        dims = seq_len(20)
    )
    print(reducedDimNames(sce))
} else {
    message("harmony is not installed; skipping Harmony integration.")
}

## ----register-embedding-------------------------------------------------------
pca_emb <- reducedDim(sce, "PCA")[, seq_len(20), drop = FALSE]
batch_vec <- as.character(sce$treatments)
batch_centres <- t(vapply(
    split(seq_len(nrow(pca_emb)), batch_vec),
    function(idx) colMeans(pca_emb[idx, , drop = FALSE]),
    numeric(ncol(pca_emb))
))
centred_emb <- pca_emb - batch_centres[batch_vec, , drop = FALSE]

sce <- RegisterIntegrationEmbedding(
    sce,
    embedding = centred_emb,
    method_name = "centred"
)
reducedDimNames(sce)

## ----choose-reduction---------------------------------------------------------
integrated_names <- grep(
    "^integrated_",
    reducedDimNames(sce),
    value = TRUE
)
main_reduction <- if ("integrated_mnn" %in% integrated_names) {
    "integrated_mnn"
} else {
    "integrated_centred"
}
main_reduction

## ----run-umap-----------------------------------------------------------------
if (requireNamespace("uwot", quietly = TRUE)) {
    set.seed(101)
    sce <- RunIntegratedUMAP(
        sce,
        reduction = main_reduction,
        n_neighbors = 20
    )
    vis_reduction <- paste0("umap_", main_reduction)
} else {
    message("uwot is not installed; plotting the embedding directly.")
    vis_reduction <- main_reduction
}
vis_reduction

## ----choose-unintegrated------------------------------------------------------
unintegrated_reduction <- if (
    "umap.unintegrated" %in% reducedDimNames(sce)
) {
    "umap.unintegrated"
} else {
    "PCA"
}
unintegrated_reduction

## ----run-clustering-----------------------------------------------------------
set.seed(101)
sce <- RunIntegratedClustering(
    sce,
    reduction = main_reduction,
    cluster_col = "integrated_clusters",
    k = 20,
    algorithm = "louvain",
    resolution = 1
)
table(sce$integrated_clusters)

## ----compare-clusters---------------------------------------------------------
table(sce$integrated_clusters, sce$clusters)

## ----plot-batch---------------------------------------------------------------
PlotReduction(
    sce,
    reduction = vis_reduction,
    colour_by = "treatments",
    point_size = 0.5
)

## ----plot-labels--------------------------------------------------------------
PlotReduction(
    sce,
    reduction = vis_reduction,
    colour_by = "clusters",
    point_size = 0.5
)

## ----plot-time----------------------------------------------------------------
sce$timepoint_factor <- factor(
    sce$timepoints,
    levels = c("0", "1", "4", "7")
)
PlotReduction(
    sce,
    reduction = vis_reduction,
    colour_by = "timepoint_factor",
    point_size = 0.5
)

## ----plot-facet---------------------------------------------------------------
PlotReduction(
    sce,
    reduction = vis_reduction,
    colour_by = "integrated_clusters",
    facet_by = "treatments",
    point_size = 0.5
)

## ----integration-overview-----------------------------------------------------
overview <- PlotIntegrationOverview(
    sce,
    unintegrated_reduction = unintegrated_reduction,
    integrated_reduction = vis_reduction,
    batch_col = "treatments",
    label_col = "clusters",
    split_by = "treatments",
    point_size = 0.5
)
names(overview)

## ----overview-before-batch----------------------------------------------------
overview$unintegrated_by_batch

## ----overview-after-batch-----------------------------------------------------
overview$integrated_by_batch

## ----overview-after-label-----------------------------------------------------
overview$integrated_by_label

## ----run-benchmark------------------------------------------------------------
bench_methods <- c("PCA", integrated_names)
set.seed(101)
bench <- BenchmarkIntegration(
    sce,
    batch_col = "treatments",
    label_col = "clusters",
    methods = bench_methods
)
bench$ranking

## ----benchmark-scores---------------------------------------------------------
head(bench$scores, 8)

## ----benchmark-overall--------------------------------------------------------
bench$plots$overall

## ----benchmark-tradeoff-------------------------------------------------------
bench$plots$tradeoff

## ----benchmark-key-metrics----------------------------------------------------
bench$plots$key_metrics

## ----assign-replicates--------------------------------------------------------
set.seed(101)
replicate_id <- integer(ncol(sce))
for (idx in split(seq_len(ncol(sce)), sce$sample)) {
    replicate_id[idx] <- sample(rep(seq_len(3), length.out = length(idx)))
}
sce$replicate <- factor(replicate_id)
sce$pb_id <- paste(sce$sample, sce$replicate, sep = "_r")
table(sce$sample, sce$replicate)

## ----aggregate-pseudobulk-----------------------------------------------------
pb <- aggregate_pseudobulk(
    sce,
    sample_cols = c(
        "pb_id", "sample", "replicate", "treatments", "timepoints"
    ),
    aggregation = "sum",
    min_cells = 10
)
dim(pb)
head(as.data.frame(colData(pb))[, c(
    "pb_id", "treatments", "timepoints", "ncells"
)], 4)

## ----normalize-pseudobulk-----------------------------------------------------
pb <- normalize_pseudobulk(pb, method = "TMM")
summary(pb$norm.factors)

## ----pseudobulk-pca-----------------------------------------------------------
pb_pca <- stats::prcomp(
    t(as.matrix(assay(pb, "logcounts"))),
    scale. = TRUE
)
reducedDim(pb, "pb_pca") <- pb_pca$x[, seq_len(4), drop = FALSE]
pb$timepoint_factor <- factor(pb$timepoints, levels = c("0", "1", "4", "7"))
PlotReduction(pb, reduction = "pb_pca", colour_by = "treatments",
    point_size = 3)

## ----pseudobulk-pca-time------------------------------------------------------
PlotReduction(pb, reduction = "pb_pca", colour_by = "timepoint_factor",
    point_size = 3)

## ----pairwise-de--------------------------------------------------------------
pairwise_time <- run_edger_pairwise_de(
    pb,
    group_col = "timepoints",
    covariates = "treatments",
    contrast_type = "reference",
    reference_level = "0",
    normalization = "TMM",
    fdr_cutoff = 0.05
)
names(pairwise_time$tables)

## ----pairwise-summary---------------------------------------------------------
table(
    pairwise_time$combined_table$contrast,
    pairwise_time$combined_table$significant
)

## ----pairwise-top-genes-------------------------------------------------------
tab_4_vs_0 <- pairwise_time$tables[["4_vs_0"]]
head(tab_4_vs_0[order(tab_4_vs_0$FDR), c("gene", "logFC", "FDR")], 10)

## ----pairwise-volcano---------------------------------------------------------
plot_pairwise_de_volcano(
    pairwise_time,
    contrast = "4_vs_0",
    fdr_cutoff = 0.05,
    lfc_cutoff = 1,
    top_n_labels = 8
)

## ----pairwise-ma--------------------------------------------------------------
plot_pairwise_de_ma(
    pairwise_time,
    contrast = "4_vs_0",
    fdr_cutoff = 0.05
)

## ----pairwise-composite-------------------------------------------------------
pairwise_dose <- run_edger_pairwise_de(
    pb,
    group_col = "treatments",
    combine_group_cols = c("treatments", "timepoints"),
    contrast_type = "manual",
    manual_pairs = data.frame(
        group1 = c("PMB2__7", "PMB0.5__7"),
        group2 = c("control__7", "control__7"),
        stringsAsFactors = FALSE
    ),
    normalization = "TMM"
)
names(pairwise_dose$tables)

## ----pairwise-composite-volcano-----------------------------------------------
plot_pairwise_de_volcano(
    pairwise_dose,
    contrast = "PMB2__7_vs_control__7",
    fdr_cutoff = 0.05,
    lfc_cutoff = 1,
    top_n_labels = 8
)

## ----spline-de----------------------------------------------------------------
spline_de <- run_edger_spline_de(
    pb,
    time_col = "timepoints",
    condition_col = "treatments",
    df = 2,
    normalization = "TMM",
    return_curve_for = "significant",
    fdr_cutoff = 0.05
)
names(spline_de$by_condition)

## ----spline-summary-----------------------------------------------------------
table(
    spline_de$combined_table$treatments,
    spline_de$combined_table$significant
)

## ----spline-top-genes---------------------------------------------------------
spline_top <- spline_de$combined_table[
    order(spline_de$combined_table$FDR),
    c("gene", "treatments", "FDR")
]
head(spline_top, 10)

## ----spline-curves------------------------------------------------------------
curve_genes <- head(unique(spline_top$gene), 6)
plot_time_series_deg_curves(spline_de, genes = curve_genes, ncol = 3)

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