## ----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, message = FALSE-------------------------------------------
library(scProkaR)
library(SingleCellExperiment)
library(SummarizedExperiment)

## ----knn-backends-------------------------------------------------------------
c(
    FNN = requireNamespace("FNN", quietly = TRUE),
    BiocNeighbors = requireNamespace("BiocNeighbors", quietly = TRUE)
)

## ----simulate-multidrug-------------------------------------------------------
set.seed(101)

sce <- simulate_tata_multidrug_sce(
    n_cells = 3000,
    n_features = 80,
    n_pcs = 15
)

sce$timepoint_factor <- factor(
    sce$timepoint,
    levels = sort(unique(sce$timepoint))
)

sce

## ----simulation-design--------------------------------------------------------
with(
    as.data.frame(colData(sce)),
    table(condition, timepoint)
)

## ----plot-simulation-time-----------------------------------------------------
PlotReduction(sce, reduction = "UMAP", colour_by = "timepoint_factor")

## ----plot-simulation-condition------------------------------------------------
PlotReduction(sce, reduction = "UMAP", colour_by = "condition")

## ----plot-simulation-branch---------------------------------------------------
PlotReduction(sce, reduction = "UMAP", colour_by = "simulated_branch")

## ----build-knn----------------------------------------------------------------
knn <- build_knn_graph(
    x = sce,
    dimred = "PCA",
    k = 30
)

c(
    n_cells = igraph::vcount(knn$graph),
    n_edges = igraph::ecount(knn$graph),
    k = knn$k
)

## ----knn-engine---------------------------------------------------------------
knn$knn_engine

## ----cluster-states-----------------------------------------------------------
sce$tata_cluster_step <- cluster_graph_states(
    cell_graph = knn$graph,
    method = "louvain"
)

table(sce$tata_cluster_step)

## ----plot-cluster-states------------------------------------------------------
PlotReduction(sce, reduction = "UMAP", colour_by = "tata_cluster_step")

## ----cluster-silhouette-------------------------------------------------------
if (requireNamespace("cluster", quietly = TRUE)) {
    set.seed(101)
    idx <- sample(seq_len(ncol(sce)), 600)
    pcs <- reducedDim(sce, "PCA")[idx, seq_len(5)]
    sil <- cluster::silhouette(
        as.integer(sce$tata_cluster_step[idx]),
        stats::dist(pcs)
    )
    round(mean(sil[, "sil_width"]), 3)
} else {
    message("Install 'cluster' for the optional silhouette diagnostic.")
}

## ----topology-weights---------------------------------------------------------
topology_table <- compute_topology_weights(
    adjacency = knn$adjacency,
    clusters = sce$tata_cluster_step
)

head(
    topology_table[
        order(-topology_table$topology_weight),
        c(
            "cluster_a", "cluster_b", "observed_edges",
            "expected_edges", "topology_weight"
        )
    ],
    6
)

## ----time-weights-------------------------------------------------------------
time_table <- compute_time_weights(
    adjacency = knn$adjacency,
    clusters = sce$tata_cluster_step,
    timepoint = sce$timepoint
)

head(
    time_table[
        order(-time_table$time_consistency),
        c(
            "cluster_a", "cluster_b", "n_time_edges",
            "flow_score", "time_weight", "time_consistency"
        )
    ],
    6
)

## ----plot-edge-score-scatter--------------------------------------------------
edge_metric_table <- merge(
    topology_table,
    time_table,
    by = c("cluster_a", "cluster_b"),
    all.x = TRUE,
    sort = FALSE
)

ggplot2::ggplot(
    edge_metric_table,
    ggplot2::aes(
        x = topology_weight,
        y = flow_score,
        size = observed_edges,
        colour = time_consistency
    )
) +
    ggplot2::geom_hline(
        yintercept = c(-0.10, 0.10),
        linetype = 2,
        colour = "grey60"
    ) +
    ggplot2::geom_point(alpha = 0.8) +
    ggplot2::scale_colour_viridis_c(option = "plasma") +
    ggplot2::theme_classic() +
    ggplot2::labs(
        title = "Topology and time consistency for cluster pairs",
        x = "topology weight",
        y = "temporal flow score",
        colour = "time consistency",
        size = "observed edges"
    )

