feat: add ClonalDiversityPlot

NAMESPACE

@@ -4,6 +4,7 @@ export(CellDimPlot)
 export(CellStatPlot)
 export(ClonalAbundancePlot)
 export(ClonalCompositionPlot)
+export(ClonalDiversityPlot)
 export(ClonalLengthPlot)
 export(ClonalOverlapPlot)
 export(ClonalResidencyPlot)
@@ -35,6 +36,7 @@ importFrom(dplyr,n)
 importFrom(dplyr,rename_with)
 importFrom(dplyr,select)
 importFrom(dplyr,slice_min)
+importFrom(dplyr,slice_sample)
 importFrom(dplyr,starts_with)
 importFrom(dplyr,summarise)
 importFrom(dplyr,ungroup)
@@ -86,12 +88,14 @@ importFrom(rlang,sym)
 importFrom(rlang,syms)
 importFrom(scRepertoire,addVariable)
 importFrom(scRepertoire,clonalAbundance)
+importFrom(scRepertoire,clonalDiversity)
 importFrom(scRepertoire,clonalHomeostasis)
 importFrom(scRepertoire,clonalLength)
 importFrom(scRepertoire,clonalOverlap)
 importFrom(scRepertoire,clonalProportion)
 importFrom(scRepertoire,clonalQuant)
 importFrom(scRepertoire,clonalScatter)
+importFrom(stats,as.dist)
 importFrom(stats,cor.test)
 importFrom(stats,quantile)
 importFrom(stringr,str_wrap)

R/clonaldivplot.R

@@ -0,0 +1,172 @@
+#' Calculate the clonal diversities.
+#'
+#' @keywords internal
+#' @importFrom utils getFromNamespace
+#' @importFrom dplyr slice_sample
+ClonalDiversity <- function(
+    input.data, cloneCall = "gene", chain = "both",
+    method = c("shannon", "gini.coeff", "inv.simpson", "norm.entropy", "gini.simpson", "chao1", "ACE"),
+    group_by = NULL, n_boots = 100) {
+    method <- match.arg(method)
+    if (method == "gini.coeff") {
+        div_fn <- function(d) {
+            n <- length(d)
+            1 / n * (n + 1 - 2 * sum((n + 1 - 1:n) * d) / sum(d))
+        }
+    } else {
+        div_fn <- getFromNamespace(paste0(".", method), "scRepertoire")
+    }
+
+    is_seurat_object <- getFromNamespace("is_seurat_object", "scRepertoire")
+    is_se_object <- getFromNamespace("is_se_object", "scRepertoire")
+    .data.wrangle <- getFromNamespace(".data.wrangle", "scRepertoire")
+    .theCall <- getFromNamespace(".theCall", "scRepertoire")
+    .groupList <- getFromNamespace(".groupList", "scRepertoire")
+    .short.check <- getFromNamespace(".short.check", "scRepertoire")
+
+    sco <- is_seurat_object(input.data) | is_se_object(input.data)
+    input.data <- .data.wrangle(input.data, group_by, .theCall(input.data,
+        cloneCall,
+        check.df = FALSE
+    ), chain)
+    cloneCall <- .theCall(input.data, cloneCall)
+    mat <- NULL
+    sample <- c()
+    if (!is.null(group_by) & !sco) {
+        input.data <- .groupList(input.data, group_by)
+    }
+    min <- .short.check(input.data, cloneCall)
+    for (i in seq_along(input.data)) {
+        data <- as.data.frame(table(input.data[[i]][, cloneCall]))
+        mat_a <- NULL
+        sample <- c()
+        if (n_boots > 0) {
+            sample <- div_fn(data$Freq)
+            mat_a <- rbind(mat_a, sample)
+            mat_a[is.na(mat_a)] <- 0
+            mat <- rbind(mat, mat_a)
+            mat <- as.data.frame(mat)
+        } else {
+            for (j in seq(seq_len(n_boots))) {
+                x <- slice_sample(data, n = min)
+                sample <- div_fn(x$Freq)
+                mat_a <- rbind(mat_a, sample)
+            }
+            mat_a[is.na(mat_a)] <- 0
+            mat_b <- colMeans(mat_a)
+            mat_b <- as.data.frame(t(mat_b))
+            mat <- rbind(mat, mat_b)
+        }
+    }
+    if (is.null(group_by)) {
+        group_by <- "Group"
+    }
+    colnames(mat) <- method
+    mat[, group_by] <- names(input.data)
+    mat
+}
+
+#' ClonalDiversityPlot
+#'
+#' Plot the clonal diversities of the samples/groups.
+#'
+#' @param data The product of [scRepertoire::combineTCR], [scRepertoire::combineTCR], or
+#'  [scRepertoire::combineExpression].
+#' @param clone_call How to call the clone - VDJC gene (gene), CDR3 nucleotide (nt),
+#'  CDR3 amino acid (aa), VDJC gene + CDR3 nucleotide (strict) or a custom variable
+#'  in the data
+#' @param chain indicate if both or a specific chain should be used - e.g. "both",
+#'  "TRA", "TRG", "IGH", "IGL"
+#' @param method The method to calculate the diversity. Options are "shannon" (default),
+#'  "inv.simpson", "norm.entropy", "gini.simpson", "chao1" and "ACE".
+#'  See [scRepertoire::clonalDiversity] for details.
+#' @param plot_type The type of plot. Options are "bar", "box" and "violin".
+#' @param position The position adjustment for the bars. Default is "dodge".
+#' @param group_by A character vector of column names to group the samples. Default is NULL.
+#' @param facet_by A character vector of column names to facet the plots. Default is NULL.
+#' @param split_by A character vector of column names to split the plots. Default is NULL.
+#' @param xlab The x-axis label. Default is NULL.
+#' @param ylab The y-axis label. Default is NULL.
+#' @param ... Other arguments passed to the specific plot function.
+#'  * For "bar", [plotthis::BarPlot]
+#'  * For "box", [plotthis::BoxPlot]
+#'  * For "violin", [plotthis::ViolinPlot]
+#' @return A ggplot object or a list if `combine` is FALSE
+#' @export
+#' @importFrom tidyr separate
+#' @importFrom scRepertoire clonalDiversity
+#' @importFrom plotthis BarPlot BoxPlot ViolinPlot
+#' @examples
+#' set.seed(8525)
+#' data(contig_list, package = "scRepertoire")
+#' data <- scRepertoire::combineTCR(contig_list,
+#'     samples = c("P17B", "P17L", "P18B", "P18L", "P19B","P19L", "P20B", "P20L"))
+#' data <- scRepertoire::addVariable(data,
+#'     variable.name = "Type",
+#'     variables = rep(c("B", "L"), 4)
+#' )
+#' data <- scRepertoire::addVariable(data,
+#'     variable.name = "Subject",
+#'     variables = rep(c("P17", "P18", "P19", "P20"), each = 2)
+#' )
+#'
+#' ClonalDiversityPlot(data)
+#' ClonalDiversityPlot(data, group_by = "Type")
+#' ClonalDiversityPlot(data, group_by = "Type", plot_type = "box")
+#' ClonalDiversityPlot(data, group_by = "Type", plot_type = "violin")
+#' ClonalDiversityPlot(data, group_by = "Type", plot_type = "violin",
+#'   method = "gini.coeff", add_box = TRUE)
+ClonalDiversityPlot <- function(
+    data, clone_call = "gene", chain = "both",
+    method = c("shannon", "gini.coeff", "inv.simpson", "norm.entropy", "gini.simpson", "chao1", "ACE"),
+    plot_type = c("bar", "box", "violin"), position = "dodge",
+    group_by = NULL, facet_by = NULL, split_by = NULL,
+    xlab = NULL, ylab = NULL,
+    ...) {
+    method <- match.arg(method)
+    plot_type <- match.arg(plot_type)
+
+    if (plot_type %in% c("box", "violin") && is.null(group_by)) {
+        stop("'group_by' must be provided for box/violin ClonalDiversityPlot")
+    }
+    all_groupings <- unique(c("Sample", group_by, facet_by, split_by))
+    method_name <- switch(method,
+        shannon = "Shannon Index",
+        gini.coeff = "Gini Coefficient",
+        inv.simpson = "Inverse Simpson Index",
+        norm.entropy = "Normalized Entropy",
+        gini.simpson = "Gini-Simpson Index",
+        chao1 = "Chao1 Index",
+        ACE = "ACE Index"
+    )
+
+    data <- MergeClonalGroupings(data, all_groupings)
+    data <- ClonalDiversity(data,
+        cloneCall = clone_call, chain = chain, method = method,
+        group_by = ".group"
+    )
+    data <- separate(data, ".group", into = all_groupings, sep = " // ")
+
+    if (plot_type == "bar") {
+        x <- group_by %||% "Sample"
+        group_by <- if(is.null(group_by)) NULL else "Sample"
+        BarPlot(data,
+            x = x, y = method, group_by = group_by, position = position,
+            xlab = xlab %||% group_by, ylab = ylab %||% method_name,
+            split_by = split_by, facet_by = facet_by, ...
+        )
+    } else if (plot_type == "box") {
+        BoxPlot(data,
+            x = group_by, y = method,
+            xlab = xlab %||% group_by, ylab = ylab %||% method_name,
+            split_by = split_by, facet_by = facet_by, ...
+        )
+    } else if (plot_type == "violin") {
+        ViolinPlot(data,
+            x = group_by, y = method,
+            xlab = xlab %||% group_by, ylab = ylab %||% method_name,
+            split_by = split_by, facet_by = facet_by, ...
+        )
+    }
+}
+

