diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..4b43dad --- /dev/null +++ b/.gitignore @@ -0,0 +1,15 @@ +CLAUDE.md + +# macOS resource forks and metadata +.DS_Store +._.DS_Store +._* + +# vignette build artifacts (rendered HTMLs in vignettes/ dir) +vignettes/*.html + +# Claude Code internal directory +.claude/ + +# scratch / legacy +old/ diff --git a/DESCRIPTION b/DESCRIPTION index dc562fa..204f40b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: compact Title: compact -Version: 0.0.3 +Version: 0.1.0 Authors@R: c( person("Sam", "Morabito", , "smorabit@uci.edu", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7768-4856")), @@ -28,6 +28,7 @@ Imports: xgboost (>= 1.7.11.1), SHAPforxgboost (>= 0.1.3), ggplot2, + ggrepel, data.table, ggrastr, pbapply, diff --git a/NAMESPACE b/NAMESPACE index 91f56ed..109cd7c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -16,11 +16,15 @@ export(FindShapKeyDriver) export(HeatmapDistance) export(LoadAlertSystem) export(LoadAlertSystemPathways) +export(MacrostateTransitions) export(ModelZINB) export(ModulePerturbation) export(PerturbationLog2FC) export(PerturbationTransitions) export(PerturbationVectors) +export(PlotMacrostateTransitions) +export(PlotMarkovEmbedding) +export(PlotMarkovScatter) export(PlotTransitionVectors) export(PredictAttractors) export(PredictCommitment) @@ -52,13 +56,28 @@ importFrom(dplyr,rename) importFrom(ggplot2,aes) importFrom(ggplot2,coord_fixed) importFrom(ggplot2,element_blank) +importFrom(ggplot2,element_line) +importFrom(ggplot2,element_rect) importFrom(ggplot2,element_text) +importFrom(ggplot2,geom_hline) +importFrom(ggplot2,geom_point) +importFrom(ggplot2,geom_smooth) +importFrom(ggplot2,geom_text) importFrom(ggplot2,geom_tile) +importFrom(ggplot2,geom_vline) importFrom(ggplot2,ggplot) importFrom(ggplot2,ggtitle) +importFrom(ggplot2,labs) +importFrom(ggplot2,scale_color_discrete) +importFrom(ggplot2,scale_color_gradientn) +importFrom(ggplot2,scale_color_manual) importFrom(ggplot2,scale_fill_gradientn) +importFrom(ggplot2,scale_x_discrete) +importFrom(ggplot2,scale_y_discrete) importFrom(ggplot2,theme) +importFrom(ggplot2,theme_bw) importFrom(ggplot2,theme_minimal) +importFrom(ggrepel,geom_label_repel) importFrom(igraph,components) importFrom(igraph,graph_from_adjacency_matrix) importFrom(methods,as) diff --git a/NEWS.md b/NEWS.md index 28e3af1..c4180ee 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,14 @@ +# 0.1.0 (05-20-2026) +## Added +- `MacrostateTransitions` and `PlotMacrostateTransitions`: coarse-grain the cell-level transition probability matrix into a K × K group-level summary, with a stability index on the diagonal and a heatmap visualization. +- `CombinePerturbAssays`: combine multiple perturbation assays by additive delta superposition into a single composite assay. +- New vignette: **Basics — Simulated Branching Trajectory** (`simulation_tutorial`) covering the full `ModulePerturbation` workflow, graph connectivity checks, expression inspection, vector fields, macrostate transitions, and Markov chain analyses. +- New vignette: **Advanced Use Cases — NSCLC CD8+ T Cells** (`advanced_NSCLC`) demonstrating step-by-step use of `ApplyPerturbation`, `ApplyPropagation`, and `PerturbationTransitions` with a custom gene set, a pseudobulk Pearson correlation network, `CustomPerturbation`, and `CombinePerturbAssays`. +- pkgdown documentation site updated with both new vignettes. + +## Changes +- None + # 0.0.3 (02-20-2025) ## Added - New functions for cell fate analysis: `PredictPerturbationTime`, `PredictAttractors`, `PredictFates`, and `PredictCommitment`. diff --git a/R/MacrostateTransitions.R b/R/MacrostateTransitions.R new file mode 100644 index 0000000..7777d91 --- /dev/null +++ b/R/MacrostateTransitions.R @@ -0,0 +1,510 @@ +#' MacrostateTransitions +#' +#' @description +#' Coarse-grains the cell-level perturbation transition probability (TP) matrix +#' into a group-level summary matrix (K × K), analogous to the macrostate +#' transition probability heatmap in CellRank (Lange et al. 2022, Nature Methods, +#' Fig. 2c). For each pair of user-defined cell groups, computes the average +#' transition probability flowing from cells in the source group to cells in the +#' destination group under the specified perturbation. +#' +#' The diagonal of the resulting matrix gives a \emph{stability index} per group: +#' the mean probability that a cell in that group transitions to another cell +#' \emph{within the same group}. High stability (close to 1) indicates +#' attractor-like behavior — the perturbation does not move cells out of this +#' group. Low stability indicates a transient source state from which cells are +#' being redirected toward other groups. +#' +#' Unlike \code{PlotTransitionVectors}, this analysis is entirely independent of +#' the low-dimensional embedding. All computations operate on the full-dimensional +#' transition matrix, avoiding the distortions introduced by UMAP/t-SNE +#' projections. +#' +#' @param seurat_obj A Seurat object with a computed perturbation transition +#' probability matrix stored in \code{seurat_obj@graphs} (produced by +#' \code{ModulePerturbation}, \code{TFPerturbation}, or +#' \code{PerturbationTransitions}). +#' @param perturbation_name Character. The base name of the perturbation. The +#' function looks up the TP graph as \code{paste0(perturbation_name, '_tp')}. +#' @param graph Character. Name of the KNN graph stored in +#' \code{seurat_obj@graphs} (e.g., \code{"RNA_nn"}). Used to mask transition +#' probabilities to the biological manifold, identical to the masking applied +#' in \code{PredictAttractors} and \code{PredictPerturbationTime}. +#' @param group.by Character. Metadata column containing group labels used for +#' coarse-graining (e.g., \code{"functional_anno"}, \code{"seurat_clusters"}). +#' @param normalize_rows Logical. If \code{TRUE} (default), each row of the +#' coarse-grained matrix Q is normalized to sum to 1, making it row-stochastic. +#' If \code{FALSE}, rows sum to the fraction of total outflow captured within +#' the non-dropped groups (useful for diagnosing how much flow escapes excluded +#' small groups). +#' @param min_group_size Integer. Groups with fewer than this many cells are +#' excluded before computing Q, with a warning. Default: \code{5}. +#' @param store_result Logical. If \code{TRUE} (default), the result list is +#' stored in \code{seurat_obj@misc$MacrostateTransitions[[result_name]]}. +#' @param result_name Character. Key used for \code{@misc} storage. Defaults to +#' \code{perturbation_name}. +#' @param verbose Logical. Print progress messages. Default: \code{TRUE}. +#' +#' @return The Seurat object (if \code{store_result = TRUE}) with the following +#' list stored in \code{seurat_obj@misc$MacrostateTransitions[[result_name]]}: +#' \describe{ +#' \item{\code{Q}}{K × K numeric matrix. Coarse-grained (row-stochastic) +#' transition probability matrix. \code{Q[k1, k2]} is the average +#' probability of a cell in group \code{k1} transitioning to a cell in +#' group \code{k2}.} +#' \item{\code{stability}}{Named numeric vector of length K. Diagonal of Q; +#' the self-transition probability (stability index) for each group.} +#' \item{\code{group_sizes}}{Named integer vector. Number of cells per group +#' after small-group filtering.} +#' \item{\code{group.by}}{Character. The metadata column used.} +#' \item{\code{perturbation_name}}{Character. The perturbation name.} +#' } +#' If \code{store_result = FALSE}, returns the result list directly. +#' +#' @details +#' \strong{Mathematical summary:} +#' +#' Let P be the n × n row-stochastic cell-level transition matrix (after +#' transposing the stored column-stochastic TP, applying KNN masking, and +#' row-normalizing). Let M be the n × K sparse indicator matrix where +#' \code{M[i, k] = 1} if cell i belongs to group k. Then: +#' +#' \deqn{Q = \mathrm{diag}(1/n_k) \cdot M^\top \cdot P \cdot M} +#' +#' where \eqn{n_k} is the number of cells in group k. Q is row-stochastic +#' (each row sums to 1) because P is row-stochastic and M sums each row over +#' all K groups. +#' +#' \strong{Complexity:} O(n · k_nn · K). The bottleneck is \code{P \%*\% M} +#' (sparse n × n times sparse n × K). For n = 150,000, k_nn = 20, K = 20 this +#' requires ~3M floating-point operations and ~24 MB of sparse storage. No dense +#' n × n matrix is ever materialized. +#' +#' @seealso \code{\link{PlotMacrostateTransitions}}, +#' \code{\link{PredictAttractors}}, \code{\link{PredictFates}}, +#' \code{\link{VectorFieldCoherence}} +#' +#' @references +#' Lange, M. et al. CellRank for directed single-cell fate mapping. +#' \emph{Nature Methods} \strong{19}, 159–170 (2022). +#' \doi{10.1038/s41592-021-01346-6} +#' +#' @export +MacrostateTransitions <- function( + seurat_obj, + perturbation_name, + graph, + group.by, + normalize_rows = TRUE, + min_group_size = 5, + store_result = TRUE, + result_name = NULL, + verbose = TRUE +) { + + # ------------------------------------------------------------------------- + # input validation + # ------------------------------------------------------------------------- + + tp_graph_name <- paste0(perturbation_name, '_tp') + + if (!(graph %in% names(seurat_obj@graphs))) { + stop(paste0( + "Error: KNN graph '", graph, + "' not found in seurat_obj@graphs. ", + "Available graphs: ", paste(names(seurat_obj@graphs), collapse = ", ") + )) + } + + if (!(tp_graph_name %in% names(seurat_obj@graphs))) { + stop(paste0( + "Error: Transition probability graph '", tp_graph_name, + "' not found in seurat_obj@graphs. ", + "Did you run ModulePerturbation() or PerturbationTransitions()?", + " Available graphs: ", paste(names(seurat_obj@graphs), collapse = ", ") + )) + } + + if (!(group.by %in% colnames(seurat_obj@meta.data))) { + stop(paste0( + "Error: Metadata column '", group.by, + "' not found in seurat_obj@meta.data. ", + "Available columns: ", + paste(colnames(seurat_obj@meta.data), collapse = ", ") + )) + } + + if (is.null(result_name)) { + result_name <- perturbation_name + } + + # ------------------------------------------------------------------------- + # prepare group labels and filter small groups + # ------------------------------------------------------------------------- + + group_labels <- seurat_obj@meta.data[[group.by]] + + # convert to character to avoid factor level issues + group_labels <- as.character(group_labels) + + # compute group sizes + group_tab <- table(group_labels) + all_groups <- names(group_tab) + + # identify and warn about small groups + small_groups <- names(group_tab[group_tab < min_group_size]) + if (length(small_groups) > 0) { + warning(paste0( + "The following groups have fewer than ", min_group_size, + " cells and will be excluded: ", + paste(small_groups, collapse = ", ") + )) + } + + keep_groups <- names(group_tab[group_tab >= min_group_size]) + if (length(keep_groups) == 0) { + stop(paste0( + "Error: No groups with >= ", min_group_size, + " cells remain after filtering. ", + "Reduce min_group_size or check your group.by column." + )) + } + + # filter cells to kept groups + keep_cells <- which(group_labels %in% keep_groups) + group_labels_keep <- group_labels[keep_cells] + group_levels <- sort(unique(group_labels_keep)) + K <- length(group_levels) + + if (verbose) { + message(paste0( + "MacrostateTransitions: ", K, " groups retained from '", group.by, + "' (", length(keep_cells), " of ", ncol(seurat_obj), " cells)." + )) + } + + # ------------------------------------------------------------------------- + # prepare the row-stochastic transition matrix P + # ------------------------------------------------------------------------- + + if (verbose) message("Loading and preparing transition probability matrix...") + + # load KNN graph; set self-loops to 1 (consistent with PredictAttractors) + cell_graph <- SeuratObject::Graphs(seurat_obj, slot = graph) + cell_graph <- Matrix::Matrix(cell_graph, sparse = TRUE) + diag(cell_graph) <- 1 + + # load the stored TP (column-stochastic: tp[i,j] = P(j -> i)) + tp <- SeuratObject::Graphs(seurat_obj, slot = tp_graph_name) + tp <- Matrix::Matrix(tp, sparse = TRUE) + tp[is.na(tp)] <- 0 + + # transpose to row-stochastic: P[i,j] = P(i -> j) + tp <- Matrix::t(tp) + + # mask to KNN neighborhood + tp <- tp * cell_graph + + # row-normalize + row_sums <- Matrix::rowSums(tp) + row_sums[row_sums == 0] <- 1 + P <- tp / row_sums + + # subset to cells in kept groups (rows and columns) + P <- P[keep_cells, keep_cells] + + # re-normalize rows after subsetting (some outflow to excluded cells is lost) + row_sums2 <- Matrix::rowSums(P) + row_sums2[row_sums2 == 0] <- 1 + P <- P / row_sums2 + + # ------------------------------------------------------------------------- + # build sparse membership matrix M (n_keep × K) + # ------------------------------------------------------------------------- + + if (verbose) message("Building membership matrix and computing coarse-grained Q...") + + group_idx <- as.integer(factor(group_labels_keep, levels = group_levels)) + n_keep <- length(keep_cells) + + M <- Matrix::sparseMatrix( + i = seq_len(n_keep), + j = group_idx, + x = 1, + dims = c(n_keep, K), + dimnames = list(NULL, group_levels) + ) + + # ------------------------------------------------------------------------- + # compute Q = diag(1/n_k) * M^T * P * M + # ------------------------------------------------------------------------- + + group_sizes <- Matrix::colSums(M) # K-vector of group sizes + M_norm <- Matrix::Diagonal(x = 1 / group_sizes) %*% Matrix::t(M) # K × n + + # two sparse multiplications: P %*% M is the expensive step (n × K) + # then M_norm %*% (P %*% M) is K × n times n × K = K × K + PM <- P %*% M # n × K sparse + Q <- as.matrix(M_norm %*% PM) # K × K dense (tiny) + + rownames(Q) <- group_levels + colnames(Q) <- group_levels + + # optional row re-normalization after subsetting may have left rows not + # summing exactly to 1; apply normalize_rows if requested + if (normalize_rows) { + Q_row_sums <- rowSums(Q) + Q_row_sums[Q_row_sums == 0] <- 1 + Q <- Q / Q_row_sums + } + + # ------------------------------------------------------------------------- + # compute stability index and package result + # ------------------------------------------------------------------------- + + stability <- diag(Q) + names(stability) <- group_levels + + group_sizes_int <- as.integer(group_sizes) + names(group_sizes_int) <- group_levels + + result <- list( + Q = Q, + stability = stability, + group_sizes = group_sizes_int, + group.by = group.by, + perturbation_name = perturbation_name + ) + + if (verbose) { + message("Done. Stability index summary:") + message(paste( + sprintf(" %-30s %.3f (n = %d)", group_levels, stability, group_sizes_int), + collapse = "\n" + )) + } + + # ------------------------------------------------------------------------- + # store result and return + # ------------------------------------------------------------------------- + + if (store_result) { + if (is.null(seurat_obj@misc$MacrostateTransitions)) { + seurat_obj@misc$MacrostateTransitions <- list() + } + seurat_obj@misc$MacrostateTransitions[[result_name]] <- result + return(seurat_obj) + } else { + return(result) + } +} + + +#' PlotMacrostateTransitions +#' +#' @description +#' Visualizes the coarse-grained group-level transition probability matrix +#' computed by \code{\link{MacrostateTransitions}} as a heatmap, analogous to +#' the macrostate transition probability heatmap in CellRank (Lange et al. 2022, +#' Fig. 2c). +#' +#' Each tile \code{[source, dest]} is colored by the average transition +#' probability flowing from the source group to the destination group under the +#' specified perturbation. The diagonal tiles represent self-transitions +#' (stability index) and are optionally annotated with the numeric stability +#' value and a heavier border for visual emphasis. +#' +#' @param seurat_obj A Seurat object with results stored in +#' \code{seurat_obj@misc$MacrostateTransitions} (produced by +#' \code{\link{MacrostateTransitions}}). +#' @param perturbation_name Character. Name of the perturbation to visualize. +#' Used to look up \code{result_name} in \code{@misc$MacrostateTransitions}. +#' @param result_name Character. Key to retrieve from +#' \code{seurat_obj@misc$MacrostateTransitions}. Defaults to +#' \code{perturbation_name}. +#' @param group_order Character vector. Order of groups on both axes. Defaults +#' to alphabetical order (the order stored in the result). Both the x-axis +#' (destination) and the y-axis (source, reversed) follow this order. +#' @param color_scale Character vector. Two or more colors defining the +#' low-to-high color gradient, passed to +#' \code{ggplot2::scale_fill_gradientn}. Default: +#' \code{c("white", "#2166AC")} (white to dark blue). For a diverging scale +#' (e.g., when plotting ΔQ between perturbations), supply three colors with +#' the midpoint as the second element. +#' @param show_stability Logical. If \code{TRUE} (default), overlay the +#' numeric stability index on each diagonal tile. +#' @param stability_size Numeric. Font size for the stability text annotation. +#' Default: \code{3}. +#' @param diagonal_border Logical. If \code{TRUE} (default), draw a heavier +#' black border around the diagonal tiles to visually highlight +#' self-transitions. +#' @param text_size Numeric. Base text size for axis labels and theme. Default: +#' \code{10}. +#' @param legend_title Character. Title for the fill legend. Default: +#' \code{"Transition\nProbability"}. +#' @param title Character. Plot title. Defaults to \code{perturbation_name}. +#' +#' @return A \code{ggplot2} object. +#' +#' @seealso \code{\link{MacrostateTransitions}}, \code{\link{PredictAttractors}}, +#' \code{\link{VectorFieldCoherence}} +#' +#' @references +#' Lange, M. et al. CellRank for directed single-cell fate mapping. +#' \emph{Nature Methods} \strong{19}, 159–170 (2022). +#' \doi{10.1038/s41592-021-01346-6} +#' +#' @importFrom ggplot2 ggplot aes geom_tile geom_text scale_fill_gradientn +#' scale_x_discrete scale_y_discrete labs theme_bw theme element_text +#' element_blank element_line element_rect +#' +#' @export +PlotMacrostateTransitions <- function( + seurat_obj, + perturbation_name, + result_name = NULL, + group_order = NULL, + color_scale = c("white", "#2166AC"), + show_stability = TRUE, + stability_size = 3, + diagonal_border = TRUE, + text_size = 10, + legend_title = "Transition\nProbability", + title = NULL +) { + + # ------------------------------------------------------------------------- + # retrieve stored result + # ------------------------------------------------------------------------- + + if (is.null(result_name)) { + result_name <- perturbation_name + } + + if (is.null(seurat_obj@misc$MacrostateTransitions)) { + stop(paste0( + "Error: No MacrostateTransitions results found in seurat_obj@misc. ", + "Run MacrostateTransitions() first." + )) + } + + if (!(result_name %in% names(seurat_obj@misc$MacrostateTransitions))) { + stop(paste0( + "Error: Result '", result_name, "' not found in ", + "seurat_obj@misc$MacrostateTransitions. ", + "Available results: ", + paste(names(seurat_obj@misc$MacrostateTransitions), collapse = ", ") + )) + } + + res <- seurat_obj@misc$MacrostateTransitions[[result_name]] + Q <- res$Q + stability <- res$stability + + # ------------------------------------------------------------------------- + # group ordering + # ------------------------------------------------------------------------- + + all_groups <- rownames(Q) + + if (is.null(group_order)) { + group_order <- all_groups + } else { + missing_groups <- setdiff(group_order, all_groups) + if (length(missing_groups) > 0) { + warning(paste0( + "The following groups in 'group_order' are not in the result and will be ignored: ", + paste(missing_groups, collapse = ", ") + )) + } + group_order <- intersect(group_order, all_groups) + if (length(group_order) == 0) { + stop("Error: No valid groups remain after filtering 'group_order'. Check group names.") + } + # subset Q and stability to the requested groups (in case group_order is a subset) + Q <- Q[group_order, group_order, drop = FALSE] + stability <- stability[group_order] + } + + # ------------------------------------------------------------------------- + # reshape Q to long format for ggplot2 + # ------------------------------------------------------------------------- + + df_long <- expand.grid( + source = group_order, + dest = group_order, + stringsAsFactors = FALSE + ) + df_long$probability <- Q[cbind(df_long$source, df_long$dest)] + + # enforce factor ordering + df_long$source <- factor(df_long$source, levels = group_order) + df_long$dest <- factor(df_long$dest, levels = group_order) + + # diagonal data frame for annotation + diag_df <- data.frame( + source = factor(group_order, levels = group_order), + dest = factor(group_order, levels = group_order), + label = sprintf("%.2f", stability[group_order]), + stringsAsFactors = FALSE + ) + + # ------------------------------------------------------------------------- + # plot title + # ------------------------------------------------------------------------- + + if (is.null(title)) { + title <- perturbation_name + } + + # ------------------------------------------------------------------------- + # build heatmap + # ------------------------------------------------------------------------- + + p <- ggplot2::ggplot(df_long, ggplot2::aes(x = dest, y = source, fill = probability)) + + ggplot2::geom_tile(color = "grey85", linewidth = 0.3) + + ggplot2::scale_fill_gradientn( + colors = color_scale, + limits = c(0, max(Q, na.rm = TRUE)), + name = legend_title + ) + + ggplot2::scale_y_discrete(limits = rev(group_order)) + + ggplot2::scale_x_discrete(limits = group_order) + + ggplot2::labs( + x = "Destination group", + y = "Source group", + title = title + ) + + ggplot2::theme_bw(base_size = text_size) + + ggplot2::theme( + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1), + panel.grid = ggplot2::element_blank(), + plot.title = ggplot2::element_text(hjust = 0.5, face = "bold") + ) + + # stability text on diagonal tiles + if (show_stability) { + p <- p + ggplot2::geom_text( + data = diag_df, + mapping = ggplot2::aes(x = dest, y = source, label = label), + size = stability_size, + color = "white", + fontface = "bold", + inherit.aes = FALSE + ) + } + + # heavier border around diagonal tiles + if (diagonal_border) { + p <- p + ggplot2::geom_tile( + data = diag_df, + mapping = ggplot2::aes(x = dest, y = source), + fill = NA, + color = "black", + linewidth = 0.8, + inherit.aes = FALSE + ) + } + + return(p) +} diff --git a/R/MarkovPlots.R b/R/MarkovPlots.R new file mode 100644 index 0000000..b3ad7fd --- /dev/null +++ b/R/MarkovPlots.R @@ -0,0 +1,402 @@ + +#' Plot a Markov Score on a Dimensionality Reduction Embedding +#' +#' Visualizes a Markov analysis score (from \code{\link{PredictAttractors}}, +#' \code{\link{PredictPerturbationTime}}, \code{\link{PredictFates}}, or +#' \code{\link{PredictCommitment}}) as a colored scatter plot on any 2D +#' reduction stored in the Seurat object. Optionally marks source and/or sink +#' cell populations with a centroid diamond and greys out unconverged cells so +#' they do not contaminate the color scale. +#' +#' @param seurat_obj A Seurat object containing the dimensionality reduction +#' and the score in \code{@meta.data}. +#' @param feature Character. Metadata column name to color cells by (e.g., +#' \code{"attractor_score"}, \code{"perturbation_pseudotime"}, +#' \code{"commitment_score"}). +#' @param reduction Character. Name of the reduction stored in +#' \code{seurat_obj@reductions}. Default: \code{"umap"}. +#' @param source_cells Character vector or NULL. Cell barcodes belonging to the +#' source population. Their mean embedding coordinates are computed and drawn +#' as a single centroid diamond (shape 18). NULL = no source centroid. +#' @param sink_cells Character vector or NULL. Cell barcodes belonging to the +#' sink population. Their mean embedding coordinates are computed and drawn as +#' a single centroid diamond (shape 18). NULL = no sink centroid. +#' @param source_label Character. Text label shown next to the source centroid +#' via \code{ggrepel::geom_label_repel}. Only used when +#' \code{label_centroids = TRUE} and \code{source_cells} is not NULL. +#' Default: \code{"Source"}. +#' @param sink_label Character. Text label shown next to the sink centroid via +#' \code{ggrepel::geom_label_repel}. Only used when +#' \code{label_centroids = TRUE} and \code{sink_cells} is not NULL. +#' Default: \code{"Sink"}. +#' @param centroid_size Numeric. Size of the centroid diamond marker. +#' Default: \code{5}. +#' @param centroid_color Character. Color of the centroid diamond marker and its +#' label. Default: \code{"black"}. +#' @param label_centroids Logical. If \code{TRUE}, add a text label to each +#' centroid using \code{ggrepel::geom_label_repel}. Default: \code{TRUE}. +#' @param label_size Numeric. Font size for centroid labels. Default: \code{3}. +#' @param unconverged_cells Character vector or NULL. Cell barcodes that hit +#' \code{max_iter} in \code{\link{PredictPerturbationTime}}. These cells are +#' drawn in \code{unconverged_color} below all other layers and excluded from +#' the color scale limits. NULL = no special handling. +#' @param unconverged_color Character. Color for unconverged cells. +#' Default: \code{"grey70"}. +#' @param color_scale Character vector. Two or more colors for the sequential +#' color gradient passed to \code{scale_color_gradientn}. +#' Default: \code{c("lightgrey", "#2166AC")}. +#' @param pt_size Numeric. Point size for the main cell layer. Default: \code{0.8}. +#' @param alpha Numeric. Transparency for the main cell layer. Default: \code{0.8}. +#' @param legend_title Character or NULL. Legend title. If NULL, defaults to +#' \code{feature}. +#' @param title Character or NULL. Plot title. If NULL, defaults to +#' \code{feature}. +#' @param ... Additional arguments (currently unused; reserved for future use). +#' +#' @return A \code{ggplot2} object. +#' +#' @details +#' Layers are drawn bottom to top: unconverged cells (grey) → converged cells +#' (colored gradient) → source centroid diamond → sink centroid diamond → +#' centroid labels (ggrepel). A single centroid per group is computed from the +#' mean of the source/sink cells' embedding coordinates, so the marker position +#' represents the group center rather than cluttering the plot with one marker +#' per cell. +#' +#' The color scale limits are computed from converged cells only, so +#' unconverged cells (which often carry a meaningless \code{max_iter} value) +#' do not distort the range. +#' +#' @seealso \code{\link{PredictAttractors}}, \code{\link{PredictPerturbationTime}}, +#' \code{\link{PredictFates}}, \code{\link{PredictCommitment}}, +#' \code{\link{PlotMarkovScatter}} +#' +#' @importFrom ggplot2 ggplot aes geom_point scale_color_gradientn labs +#' @importFrom ggrepel geom_label_repel +#' +#' @export +PlotMarkovEmbedding <- function( + seurat_obj, + feature, + reduction = "umap", + source_cells = NULL, + sink_cells = NULL, + source_label = "Source", + sink_label = "Sink", + centroid_size = 5, + centroid_color = "black", + label_centroids = TRUE, + label_size = 3, + unconverged_cells = NULL, + unconverged_color = "grey70", + color_scale = c("lightgrey", "#2166AC"), + pt_size = 0.8, + alpha = 0.8, + legend_title = NULL, + title = NULL, + ... +) { + # validate feature + if (!feature %in% colnames(seurat_obj@meta.data)) { + stop(sprintf( + "feature '%s' not found in seurat_obj@meta.data.", + feature + )) + } + + # validate reduction + if (!reduction %in% names(seurat_obj@reductions)) { + stop(sprintf( + "reduction '%s' not found in seurat_obj@reductions. Available: %s", + reduction, paste(names(seurat_obj@reductions), collapse = ", ") + )) + } + + legend_title <- if (is.null(legend_title)) feature else legend_title + title <- if (is.null(title)) feature else title + + # extract embedding coordinates + emb <- as.data.frame(Seurat::Embeddings(seurat_obj, reduction = reduction)) + colnames(emb)[1:2] <- c("dim1", "dim2") + emb$cell <- rownames(emb) + emb$score <- seurat_obj@meta.data[rownames(emb), feature] + + # split unconverged from converged cells + if (!is.null(unconverged_cells)) { + emb_unconverged <- emb[emb$cell %in% unconverged_cells, ] + emb_converged <- emb[!emb$cell %in% unconverged_cells, ] + } else { + emb_unconverged <- emb[0, ] + emb_converged <- emb + } + + # color scale limits exclude unconverged cells + score_limits <- range(emb_converged$score, na.rm = TRUE) + + # base plot — no data yet so we can add layers in the correct order + p <- ggplot2::ggplot() + + # layer 1: unconverged cells (grey, bottom-most) + if (nrow(emb_unconverged) > 0) { + p <- p + ggplot2::geom_point( + data = emb_unconverged, + ggplot2::aes(x = dim1, y = dim2), + color = unconverged_color, + size = pt_size, + alpha = alpha + ) + } + + # layer 2: converged cells colored by score — sorted so high-value cells render on top + emb_converged <- emb_converged[order(emb_converged$score, na.last = FALSE), ] + + p <- p + + ggplot2::geom_point( + data = emb_converged, + ggplot2::aes(x = dim1, y = dim2, color = score), + size = pt_size, + alpha = alpha + ) + + ggplot2::scale_color_gradientn( + colors = color_scale, + limits = score_limits, + name = legend_title + ) + + # layers 3–4: source and sink centroid diamonds + centroid_rows <- list() + if (!is.null(source_cells)) { + src_in_obj <- source_cells[source_cells %in% rownames(emb)] + centroid_rows[["source"]] <- data.frame( + dim1 = mean(emb[src_in_obj, "dim1"]), + dim2 = mean(emb[src_in_obj, "dim2"]), + label = source_label, + stringsAsFactors = FALSE + ) + } + if (!is.null(sink_cells)) { + snk_in_obj <- sink_cells[sink_cells %in% rownames(emb)] + centroid_rows[["sink"]] <- data.frame( + dim1 = mean(emb[snk_in_obj, "dim1"]), + dim2 = mean(emb[snk_in_obj, "dim2"]), + label = sink_label, + stringsAsFactors = FALSE + ) + } + + if (length(centroid_rows) > 0) { + centroid_df <- do.call(rbind, centroid_rows) + + p <- p + ggplot2::geom_point( + data = centroid_df, + ggplot2::aes(x = dim1, y = dim2), + shape = 18, + size = centroid_size, + color = centroid_color, + inherit.aes = FALSE + ) + + if (label_centroids) { + p <- p + ggrepel::geom_label_repel( + data = centroid_df, + ggplot2::aes(x = dim1, y = dim2, label = label), + size = label_size, + color = centroid_color, + fill = "white", + box.padding = 0.5, + inherit.aes = FALSE + ) + } + } + + p <- p + + ggplot2::labs(title = title) + + hdWGCNA::umap_theme() + + p +} + + +#' Scatter Plot of Two Markov or Biological Features +#' +#' Plots any two metadata columns as a scatter plot, colored by an optional +#' grouping variable. Designed for visualizing relationships between Markov +#' analysis scores (e.g., attractor score vs. hitting time, commitment +#' probability vs. module eigengene) but works with any numeric metadata pair. +#' +#' @param seurat_obj A Seurat object with all requested features in +#' \code{@meta.data}. +#' @param x_feature Character. Metadata column for the x-axis. +#' @param y_feature Character. Metadata column for the y-axis. +#' @param color.by Character or NULL. Metadata column for point color. If +#' NULL, all points are drawn in a single color (\code{"#2166AC"}). +#' @param color_palette Named character vector or NULL. Named color vector for +#' categorical \code{color.by} groups passed to \code{scale_color_manual}. +#' If NULL, \code{ggplot2} default colors are used. Ignored when +#' \code{color.by} is numeric. +#' @param pt_size Numeric. Point size. Default: \code{0.8}. +#' @param alpha Numeric. Point transparency. Default: \code{0.6}. +#' @param add_smooth Logical. Add a smoothing line via \code{geom_smooth}. +#' Default: \code{FALSE}. +#' @param smooth_method Character. Smoothing method passed to +#' \code{geom_smooth} (e.g., \code{"loess"}, \code{"lm"}). +#' Default: \code{"loess"}. +#' @param smooth_se Logical. Show confidence interval band around the smooth +#' line. Default: \code{TRUE}. +#' @param hline Numeric or NULL. Y-value for a horizontal reference line. +#' NULL = no line. Drawn under the point layer. Default: \code{NULL}. +#' @param hline_color Character. Color for the horizontal reference line. +#' Default: \code{"grey40"}. +#' @param hline_linetype Character. Linetype for the horizontal reference line. +#' Default: \code{"dashed"}. +#' @param vline Numeric or NULL. X-value for a vertical reference line. +#' NULL = no line. Drawn under the point layer. Default: \code{NULL}. +#' @param vline_color Character. Color for the vertical reference line. +#' Default: \code{"grey40"}. +#' @param vline_linetype Character. Linetype for the vertical reference line. +#' Default: \code{"dashed"}. +#' @param x_label Character or NULL. X-axis label. If NULL, defaults to +#' \code{x_feature}. +#' @param y_label Character or NULL. Y-axis label. If NULL, defaults to +#' \code{y_feature}. +#' @param title Character or NULL. Plot title. +#' @param legend_title Character or NULL. Legend title. If NULL, defaults to +#' \code{color.by}. +#' +#' @return A \code{ggplot2} object. +#' +#' @details +#' Color handling adapts to the type of \code{color.by}: +#' \itemize{ +#' \item NULL: fixed single color. +#' \item Numeric column: sequential gradient via \code{scale_color_gradientn}. +#' \item Factor/character column: discrete colors via \code{scale_color_manual} +#' (if \code{color_palette} provided) or \code{scale_color_discrete}. +#' } +#' +#' Reference lines (\code{hline}, \code{vline}) are drawn before the point +#' layer so they sit under the points. +#' +#' @seealso \code{\link{PredictAttractors}}, \code{\link{PredictPerturbationTime}}, +#' \code{\link{PredictCommitment}}, \code{\link{PlotMarkovEmbedding}} +#' +#' @importFrom ggplot2 ggplot aes geom_point geom_hline geom_vline geom_smooth +#' scale_color_gradientn scale_color_manual scale_color_discrete labs theme_bw +#' +#' @export +PlotMarkovScatter <- function( + seurat_obj, + x_feature, + y_feature, + color.by = NULL, + color_palette = NULL, + pt_size = 0.8, + alpha = 0.6, + add_smooth = FALSE, + smooth_method = "loess", + smooth_se = TRUE, + hline = NULL, + hline_color = "grey40", + hline_linetype = "dashed", + vline = NULL, + vline_color = "grey40", + vline_linetype = "dashed", + x_label = NULL, + y_label = NULL, + title = NULL, + legend_title = NULL +) { + # validate x_feature + if (!x_feature %in% colnames(seurat_obj@meta.data)) { + stop(sprintf( + "x_feature '%s' not found in seurat_obj@meta.data.", + x_feature + )) + } + + # validate y_feature + if (!y_feature %in% colnames(seurat_obj@meta.data)) { + stop(sprintf( + "y_feature '%s' not found in seurat_obj@meta.data.", + y_feature + )) + } + + x_label <- if (is.null(x_label)) x_feature else x_label + y_label <- if (is.null(y_label)) y_feature else y_label + legend_title <- if (is.null(legend_title) && !is.null(color.by)) color.by else legend_title + + # build plotting data frame + df <- data.frame( + x = seurat_obj@meta.data[[x_feature]], + y = seurat_obj@meta.data[[y_feature]], + cell = colnames(seurat_obj), + stringsAsFactors = FALSE + ) + if (!is.null(color.by)) { + df$color_var <- seurat_obj@meta.data[[color.by]] + } + df <- df[!is.na(df$x) & !is.na(df$y), ] + + p <- ggplot2::ggplot(df, ggplot2::aes(x = x, y = y)) + + # reference lines drawn under points + if (!is.null(hline)) { + p <- p + ggplot2::geom_hline( + yintercept = hline, + color = hline_color, + linetype = hline_linetype + ) + } + if (!is.null(vline)) { + p <- p + ggplot2::geom_vline( + xintercept = vline, + color = vline_color, + linetype = vline_linetype + ) + } + + # point layer and color scale + if (is.null(color.by)) { + p <- p + ggplot2::geom_point( + color = "#2166AC", + size = pt_size, + alpha = alpha + ) + } else if (is.numeric(df$color_var)) { + p <- p + + ggplot2::geom_point( + ggplot2::aes(color = color_var), + size = pt_size, + alpha = alpha + ) + + ggplot2::scale_color_gradientn( + colors = c("lightgrey", "#2166AC"), + name = legend_title + ) + } else { + p <- p + ggplot2::geom_point( + ggplot2::aes(color = color_var), + size = pt_size, + alpha = alpha + ) + if (!is.null(color_palette)) { + p <- p + ggplot2::scale_color_manual( + values = color_palette, + name = legend_title + ) + } else { + p <- p + ggplot2::scale_color_discrete(name = legend_title) + } + } + + if (add_smooth) { + p <- p + ggplot2::geom_smooth(method = smooth_method, se = smooth_se) + } + + p <- p + + ggplot2::labs(x = x_label, y = y_label, title = title) + + ggplot2::theme_bw() + + p +} diff --git a/R/PlotTransitionVectors.R b/R/PlotTransitionVectors.R index 29a423a..a58132a 100644 --- a/R/PlotTransitionVectors.R +++ b/R/PlotTransitionVectors.R @@ -4,8 +4,9 @@ #' This function visualizes cell-state transitions in a Seurat object by plotting transition vectors on a #' dimensionality-reduced embedding (e.g., UMAP). It computes transition vectors based on a specified #' perturbation and overlays them on a scatter plot of the cells, allowing insights into cell-state shifts -#' following the perturbation. Three visualization modes are available: classic grid arrows -#' (\code{"arrows"}), streamlines (\code{"streamlines"}), or both overlaid (\code{"both"}). +#' following the perturbation. Four visualization modes are available: classic grid arrows +#' (\code{"arrows"}), streamlines (\code{"streamlines"}), both overlaid (\code{"both"}), or one arrow +#' per cell (\code{"cells"}). #' #' @param seurat_obj A Seurat object containing single-cell data, including metadata and embeddings. #' @param perturbation_name A character string specifying the name of the perturbation. This is used @@ -32,11 +33,15 @@ #' shorter arrows are more transparent (default: \code{TRUE}). #' @param use_velocyto Logical. Use the velocyto.R C backend for embedding-arrow calculation #' instead of the built-in \code{SparseEmbArrows} (default: \code{FALSE}). -#' @param plot_mode Character. One of \code{"arrows"} (default), \code{"streamlines"}, or -#' \code{"both"}. \code{"arrows"} reproduces the classic grid-arrow plot. \code{"streamlines"} -#' replaces arrows with continuous streamlines traced through the vector field. -#' \code{"both"} overlays streamlines on top of the arrow plot. Streamline rendering -#' requires the \pkg{ggvfields} package. +#' @param plot_mode Character. One of \code{"arrows"} (default), \code{"streamlines"}, +#' \code{"both"}, or \code{"cells"}. \code{"arrows"} reproduces the classic grid-arrow +#' plot. \code{"streamlines"} replaces arrows with continuous streamlines traced through +#' the vector field. \code{"both"} overlays streamlines on top of the arrow plot. +#' \code{"cells"} draws one arrow per cell at its exact embedding position, bypassing +#' grid aggregation entirely — the arrows themselves are coloured by \code{color.by} and +#' no background scatter plot is drawn. Recommended for small datasets (< ~2,000 cells) +#' where grid bins would be too sparse. Streamline rendering requires the \pkg{ggvfields} +#' package. #' @param stream_n Integer. Number of streamline seed points along each axis of the #' embedding grid passed to \code{ggvfields::geom_stream_field()} (default: 15). #' @param stream_L Numeric. Target length of each streamline, in embedding coordinate @@ -79,7 +84,7 @@ PlotTransitionVectors <- function( raster_dpi = 300, arrow_alpha = TRUE, use_velocyto = FALSE, - plot_mode = c("arrows", "streamlines", "both"), + plot_mode = c("arrows", "streamlines", "both", "cells"), stream_n = 15, stream_L = 2, stream_normalize = TRUE, @@ -126,6 +131,59 @@ PlotTransitionVectors <- function( emb <- Reductions(seurat_obj, reduction)@cell.embeddings[,1:2] colnames(emb) <- paste0('emb_', 1:ncol(emb)) + # ------------------------------------------------------------------------- + # cells mode: one arrow per cell, colored by color.by, no scatter layer + # ------------------------------------------------------------------------- + if (plot_mode == "cells") { + + # arsd is already capped at max_pct and normalized to [-1, 1]. + # Scale back to embedding coordinates using the grid-cell diagonal as + # a reference length (same reference grid_vectors uses), then apply + # arrow_scale as a final multiplier so the parameter remains meaningful. + cell_emb <- emb[rownames(ars), , drop = FALSE] + emb_limits <- apply(cell_emb, 2, range) + emb_delta <- (emb_limits[2, ] - emb_limits[1, ]) / grid_resolution + cell_scale <- sqrt(sum(emb_delta^2)) / 2 * arrow_scale + + cell_df <- data.frame( + x0 = ars$x0, + y0 = ars$y0, + x1 = ars$x0 + arsd$xd * cell_scale, + y1 = ars$y0 + arsd$yd * cell_scale, + value = seurat_obj@meta.data[rownames(ars), color.by], + vector.length = sqrt((arsd$xd * cell_scale)^2 + (arsd$yd * cell_scale)^2) + ) + + if (arrow_alpha) { + p <- ggplot(cell_df, aes(x = x0, y = y0)) + + geom_segment( + aes( + x = x0, y = y0, + xend = x1, yend = y1, + color = value, + alpha = vector.length + ), + arrow = grid::arrow(length = unit(0.1, "cm")), + linewidth = arrow_size + ) + + guides(alpha = "none") + } else { + p <- ggplot(cell_df, aes(x = x0, y = y0)) + + geom_segment( + aes( + x = x0, y = y0, + xend = x1, yend = y1, + color = value + ), + arrow = grid::arrow(length = unit(0.1, "cm")), + linewidth = arrow_size + ) + } + + p <- p + hdWGCNA::umap_theme() + return(p) + } + # build the grid-aggregated vector field grid.df <- grid_vectors(emb[rownames(ars),], arsd, resolution=grid_resolution) diff --git a/R/PredictCellFates.R b/R/PredictCellFates.R index 2e03e3c..318fcd7 100644 --- a/R/PredictCellFates.R +++ b/R/PredictCellFates.R @@ -219,14 +219,19 @@ PredictPerturbationTime <- function( #' \code{seurat_obj@graphs} (e.g., "RNA_nn") used to mask the transition probabilities. #' @param output_name Character. The column name used to store the resulting attractor #' scores in \code{seurat_obj@meta.data}. Default is \code{"attractor_score"}. -#' @param quantile_threshold Numeric. The percentile threshold used to explicitly define +#' @param quantile_threshold Numeric. The percentile threshold used to explicitly define #' which cells belong to the terminal sink state. Default is \code{0.98} (top 2%). -#' @param return_seurat Logical. If \code{TRUE}, returns the updated Seurat object. -#' If \code{FALSE}, returns a list containing the numeric attractor scores and a +#' @param rank_transform Logical. If \code{TRUE}, transforms the raw stationary +#' distribution values into percentile ranks using the empirical CDF (ECDF), +#' mapping scores to \code{[0, 1]} with better dynamic range for visualization. +#' The \code{quantile_threshold} is always applied to the raw stationary +#' distribution regardless of this setting. Default is \code{FALSE}. +#' @param return_seurat Logical. If \code{TRUE}, returns the updated Seurat object. +#' If \code{FALSE}, returns a list containing the numeric attractor scores and a #' character vector of the identified sink cells. Default is \code{TRUE}. #' -#' @return If \code{return_seurat = TRUE}, a Seurat object with the calculated attractor -#' scores appended to the metadata. If \code{return_seurat = FALSE}, a list containing +#' @return If \code{return_seurat = TRUE}, a Seurat object with the calculated attractor +#' scores appended to the metadata. If \code{return_seurat = FALSE}, a list containing #' \code{attractor_score} and \code{sink_cells}. #' #' @export @@ -236,6 +241,7 @@ PredictAttractors <- function( graph, output_name = "attractor_score", quantile_threshold = 0.98, + rank_transform = FALSE, return_seurat = TRUE ){ @@ -295,15 +301,21 @@ PredictAttractors <- function( stat_dist <- stat_dist / sum(stat_dist) names(stat_dist) <- rownames(P) - # define the sinks based on the stationary distribution - threshold <- quantile(stat_dist, quantile_threshold) + # define the sinks based on the raw stationary distribution (before any transformation) + threshold <- quantile(stat_dist, quantile_threshold) auto_sink_cells <- names(which(stat_dist >= threshold)) + # apply rank transformation if requested + if (rank_transform) { + stat_dist <- ecdf(stat_dist)(stat_dist) + names(stat_dist) <- rownames(P) + } + # ------------------------------------------------------------------------- # return results # ------------------------------------------------------------------------- - if(return_seurat){ + if (return_seurat) { # add to the seurat object metadata seurat_obj <- AddMetaData( object = seurat_obj, @@ -311,7 +323,7 @@ PredictAttractors <- function( col.name = output_name ) return(seurat_obj) - + } else { # return a list containing both the continuous scores and the discrete sink cells return(list( diff --git a/_pkgdown.yml b/_pkgdown.yml index 6db7817..0e8bf35 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -26,12 +26,27 @@ navbar: structure: left: [articles, reference, news] right: [search, github] -components: - news: - text: News - href: news/index.html - articles: - text: Vignettes - menu: - - text: Quick Start - href: articles/quickstart.html \ No newline at end of file + components: + news: + text: News + href: news/index.html + articles: + text: Vignettes + menu: + - text: Quick Start + href: articles/quickstart.html + - text: Basics — Simulated Branching Trajectory + href: articles/simulation_tutorial.html + - text: Advanced Use Cases — NSCLC CD8+ T Cells + href: articles/advanced_NSCLC.html + +articles: +- title: Getting Started + desc: Introductory tutorials covering the core compact workflow. + contents: + - quickstart + - simulation_tutorial +- title: Advanced Use Cases + desc: Advanced workflows using custom gene sets, custom networks, and multi-perturbation analyses. + contents: + - advanced_NSCLC diff --git a/docs/._index.html b/docs/._index.html new file mode 100644 index 0000000..b4dfb61 Binary files /dev/null and b/docs/._index.html differ diff --git a/docs/404.html b/docs/404.html index 3f45b50..4ed1282 100644 --- a/docs/404.html +++ b/docs/404.html @@ -1,8 +1,7 @@ - - + Page not found (404) • compact @@ -20,7 +19,7 @@ compact - 0.0.3 + 0.1.0 -