## ----build-graph-manual-------------------------------------------------------
tata_manual <- build_tata_graph(
    topology_table = topology_table,
    time_table = time_table,
    clusters = sce$tata_cluster_step,
    embedding = reducedDim(sce, "PCA"),
    timepoint = sce$timepoint,
    alpha = 1,
    beta = 1,
    direction_threshold = 0.10,
    prune_threshold = 0.001,
    time_weight_mode = "directional_confidence"
)

head(
    tata_manual$edge_table[
        order(-tata_manual$edge_table$tata_weight),
        c(
            "cluster_a", "cluster_b", "topology_weight", "flow_score",
            "time_consistency", "tata_weight", "direction_label", "kept"
        )
    ],
    8
)

## ----plot-manual-graph-embedded-----------------------------------------------
plot_tata_cluster_graph(
    list(
        cluster_graph = tata_manual$graph,
        sce = sce,
        parameters = list(cluster_col = "tata_cluster_step")
    ),
    dimred = "UMAP",
    layout_mode = "embedding",
    show_cells = TRUE,
    colour_by = "tata_cluster_step",
    node_colour_by = "cluster",
    main = "Manual TATA graph aligned to UMAP"
)

## ----plot-manual-graph-abstract-----------------------------------------------
plot_tata_cluster_graph(
    list(
        cluster_graph = tata_manual$graph,
        sce = sce,
        parameters = list(cluster_col = "tata_cluster_step")
    ),
    layout_mode = "graph",
    node_colour_by = "median_time",
    main = "Manual TATA abstract graph"
)

## ----manual-pseudotime--------------------------------------------------------
sce$tata_cluster_step <- factor(sce$tata_cluster_step)

pt_manual <- compute_tata_pseudotime(
    cluster_graph = tata_manual$graph,
    sce = sce,
    cluster_col = "tata_cluster_step",
    time_col = "timepoint"
)

pt_manual$root_cluster
summary(pt_manual$cell_pseudotime_scaled)

## ----manual-root-check--------------------------------------------------------
tapply(
    sce$timepoint,
    sce$tata_cluster_step,
    median
)[pt_manual$root_cluster]

## ----run-tata-----------------------------------------------------------------
set.seed(101)

tata <- run_tata(
    sce = sce,
    dimred = "PCA",
    time_col = "timepoint",
    cluster_col = "tata_cluster",
    k = 30,
    cluster_method = "louvain",
    alpha = 1,
    beta = 1,
    direction_threshold = 0.10,
    prune_threshold = 0.001,
    time_weight_mode = "directional_confidence",
    do_pseudotime = TRUE,
    expected_branches = 3
)

names(tata)

## ----inspect-coldata----------------------------------------------------------
sce_tata <- tata$sce
sce_tata$tata_cluster <- factor(sce_tata$tata_cluster)

grep("^tata_", colnames(colData(sce_tata)), value = TRUE)
reducedDimNames(sce_tata)

## ----plot-tata-clusters-------------------------------------------------------
PlotReduction(sce_tata, reduction = "UMAP", colour_by = "tata_cluster")

## ----plot-tata-pseudotime-----------------------------------------------------
PlotReduction(
    sce_tata,
    reduction = "UMAP",
    colour_by = "tata_pseudotime_scaled"
)

## ----plot-tata-graph-embedded-------------------------------------------------
plot_tata_cluster_graph(
    tata,
    dimred = "UMAP",
    layout_mode = "embedding",
    show_cells = TRUE,
    colour_by = "timepoint",
    node_colour_by = "cluster",
    main = "TATA abstract graph aligned to UMAP"
)

## ----plot-tata-graph-abstract-------------------------------------------------
plot_tata_cluster_graph(
    tata,
    layout_mode = "graph",
    node_colour_by = "median_time",
    main = "TATA abstract graph"
)

## ----plot-trajectory-all-paths------------------------------------------------
plot_tata_trajectory_embedding(
    tata,
    dimred = "UMAP",
    colour_by = "timepoint",
    smooth_paths = TRUE
)

## ----plot-trajectory-top-paths------------------------------------------------
plot_tata_trajectory_embedding(
    tata,
    dimred = "UMAP",
    colour_by = "timepoint",
    max_paths = 3,
    smooth_paths = TRUE
)

## ----run-tata-topology-only---------------------------------------------------
set.seed(101)

tata_topology <- run_tata(
    sce = sce,
    dimred = "PCA",
    use_time = FALSE,
    time_col = NULL,
    cluster_col = "tata_cluster_topology",
    k = 30,
    cluster_method = "louvain",
    alpha = 1,
    beta = 1,
    direction_threshold = 0.10,
    prune_threshold = 0.001,
    do_pseudotime = TRUE,
    expected_branches = 3
)

## ----orient-pseudotime--------------------------------------------------------
pt_tata <- orient_pseudotime_to_time(
    pseudotime = tata$cell_pseudotime,
    timepoint = tata$sce$timepoint
)

pt_topology <- orient_pseudotime_to_time(
    pseudotime = tata_topology$cell_pseudotime,
    timepoint = tata_topology$sce$timepoint
)

## ----benchmark-local-order----------------------------------------------------
local_tata <- local_temporal_order_metrics(
    pseudotime = pt_tata,
    adjacency = tata$adjacency_matrix,
    timepoint = tata$sce$timepoint
)

local_topology <- local_temporal_order_metrics(
    pseudotime = pt_topology,
    adjacency = tata_topology$adjacency_matrix,
    timepoint = tata_topology$sce$timepoint
)

rbind(TATA = local_tata, topology_only = local_topology)

## ----benchmark-adjacent-auc---------------------------------------------------
auc_tata <- adjacent_timepoint_auc(
    pseudotime = pt_tata,
    timepoint = tata$sce$timepoint,
    group = tata$sce$condition
)

auc_topology <- adjacent_timepoint_auc(
    pseudotime = pt_topology,
    timepoint = tata_topology$sce$timepoint,
    group = tata_topology$sce$condition
)

rbind(TATA = auc_tata$summary, topology_only = auc_topology$summary)

## ----benchmark-auc-pairs------------------------------------------------------
head(auc_tata$pair_table)

## ----benchmark-reachability---------------------------------------------------
reach_tata <- directed_reachability_concordance(
    cluster_graph = tata$cluster_graph,
    clusters = tata$sce$tata_cluster,
    timepoint = tata$sce$timepoint,
    min_time_gap = 24
)

reach_tata$summary
head(reach_tata$pair_table)

## ----benchmark-anticausal-----------------------------------------------------
anticausal_tata <- anticausal_edge_mass(
    cluster_graph = tata$cluster_graph,
    clusters = tata$sce$tata_cluster,
    timepoint = tata$sce$timepoint,
    tolerance = 24
)

anticausal_tata$summary

## ----benchmark-anticausal-edges-----------------------------------------------
head(
    anticausal_tata$edge_table[
        order(anticausal_tata$edge_table$time_delta),
        c("from", "to", "weight", "time_delta", "temporal_class")
    ]
)

## ----benchmark-branch-onset---------------------------------------------------
true_onset <- c(drug_A = 24, drug_B = 96, drug_C = 48)

onset_tata <- branch_onset_metrics(
    pseudotime = pt_tata,
    branch = tata$sce$simulated_branch,
    branch_activation = tata$sce$branch_activation,
    true_onset = true_onset
)

onset_topology <- branch_onset_metrics(
    pseudotime = pt_topology,
    branch = tata_topology$sce$simulated_branch,
    branch_activation = tata_topology$sce$branch_activation,
    true_onset = true_onset
)

rbind(TATA = onset_tata$summary, topology_only = onset_topology$summary)

## ----benchmark-onset-table----------------------------------------------------
onset_tata$onset_table

## ----plot-calibration-tata----------------------------------------------------
plot_tata_time_calibration(
    tata,
    branch_col = "condition",
    summary = "median_iqr"
)

## ----plot-calibration-topology------------------------------------------------
plot_tata_time_calibration(
    tata_topology,
    time_col = "timepoint",
    branch_col = "condition",
    summary = "median_iqr"
)