R/clonalstatplot.R

@@ -1076,6 +1076,7 @@ ClonalCompositionPlot <- function(
 #' @param show_column_names Whether to show the column names. Default is TRUE.
 #' @param ... Other arguments passed to the specific plot function [plotthis::Heatmap].
 #' @return A ComplexHeatmap object or a list if `combine` is FALSE
+#' @importFrom stats as.dist
 #' @importFrom rlang syms
 #' @importFrom dplyr %>% filter
 #' @importFrom tidyr separate pivot_longer unite

R/utils.R

@@ -5,3 +5,15 @@ check_columns <- getFromNamespace("check_columns", "plotthis")
 #' @keywords internal
 #' @importFrom utils getFromNamespace
 combine_plots <- getFromNamespace("combine_plots", "plotthis")
+
+# #' Monkey patch a function from a namespace
+# #' @keywords internal
+# #' @param namespace The namespace where the function is located
+# #' @param function_name The name of the function to be patched
+# #' @param new_function The new function to be patched
+# monkey_patch <- function(namespace, function_name, new_function) {
+#     ns <- getNamespace(namespace)
+#     unlockBinding(function_name, ns)
+#     assign(function_name, new_function, envir = ns)
+#     lockBinding(function_name, ns)
+# }

_pkgdown.yml

@@ -20,6 +20,7 @@ reference:
     - ClonalResidencyPlot
     - ClonalCompositionPlot
     - ClonalOverlapPlot
+    - ClonalDiversityPlot
   - title: data
     desc: Built-in data sets
     contents: