## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    echo = TRUE,
    message = FALSE,
    warning = FALSE,
    fig.cap = "",
    tidy = TRUE
)
options(timeout = 600)

## ----setup--------------------------------------------------------------------
set.seed(42)
library(dominoSignal)
library(dplyr)
library(kableExtra)

# BiocFileCache helps with file management across sessions
bfc <- BiocFileCache::BiocFileCache(ask = FALSE)
data_url <- "https://zenodo.org/records/10951634/files/pbmc_domino_built.rds"
tmp_path <- BiocFileCache::bfcrpath(bfc, data_url)
dom <- readRDS(tmp_path)

## ----extract inputs-----------------------------------------------------------
rl_map <- dom_database(dom, name_only = FALSE)
features <- dom_tf_activation(dom)
# Access counts slot directly to maintain sparse matrix format
counts <- slot(dom, "counts")
zscores <- dom_zscores(dom)
tf_targets <- dom_linkages(dom, link_type = "tf-target")
clusters <- dom_clusters(dom, labels = TRUE)

## ----create replicates, results = "hide"--------------------------------------
build_reps <- function(clusters, seed, keep_frac, shuffle_frac) {
    # Local seed for reproducibility without affecting global seed
    withr::with_seed(seed, {
        n_cells <- length(clusters)
        keep_cells <- sort(sample.int(n_cells, size = round(n_cells * keep_frac), replace = FALSE))

        resampled_clusters <- clusters[keep_cells]
        features <- features[ , keep_cells, drop = FALSE]
        counts <- counts[ , keep_cells, drop = FALSE]
        zscores <- zscores[ , keep_cells, drop = FALSE]
        shuffled_idx <- sample(seq_along(resampled_clusters),
            size = round(shuffle_frac * length(resampled_clusters)), replace = FALSE
        )

        resampled_clusters[shuffled_idx] <- sample(resampled_clusters[shuffled_idx])

        dom_shuffle <- create_domino(
            rl_map = rl_map,
            features = features,
            counts = counts,
            z_scores = zscores,
            clusters = resampled_clusters,
            tf_targets = tf_targets,
            use_complexes = TRUE,
            remove_rec_dropout = FALSE,
            verbose = FALSE
        )

        dom_shuffle <- build_domino(
            dom = dom_shuffle,
            max_tf_pval = 0.05,
            max_tf_per_clust = 25,
            max_rec_per_tf = 25,
            rec_tf_cor_threshold = 0.1,
            min_rec_percentage = 0.05
        )
    })

    return(dom_shuffle)
}

params <- data.frame(
    set = c(161803, 314159, 6631034, 6021023, 831446, 98067),
    keep = c(0.95, 0.95, 0.95, 0.95, 0.9, 0.9),
    shuffle = c(0.15, 0.1, 0.15, 0.1, 0.15, 0.1)
)

# This does take a few minutes to run
dom_list <- purrr::pmap(params, function(set, keep, shuffle) {
    build_reps(clusters, seed = set, keep_frac = keep, shuffle_frac = shuffle)
})

names(dom_list) <- paste0("P", seq_along(dom_list))

## ----summarize linkages, message = TRUE---------------------------------------
subject_metas <- data.frame(
    patient = c("P1", "P2", "P3", "P4", "P5", "P6"),
    time = c("pre", "post", "pre", "post", "pre", "post"),
    condition = c("control", "treatment", "treatment", "control", "control", "control"),
    batch = c("A", "A", "B", "B", "C", "C"),
    stringsAsFactors = FALSE
)

link_summary <- summarize_linkages(dom_list, subject_metas,
    subject_names = subject_metas$patient
)

print(link_summary)

## ----subset treatment, message = TRUE-----------------------------------------
link_treat <- subset(link_summary, subset = condition == "treatment")
print(link_treat)

## ----subset examples, message = TRUE------------------------------------------
subset(link_summary, subset = subject_names %in% c("P1", "P3"))
subset(link_summary, subset = condition == "treatment" | time == "pre")

## ----all linkages-------------------------------------------------------------
count_all <- count_linkage(link_summary,
    cluster = "CD8_T_cell",
    linkage = "rec_lig"
) %>% slice_sample(n = 5)

## ----eval = FALSE-------------------------------------------------------------
# kableExtra::kable_styling(
#     knitr::kable(count_all,
#         caption = "Total counts of receptor-ligand linkages in CD8 T Cells"
#     ),
#     bootstrap_options = c("striped", "hover")
# )

## ----echo = FALSE-------------------------------------------------------------
kableExtra::kable_styling(knitr::kable(count_all,
        caption = "Total counts of receptor-ligand linkages in CD8 T Cells"),
    bootstrap_options = c("striped", "hover"))

## ----subset subjects----------------------------------------------------------
count_subjects <- count_linkage(link_summary,
    cluster = "CD8_T_cell",
    linkage = "tfs", subject_names = c("P1", "P2")
)

## ----eval = FALSE-------------------------------------------------------------
# kableExtra::kable_styling(
#     knitr::kable(count_subjects[1:5, ],
#         caption = "Subset of counts of TF linkages in patients P1 & P2 in CD8 T Cells"
#     ),
#     bootstrap_options = c("striped", "hover")
# )

## ----echo = FALSE-------------------------------------------------------------
kableExtra::kable_styling(knitr::kable(count_subjects[1:5, ],
        caption = "Subset of counts of TF linkages in patients P1 & P2 in CD8 T Cells"),
    bootstrap_options = c("striped", "hover"))

## ----time treatment counts----------------------------------------------------
counts_time_trt <- count_linkage(link_treat,
    cluster = "CD14_monocyte",
    linkage = "tfs_rec", group.by = "time"
)

## ----eval = FALSE-------------------------------------------------------------
# kableExtra::kable_styling(
#     knitr::kable(slice_sample(counts_time_trt, n = 5),
#         caption = "Counts of TF-Receptor linkages in CD14 Monocytes grouped by time for treatment subjects"
#     ),
#     bootstrap_options = c("striped", "hover")
# )

## ----echo = FALSE-------------------------------------------------------------
kableExtra::kable_styling(
    knitr::kable(slice_sample(counts_time_trt, n = 5),
        caption = "Counts of TF-Receptor linkages in CD14 Monocytes grouped by time for treatment subjects"
    ),
    bootstrap_options = c("striped", "hover")
)

## ----time subject counts------------------------------------------------------
counts_time_subj <- count_linkage(link_summary,
    cluster = "CD14_monocyte",
    linkage = "tfs_rec", group.by = "time", subject_names = c("P2", "P3")
)
identical(counts_time_trt, counts_time_subj)

## ----condition counts---------------------------------------------------------
counts_condition <- count_linkage(link_summary,
    cluster = "naive_CD4_T_cell",
    linkage = "incoming_lig", group.by = "condition"
) %>% slice_sample(n = 5)

## ----eval = FALSE-------------------------------------------------------------
# kableExtra::kable_styling(
#     knitr::kable(counts_condition,
#         caption = "Counts of incoming ligand linkages in naive CD4 T Cells grouped by condition"
#     ),
#     bootstrap_options = c("striped", "hover")
# )

## ----echo = FALSE-------------------------------------------------------------
kableExtra::kable_styling(knitr::kable(counts_condition,
        caption = "Counts of incoming ligand linkages in naive CD4 T Cells grouped by condition"),
    bootstrap_options = c("striped", "hover"))

## ----batch counts-------------------------------------------------------------
counts_batch <- count_linkage(link_summary,
    cluster = "CD14_monocyte",
    linkage = "rec", group.by = "batch"
) %>% slice_sample(n = 5)

## ----eval = FALSE-------------------------------------------------------------
# kableExtra::kable_styling(
#     knitr::kable(counts_batch,
#         caption = "Counts of receptor linkages in CD14 Monocytes grouped by batch"
#     ),
#     bootstrap_options = c("striped", "hover")
# )

## ----echo = FALSE-------------------------------------------------------------
kableExtra::kable_styling(knitr::kable(counts_batch,
        caption = "Counts of receptor linkages in CD14 Monocytes grouped by batch"),
    bootstrap_options = c("striped", "hover"))

## ----test time treat----------------------------------------------------------
dcst_time_trt <- test_differential_linkages(link_treat,
    cluster = "CD14_monocyte", linkage = "tfs_rec", group.by = "time"
)

## ----eval = FALSE-------------------------------------------------------------
# kableExtra::kable_styling(
#     knitr::kable(slice_sample(dcst_time_trt, n = 5),
#         caption = "DCST of TF-Receptor linkages in CD14 Monocytes grouped by time for treatment subjects"
#     ),
#     bootstrap_options = c("striped", "hover")
# )

## ----echo = FALSE-------------------------------------------------------------
kableExtra::kable_styling(knitr::kable(slice_sample(dcst_time_trt, n = 5),
        caption = "DCST of TF-Receptor linkages in CD14 Monocytes grouped by time for treatment subjects"),
    bootstrap_options = c("striped", "hover"))

## ----test time subject--------------------------------------------------------
dcst_time_subj <- test_differential_linkages(link_summary,
    cluster = "CD14_monocyte",
    linkage = "tfs_rec", group.by = "time", subject_names = c("P2", "P3")
)
identical(dcst_time_trt, dcst_time_subj)

## ----no subjects, error = TRUE------------------------------------------------
try({
test_differential_linkages(link_treat,
    cluster = "CD14_monocyte",
    linkage = "tfs_rec", group.by = "condition"
)
})

## ----dcst time----------------------------------------------------------------
dcst_time <- test_differential_linkages(link_summary,
    cluster = "CD8_T_cell",
    linkage = "rec", group.by = "time"
) %>% dplyr::arrange(desc(odds.ratio))

## ----eval = FALSE-------------------------------------------------------------
# kableExtra::kable_styling(
#     knitr::kable(slice_sample(dcst_time, n = 5),
#         caption = "DCST of receptor-ligand linkages in CD8 T Cells grouped by time"
#     ),
#     bootstrap_options = c("striped", "hover")
# )

## ----echo = FALSE-------------------------------------------------------------
kableExtra::kable_styling(knitr::kable(slice_sample(dcst_time, n = 5),
        caption = "DCST of receptor-ligand linkages in CD8 T Cells grouped by time"),
    bootstrap_options = c("striped", "hover"))

## ----dcst batch---------------------------------------------------------------
dcst_batch <- test_differential_linkages(link_summary,
    cluster = "naive_CD4_T_cell",
    linkage = "rec", group.by = "batch"
) %>% arrange(p.value)

## ----eval = FALSE-------------------------------------------------------------
# kableExtra::kable_styling(
#     knitr::kable(slice_sample(dcst_batch, n = 5),
#         caption = "DCST of receptors in naive CD4 T Cells grouped by batch"
#     ),
#     bootstrap_options = c("striped", "hover")
# )

## ----echo = FALSE-------------------------------------------------------------
kableExtra::kable_styling(knitr::kable(slice_sample(dcst_batch, n = 5),
        caption = "DCST of receptors in naive CD4 T Cells grouped by batch"),
    bootstrap_options = c("striped", "hover"))

## ----dcst patient-------------------------------------------------------------
dcst_patient <- test_differential_linkages(link_summary,
    cluster = "naive_CD4_T_cell",
    linkage = "rec_lig", group.by = "patient"
) %>% arrange(p.value)

## ----eval = FALSE-------------------------------------------------------------
# kableExtra::kable_styling(
#     knitr::kable(slice_sample(dcst_patient, n = 5),
#         caption = "DCST of receptor-ligand linkages in naive CD4 T Cells grouped by patient"
#     ),
#     bootstrap_options = c("striped", "hover")
# )

## ----echo = FALSE-------------------------------------------------------------
kableExtra::kable_styling(knitr::kable(slice_sample(dcst_patient, n = 5),
        caption = "DCST of receptor-ligand linkages in naive CD4 T Cells grouped by patient"),
    bootstrap_options = c("striped", "hover"))

## ----plot time treat----------------------------------------------------------
plot_differential_linkages(
    differential_linkages = slice_sample(dcst_time_trt, n = 15),
    test_statistic = "p.value", stat_range = c(0, 1), stat_ranking = "ascending"
)

## ----plot time----------------------------------------------------------------
plot_differential_linkages(
    differential_linkages = slice_sample(dcst_time, n = 15),
    test_statistic = "p.value", gradient = c(minimum = "#c44f4f", maximum = "#ffd771"),
    stat_ranking = "ascending"
)

## ----plot dcst time odds------------------------------------------------------
plot_differential_linkages(
    differential_linkages = slice_sample(dcst_time, n = 15),
    test_statistic = "odds.ratio", stat_range = c(0, 4), stat_ranking = "descending",
    gradient = c(minimum = "gray90", maximum = "orange")
)

## ----plot dcst batch----------------------------------------------------------
plot_differential_linkages(
    differential_linkages = slice_sample(dcst_batch, n = 15),
    test_statistic = "p.value", stat_range = c(0, 1), stat_ranking = "ascending",
    gradient = c(minimum = "#ff4c43", maximum = "#ffe6c0")
)

## ----plot dcst group colors---------------------------------------------------
plot_differential_linkages(
    differential_linkages = slice_sample(dcst_batch, n = 15),
    test_statistic = "p.value", stat_range = c(0, 1), stat_ranking = "ascending",
    gradient = c(minimum = "#7b43ff", maximum = "#c0c5ff"),
    group_palette = c(A = "#64a5e6", B = "#71fcad", C = "#c968c3")
)

## ----example workflow, eval = FALSE-------------------------------------------
# subject_metas <- data.frame(
#     patient = c("P1", "P2", "P3", "P4", "P5", "P6"),
#     time = c("pre", "post", "pre", "post", "pre", "post"),
#     condition = c("control", "treatment", "treatment", "control", "control", "control"),
#     batch = c("A", "A", "B", "B", "C", "C"),
#     stringsAsFactors = FALSE
# )
# 
# link_summary <- summarize_linkages(dom_list, subject_metas,
#     subject_names = subject_metas$patient
# )
# 
# counts_batch <- count_linkage(link_summary,
#     cluster = "CD14_monocyte",
#     linkage = "incoming_lig", group.by = "batch"
# )
# 
# dcst_batch <- test_differential_linkages(link_summary,
#     cluster = "CD14_monocyte",
#     linkage = "incoming_lig", group.by = "batch"
# )
# 
# plot_differential_linkages(dcst_batch,
#     test_statistic = "p.value", stat_ranking = "ascending",
#     gradient = c("minimum" = "#e02c2c", "maximum" = "#feffcd"),
#     group_palette = c(A = "#2a6aaa", B = "#07bd56", C = "#a244a3")
# )

## -----------------------------------------------------------------------------
Sys.Date()
sessionInfo()