## ----plot-edge-causality------------------------------------------------------
plot_tata_edge_causality(tata)

## ----branch-probabilities-----------------------------------------------------
branch_result <- compute_tata_branch_probabilities(
    tata_result = tata,
    expected_branches = 3
)

branch_result$terminal_clusters

## ----cluster-branch-table-----------------------------------------------------
prob_cols <- grep(
    "^tata_prob_",
    colnames(tata$cluster_branch_probabilities),
    value = TRUE
)

head(
    tata$cluster_branch_probabilities[
        ,
        c("cluster", prob_cols, "tata_branch_entropy",
            "tata_branch_assignment")
    ]
)

## ----plot-branch-entropy------------------------------------------------------
PlotReduction(
    sce_tata,
    reduction = "UMAP",
    colour_by = "tata_branch_entropy"
)

## ----plot-branch-assignment---------------------------------------------------
PlotReduction(
    sce_tata,
    reduction = "UMAP",
    colour_by = "tata_branch_assignment"
)

## ----plot-branch-probabilities------------------------------------------------
plot_tata_branch_probabilities(
    tata,
    dimred = "UMAP",
    ncol = 2
)

## ----select-branch-cells------------------------------------------------------
branch_one <- tata$terminal_clusters[1]

soft_selection <- select_tata_branch_cells(
    tata,
    branch = branch_one,
    selection_mode = "soft",
    probability_threshold = 0.60
)

hard_selection <- select_tata_branch_cells(
    tata,
    branch = branch_one,
    selection_mode = "hard"
)

c(
    branch = branch_one,
    soft = sum(soft_selection$selected),
    hard = sum(hard_selection$selected)
)

## ----plot-branch-selection----------------------------------------------------
plot_tata_branch_selection(
    tata,
    branch = branch_one,
    dimred = "UMAP",
    selection_mode = "soft",
    probability_threshold = 0.60
)

## ----plot-terminal-probabilities----------------------------------------------
plot_tata_terminal_probabilities(
    tata,
    n_cells = 4
)

## ----plot-branch-confusion----------------------------------------------------
plot_tata_branch_confusion(
    tata,
    truth_col = "simulated_branch",
    predicted_col = "tata_branch_assignment",
    normalize = "truth"
)

## ----select-trend-genes-------------------------------------------------------
feature_modules <- rowData(sce_tata)$feature_module

trend_genes <- unique(c(
    rownames(sce_tata)[feature_modules == "branch_A"][1:2],
    rownames(sce_tata)[feature_modules == "branch_B"][1:2],
    rownames(sce_tata)[feature_modules == "branch_C"][1:2],
    rownames(sce_tata)[feature_modules == "late"][1:2]
))

trend_genes

## ----compute-gene-trends------------------------------------------------------
tata_trends <- compute_tata_gene_trends(
    tata,
    genes = trend_genes,
    grid_length = 60
)

tata_trends$branches
head(tata_trends$trend_table)

## ----plot-gene-trends---------------------------------------------------------
plot_tata_gene_trends(
    tata_trends,
    ncol = 3
)

## ----cluster-gene-trends------------------------------------------------------
set.seed(101)

trend_clusters <- cluster_tata_gene_trends(
    tata_trends,
    n_clusters = 3
)

trend_clusters

## ----plot-gene-trend-heatmap--------------------------------------------------
set.seed(101)

plot_tata_gene_trend_heatmap(
    tata_trends,
    scale_rows = TRUE,
    n_clusters = 3
)

## ----tata-overview------------------------------------------------------------
tata_overview <- plot_tata_results(
    tata,
    dimred = "UMAP",
    n_probability_cells = 4
)

names(tata_overview)

## ----plot-overview-pseudotime-------------------------------------------------
tata_overview$pseudotime

## ----plot-overview-entropy----------------------------------------------------
tata_overview$entropy

## ----plot-overview-calibration------------------------------------------------
tata_overview$time_calibration

## ----plot-overview-branch-probabilities---------------------------------------
tata_overview$branch_probabilities

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