From c2ef10e591561e27c620e5b5b06ca2b0b65e9464 Mon Sep 17 00:00:00 2001 From: Anthony Christidis Date: Wed, 18 Sep 2024 17:10:19 -0400 Subject: [PATCH] Add new SIR space functionality. --- DESCRIPTION | 4 +- NAMESPACE | 4 + R/calculateSIRSpace.R | 99 +++++++ R/plot.calculateSIRSpace.R | 253 ++++++++++++++++++ R/projectPCA.R | 64 ++--- R/projectSIR.R | 232 ++++++++++++++++ .../VisualizationTools/calculateSIRSpace1.png | Bin 0 -> 21868 bytes .../VisualizationTools/calculateSIRSpace2.png | Bin 0 -> 27308 bytes man/calculateSIRSpace.Rd | 104 +++++++ man/conditionalMeans.Rd | 56 ++++ man/projectPCA.Rd | 24 +- man/projectSIR.Rd | 76 ++++++ tests/testthat/test-calculateSIRSpace.R | 89 ++++++ tests/testthat/test-plot.calculateSIRSpace.R | 29 ++ tests/testthat/test-projectSIR.R | 47 ++++ tests/testthat/test-regressPC.R | 28 +- vignettes/VisualizationTools.Rmd | 39 ++- 17 files changed, 1088 insertions(+), 60 deletions(-) create mode 100644 R/calculateSIRSpace.R create mode 100644 R/plot.calculateSIRSpace.R create mode 100644 R/projectSIR.R create mode 100644 inst/extdata/compressed/VisualizationTools/calculateSIRSpace1.png create mode 100644 inst/extdata/compressed/VisualizationTools/calculateSIRSpace2.png create mode 100644 man/calculateSIRSpace.Rd create mode 100644 man/conditionalMeans.Rd create mode 100644 man/projectSIR.Rd create mode 100644 tests/testthat/test-calculateSIRSpace.R create mode 100644 tests/testthat/test-plot.calculateSIRSpace.R create mode 100644 tests/testthat/test-projectSIR.R diff --git a/DESCRIPTION b/DESCRIPTION index effe8cb..dd3a274 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -36,7 +36,9 @@ Imports: transport, speedglm, cramer, - rlang + rlang, + bluster, + patchwork Suggests: AUCell, BiocStyle, diff --git a/NAMESPACE b/NAMESPACE index 9aae822..847f99a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,6 +5,7 @@ S3method(plot,calculateCellDistances) S3method(plot,calculateCellSimilarityPCA) S3method(plot,calculateDiscriminantSpace) S3method(plot,calculateNearestNeighborProbabilities) +S3method(plot,calculateSIRSpace) S3method(plot,compareCCA) S3method(plot,comparePCA) S3method(plot,comparePCASubspace) @@ -21,6 +22,7 @@ export(calculateDiscriminantSpace) export(calculateHVGOverlap) export(calculateHotellingPValue) export(calculateNearestNeighborProbabilities) +export(calculateSIRSpace) export(calculateVarImpOverlap) export(compareCCA) export(comparePCA) @@ -36,6 +38,7 @@ export(plotPairwiseDistancesDensity) export(plotQCvsAnnotation) export(plotWassersteinDistance) export(projectPCA) +export(projectSIR) export(regressPC) import(SingleCellExperiment) importFrom(SummarizedExperiment,assay) @@ -59,4 +62,5 @@ importFrom(stats,quantile) importFrom(stats,sd) importFrom(stats,setNames) importFrom(utils,combn) +importFrom(utils,head) importFrom(utils,tail) diff --git a/R/calculateSIRSpace.R b/R/calculateSIRSpace.R new file mode 100644 index 0000000..ee0f858 --- /dev/null +++ b/R/calculateSIRSpace.R @@ -0,0 +1,99 @@ +#' @title Calculate Sliced Inverse Regression (SIR) Space for Different Cell Types +#' +#' @description +#' This function calculates the SIR space projections for different cell types in the query and reference datasets. +#' +#' @details +#' The function projects the query dataset onto the SIR space of the reference dataset based on shared cell types. +#' It computes conditional means for the reference dataset, extracts the SVD components, and performs the projection +#' of both the query and reference data. It uses the `projectSIR` function to perform the actual projection and +#' allows the user to specify particular cell types for analysis. +#' +#' @param query_data A \code{\linkS4class{SingleCellExperiment}} object containing the numeric expression matrix for the query cells. +#' @param reference_data A \code{\linkS4class{SingleCellExperiment}} object containing the numeric expression matrix for the reference cells. +#' @param query_cell_type_col A character string specifying the column name in the \code{colData} of \code{query_data} that identifies the cell types. +#' @param ref_cell_type_col A character string specifying the column name in the \code{colData} of \code{reference_data} that identifies the cell types. +#' @param cell_types A character vector specifying the cell types to include in the analysis. If NULL, all common cell types between the query and reference data will be used. +#' @param multiple_cond_means Logical. Whether to compute conditional means for multiple conditions in the reference dataset. Default is TRUE. +#' @param assay_name A character string specifying the name of the assay on which to perform computations. Default is "logcounts". +#' @param cumulative_variance_threshold A numeric value specifying the cumulative variance threshold for selecting principal components. Default is 0.7. +#' @param n_neighbor A numeric value specifying the number of neighbors for computing the SIR space. Default is 1. +#' +#' @return A list containing the SIR projections, rotation matrix, and percentage of variance explained for the given cell types. +#' +#' @export +#' +#' @author Anthony Christidis, \email{anthony-alexander_christidis@hms.harvard.edu} +#' +#' @seealso \code{\link{plot.calculateSIRSpace}} +#' +#' @examples +#' # Load data +#' data("reference_data") +#' data("query_data") +#' +#' # Compute important variables for all pairwise cell comparisons +#' sir_output <- calculateSIRSpace(reference_data = reference_data, +#' query_data = query_data, +#' query_cell_type_col = "SingleR_annotation", +#' ref_cell_type_col = "expert_annotation", +#' multiple_cond_means = TRUE, +#' cumulative_variance_threshold = 0.9, +#' n_neighbor = 1) +#' +#' # Generate boxplot of SIR projections +#' plot(sir_output, plot_type = "boxplot", sir_subset = 1:6) +#' +# Function to plot cell types in SIR space +calculateSIRSpace <- function(query_data, + reference_data, + query_cell_type_col, + ref_cell_type_col, + cell_types = NULL, + multiple_cond_means = TRUE, + assay_name = "logcounts", + cumulative_variance_threshold = 0.7, + n_neighbor = 1){ + + # Check standard input arguments + argumentCheck(query_data = query_data, + reference_data = reference_data, + query_cell_type_col = query_cell_type_col, + ref_cell_type_col = ref_cell_type_col, + cell_types = cell_types, + assay_name = assay_name) + + # Check if cumulative_variance_threshold is between 0 and 1 + if (!is.numeric(cumulative_variance_threshold) || + cumulative_variance_threshold < 0 || cumulative_variance_threshold > 1) { + stop("cumulative_variance_threshold must be a numeric value between 0 and 1.") + } + + # Check if n_neighbor is a positive integer + if (!is.numeric(n_neighbor) || n_neighbor <= 0 || n_neighbor != as.integer(n_neighbor)) { + stop("n_neighbor must be a positive integer.") + } + + # Get common cell types if they are not specified by user + if(is.null(cell_types)){ + cell_types <- na.omit(unique(c(reference_data[[ref_cell_type_col]], + query_data[[query_cell_type_col]]))) + } + + # Get the projected PCA data + sir_output <- projectSIR(query_data = query_data, + reference_data = reference_data, + query_cell_type_col = query_cell_type_col, + ref_cell_type_col = ref_cell_type_col, + cell_types = cell_types, + multiple_cond_means = multiple_cond_means, + assay_name = assay_name, + cumulative_variance_threshold = cumulative_variance_threshold, + n_neighbor = n_neighbor) + + # Return SIR projections output + class(sir_output) <- c(class(sir_output), "calculateSIRSpace") + return(sir_output) +} + + diff --git a/R/plot.calculateSIRSpace.R b/R/plot.calculateSIRSpace.R new file mode 100644 index 0000000..65ac3fb --- /dev/null +++ b/R/plot.calculateSIRSpace.R @@ -0,0 +1,253 @@ +#' @title Plot Projected Data on Discriminant Spaces +#' +#' @description +#' The S3 plot method visualizes the projected reference and query data on discriminant spaces using either a scatterplot, boxplot, or varplot. +#' +#' @details +#' - **Scatterplot**: Displays projected data points, with colors used to differentiate between cell types and datasets. +#' - **Boxplot**: Shows the distribution of projected data values for each cell type, separated by datasets. +#' - **Varplot**: Highlights the top contributing variables for each discriminant axis, differentiating between positive and negative loadings. +#' +#' @param x An object of class \code{calculateSIRSpace}, which contains projected data on the discriminant space. +#' Each element should include 'ref_proj' and 'query_proj' data frames representing reference and query projections. +#' @param plot_type A character string indicating the type of plot to generate. Options are "scatterplot", "boxplot", or "varplot". Default is "scatterplot". +#' @param sir_subset A numeric vector specifying which discriminant axes (SIR components) to include in the plot. Default is the first 5 axes. +#' @param n_top_vars Number of top contributing variables to display in varplot. Default is 5. +#' @param ... Additional arguments to be passed to the plotting functions. +#' +#' @keywords internal +#' +#' @return A \code{ggplot} object representing the chosen visualization (scatterplot, boxplot, or varplot) of the projected data. +#' +#' @export +#' +#' @author Anthony Christidis, \email{anthony-alexander_christidis@hms.harvard.edu} +#' +#' @seealso \code{\link{calculateSIRSpace}} +#' +#' @rdname calculateSIRSpace +#' +#' @importFrom utils head +#' +# This function plots the SIR output either as a scatterplot for a boxplot +plot.calculateSIRSpace <- function(x, + plot_type = c("scatterplot", "boxplot", "varplot"), + sir_subset = NULL, + n_top_vars = 5, + ...){ + + # Match argument for plot type + plot_type <- match.arg(plot_type) + + # Value for SIR subset + if(is.null(sir_subset)){ + sir_subset <- seq_len(min(5, ncol(x[["rotation_mat"]]))) + } + + # Check input for SIR subset + if(any(!(sir_subset %in% seq_len(ncol(x[["rotation_mat"]]))))){ + stop("\"sir_subset\" is out of range.") + } + + # Identify cell types + cell_types <- unique(rownames(x[["cond_means"]])) + + if(plot_type == "scatterplot"){ + + # Remove cells with NA cell type + sir_projections <- na.omit(x[["sir_projections"]]) + + # Create a subset of plot names for the SIR subset + plot_names <- paste0( + "SIR", sir_subset, " (", + sprintf("%.1f%%", x[["percent_var"]][sir_subset]), ")") + + # Create all possible pairs of the specified PCs in the subset + pairs <- expand.grid(x = sir_subset, y = sir_subset) + pairs <- pairs[pairs[["x"]] < pairs[["y"]], ] + + # Create a new data frame with all possible pairs of specified PCs + data_pairs_list <- lapply(seq_len(nrow(pairs)), function(i) { + x_col <- pairs[["x"]][i] + y_col <- pairs[["y"]][i] + + # Properly subset the projection data using sir_subset indices + data_frame <- data.frame( + sir_projections[, c(x_col, y_col)], + paste(sir_projections[["dataset"]], sir_projections[["cell_type"]], sep = " ")) + + colnames(data_frame) <- c("x_value", "y_value", "cell_type_dataset") + data_frame[["x"]] <- plot_names[match(x_col, sir_subset)] + data_frame[["y"]] <- plot_names[match(y_col, sir_subset)] + + return(data_frame) + }) + + # Combine all pairs into one data frame + data_pairs <- do.call(rbind, data_pairs_list) + data_pairs[["x"]] <- factor(data_pairs[["x"]], + levels = plot_names) + data_pairs[["y"]] <- factor(data_pairs[["y"]], + levels = plot_names) + + # Define the order of cell type and dataset combinations + order_combinations <- paste(rep(c("Reference", "Query"), length(cell_types)), rep(sort(cell_types), each = 2)) + data_pairs[["cell_type_dataset"]] <- factor(data_pairs[["cell_type_dataset"]], levels = order_combinations) + cell_type_colors <- generateColors(order_combinations, paired = TRUE) + + # Set SIR identity as factor + data_pairs[["x"]] <- factor(data_pairs[["x"]], levels = unique(data_pairs[["x"]])) + data_pairs[["y"]] <- factor(data_pairs[["y"]], levels = unique(data_pairs[["y"]])) + + # Create the ggplot object (with facets if PCA) + plot_obj <- ggplot2::ggplot( + data_pairs, ggplot2::aes(x = .data[["x_value"]], + y = .data[["y_value"]], + color = .data[["cell_type_dataset"]])) + + ggplot2::geom_point(alpha = 0.5, size = 1) + + ggplot2::scale_color_manual(values = cell_type_colors, + name = "Cell Types") + + ggplot2::facet_grid(rows = ggplot2::vars(.data[["y"]]), + cols = ggplot2::vars(.data[["x"]]), + scales = "free") + + ggplot2::xlab("") + ggplot2::ylab("") + + ggplot2::theme_bw() + + ggplot2::theme( + panel.grid.minor = ggplot2::element_blank(), + panel.grid.major = ggplot2::element_line(color = "gray", + linetype = "dotted"), + plot.title = ggplot2::element_text(size = 14, + face = "bold", + hjust = 0.5), + axis.title = ggplot2::element_text(size = 12), + axis.text = ggplot2::element_text(size = 10)) + + } else if (plot_type == "boxplot"){ + + # Remove cells with NA cell type + sir_projections <- na.omit(x[["sir_projections"]]) + + # Create the long format data frame manually + sir_projections <- sir_projections[!is.na(sir_projections[["cell_type"]]),] + if(!is.null(cell_types)){ + if(all(cell_types %in% sir_projections[["cell_type"]])){ + sir_projections <- sir_projections[which(sir_projections[["cell_type"]] %in% + cell_types),] + } else{ + stop("One or more of the specified \'cell_types\' are not available.") + } + } + sir_long <- data.frame(SIR = rep(paste0("sir", sir_subset), + each = nrow(sir_projections)), + Value = unlist(c(sir_projections[, sir_subset])), + dataset = rep(sir_projections[["dataset"]], + length(sir_subset)), + cell_type = rep(sir_projections[["cell_type"]], + length(sir_subset))) + sir_long[["SIR"]] <- toupper(sir_long[["SIR"]]) + + # Create a new variable representing the combination of cell type and dataset + sir_long[["cell_type_dataset"]] <- paste(sir_long[["dataset"]], + sir_long[["cell_type"]], + sep = " ") + + # Define the order of cell type and dataset combinations + order_combinations <- paste(rep(c("Reference", "Query"), + length(unique(sir_long[["cell_type"]]))), + rep(sort(unique(sir_long[["cell_type"]])), + each = 2)) + + # Reorder the levels of cell type and dataset factor + sir_long[["cell_type_dataset"]] <- factor(sir_long[["cell_type_dataset"]], + levels = order_combinations) + + # Set SIR identity as factor + sir_long[["SIR"]] <- factor(sir_long[["SIR"]], levels = paste0("SIR", sir_subset)) + + # Define the colors for cell types + cell_type_colors <- generateColors(order_combinations, + paired = TRUE) + + # Create the ggplot + plot_obj <- ggplot2::ggplot(sir_long, ggplot2::aes( + x = .data[["cell_type"]], + y = .data[["Value"]], + fill = .data[["cell_type_dataset"]])) + + ggplot2::geom_boxplot(alpha = 0.7, outlier.shape = NA, width = 0.7) + + ggplot2::facet_wrap(~ .data[["SIR"]], scales = "free") + + ggplot2::scale_fill_manual(values = cell_type_colors, + name = "Cell Types") + + ggplot2::labs(x = "", y = "SIR Score") + + ggplot2::theme_bw() + + ggplot2::theme( + panel.grid.minor = ggplot2::element_blank(), + panel.grid.major = ggplot2::element_line(color = "gray", linetype = "dotted"), + plot.title = ggplot2::element_text(size = 14, face = "bold", hjust = 0.5), + axis.title = ggplot2::element_text(size = 12), axis.text = ggplot2::element_text(size = 10), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, size = 10)) + + } else if (plot_type == "varplot"){ + + # Rotation matrix + rotation_mat <- x[["rotation_mat"]][, sir_subset, drop = FALSE] + + # Create an empty list to store individual ggplot objects + plot_list <- list() + + # For each SIR vector (each column of rotation_mat) + for (i in seq_len(ncol(rotation_mat))) { + + sir_vector <- rotation_mat[, i] + + # Ensure that the vector has names + if (is.null(names(sir_vector))) { + names(sir_vector) <- rownames(rotation_mat) + } + + # Find the top 5 positive and top 5 negative markers + sorted_positive <- sort(sir_vector[sir_vector > 0], decreasing = TRUE) + sorted_negative <- sort(sir_vector[sir_vector < 0], decreasing = FALSE) + + top_positive <- head(sorted_positive, n_top_vars) + top_negative <- head(sorted_negative, n_top_vars) + + # Create a data frame for this SIR vector with marker names and loadings + df <- data.frame( + marker = c(names(top_positive), names(top_negative)), + loading = c(top_positive, top_negative), + SIR = rep(paste0("SIR", i), length(c(top_positive, top_negative))), + Direction = c(rep("Positive", length(top_positive)), rep("Negative", length(top_negative))) + ) + + # Ensure that 'marker' is treated as a factor, preserving the original order + df[["marker"]] <- factor(df[["marker"]], levels = df[["marker"]]) + + # Set SIR identity as factor + df[["SIR"]] <- factor(df[["SIR"]], levels = paste0("SIR", sir_subset)) + + # Create a plot for this SIR vector + plot_list[[i]] <- ggplot2::ggplot(df, ggplot2::aes( + x = .data[["marker"]], + y = .data[["loading"]], + fill = .data[["Direction"]])) + + ggplot2::geom_bar(stat = "identity") + + ggplot2::facet_wrap(~ .data[["SIR"]], scales = "free_y") + + ggplot2::coord_flip() + + ggplot2::labs(x = NULL, y = "SIR Loading", title = NULL) + + ggplot2::theme_bw() + + ggplot2::theme( + legend.position = "none", + panel.grid.minor = ggplot2::element_blank(), + panel.grid.major = ggplot2::element_line(color = "gray", linetype = "dotted"), + plot.title = ggplot2::element_text(size = 14, face = "bold", hjust = 0.5), + axis.title = ggplot2::element_text(size = 12), axis.text = ggplot2::element_text(size = 10), + axis.text.x = ggplot2::element_text(hjust = 1, size = 10)) + } + + # Combine all plots into a grid + plot_obj <- patchwork::wrap_plots(plot_list, ncol = 3) + } + + # Return plot object + return(plot_obj) +} diff --git a/R/projectPCA.R b/R/projectPCA.R index 218fb8d..c72c991 100644 --- a/R/projectPCA.R +++ b/R/projectPCA.R @@ -1,23 +1,23 @@ #' @title Project Query Data Onto PCA Space of Reference Data #' -#' @description -#' This function projects a query singleCellExperiment object onto the PCA space of a reference -#' singleCellExperiment object. The PCA analysis on the reference data is assumed to be pre-computed +#' @description +#' This function projects a query singleCellExperiment object onto the PCA space of a reference +#' singleCellExperiment object. The PCA analysis on the reference data is assumed to be pre-computed #' and stored within the object. #' -#' @details -#' This function assumes that the "PCA" element exists within the \code{reducedDims} of the reference data -#' (obtained using \code{reducedDim(reference_data)}) and that the genes used for PCA are present in both -#' the reference and query data. It performs centering and scaling of the query data based on the reference +#' @details +#' This function assumes that the "PCA" element exists within the \code{reducedDims} of the reference data +#' (obtained using \code{reducedDim(reference_data)}) and that the genes used for PCA are present in both +#' the reference and query data. It performs centering and scaling of the query data based on the reference #' data before projection. #' -#' @param query_data A \code{\linkS4class{SingleCellExperiment}} object containing numeric expression matrix +#' @param query_data A \code{\linkS4class{SingleCellExperiment}} object containing numeric expression matrix #' for the query cells. -#' @param reference_data A \code{\linkS4class{SingleCellExperiment}} object containing numeric expression matrix +#' @param reference_data A \code{\linkS4class{SingleCellExperiment}} object containing numeric expression matrix #' for the reference cells. -#' @param query_cell_type_col character. The column name in the \code{colData} of \code{query_data} +#' @param query_cell_type_col character. The column name in the \code{colData} of \code{query_data} #' that identifies the cell types. -#' @param ref_cell_type_col character. The column name in the \code{colData} of \code{reference_data} +#' @param ref_cell_type_col character. The column name in the \code{colData} of \code{reference_data} #' that identifies the cell types. #' @param pc_subset A numeric vector specifying the subset of principal components (PCs) to compare. Default is 1:10. #' @param assay_name Name of the assay on which to perform computations. @@ -32,22 +32,22 @@ #' # Load data #' data("reference_data") #' data("query_data") -#' +#' #' # Project the query data onto PCA space of reference -#' pca_output <- projectPCA(query_data = query_data, +#' pca_output <- projectPCA(query_data = query_data, #' reference_data = reference_data, -#' query_cell_type_col = "SingleR_annotation", -#' ref_cell_type_col = "expert_annotation", +#' query_cell_type_col = "SingleR_annotation", +#' ref_cell_type_col = "expert_annotation", #' pc_subset = 1:10) #' # Function to project query data onto PCA space of reference data -projectPCA <- function(query_data, - reference_data, - query_cell_type_col, - ref_cell_type_col, +projectPCA <- function(query_data, + reference_data, + query_cell_type_col, + ref_cell_type_col, pc_subset = 1:10, assay_name = "logcounts"){ - + # Check standard input arguments argumentCheck(query_data = query_data, reference_data = reference_data, @@ -61,29 +61,29 @@ projectPCA <- function(query_data, rotation_mat <- attributes(reducedDim( reference_data, "PCA"))[["rotation"]][, pc_subset, drop = FALSE] PCA_genes <- rownames(rotation_mat) - + # Check if genes used for PCA are available in query data if (!all(PCA_genes %in% rownames(assay(query_data)))) { stop("Genes in reference PCA are not found in query data.") } - + # Center and scale query data based on reference for projection centering_vec <- apply(t(as.matrix( assay(reference_data, assay_name))), 2, mean)[PCA_genes] query_mat <- scale(t(as.matrix( - assay(query_data, assay_name)))[, PCA_genes, drop = FALSE], + assay(query_data, assay_name)))[, PCA_genes, drop = FALSE], center = centering_vec, scale = FALSE) %*% rotation_mat - + # Returning output as a dataframe - return(data.frame(rbind(ref_mat, query_mat), - dataset = c(rep("Reference", nrow(ref_mat)), + return(data.frame(rbind(ref_mat, query_mat), + dataset = c(rep("Reference", nrow(ref_mat)), rep("Query", nrow(query_mat))), - cell_type = c(ifelse(rep(is.null(ref_cell_type_col), - nrow(ref_mat)), - rep(NA, nrow(ref_mat)), + cell_type = c(ifelse(rep(is.null(ref_cell_type_col), + nrow(ref_mat)), + rep(NA, nrow(ref_mat)), reference_data[[ref_cell_type_col]]), - ifelse(rep(is.null(query_cell_type_col), - nrow(query_mat)), - rep(NA, nrow(query_mat)), + ifelse(rep(is.null(query_cell_type_col), + nrow(query_mat)), + rep(NA, nrow(query_mat)), query_data[[query_cell_type_col]])))) } diff --git a/R/projectSIR.R b/R/projectSIR.R new file mode 100644 index 0000000..8910ef7 --- /dev/null +++ b/R/projectSIR.R @@ -0,0 +1,232 @@ +#' @title Project Query Data Onto SIR Space of Reference Data +#' +#' @description +#' This function projects a query \code{SingleCellExperiment} object onto the SIR (supervised independent +#' component) space of a reference \code{SingleCellExperiment} object. The SVD of the reference data is +#' computed on conditional means per cell type, and the query data is projected based on these reference +#' components. +#' +#' @details +#' The genes used for the projection (SVD) must be present in both the reference and query datasets. +#' The function first computes conditional means for each cell type in the reference data, then performs +#' SVD on these conditional means to obtain the rotation matrix used for projecting both the reference +#' and query datasets. The query data is centered and scaled based on the reference data. +#' +#' @param query_data A \code{\linkS4class{SingleCellExperiment}} object containing numeric expression matrix for the query cells. +#' @param reference_data A \code{\linkS4class{SingleCellExperiment}} object containing numeric expression matrix for the reference cells. +#' @param query_cell_type_col A character string specifying the column in the \code{colData} of \code{query_data} +#' that identifies the cell types. +#' @param ref_cell_type_col A character string specifying the column in the \code{colData} of \code{reference_data} +#' that identifies the cell types. +#' @param cell_types A character vector of cell types for which to compute conditional means in the reference data. +#' @param multiple_cond_means A logical value indicating whether to compute multiple conditional means per cell type +#' (through PCA and clustering). Defaults to \code{TRUE}. +#' @param assay_name A character string specifying the assay name on which to perform computations. Defaults to \code{"logcounts"}. +#' @param cumulative_variance_threshold A numeric value between 0 and 1 specifying the variance threshold for PCA +#' when computing multiple conditional means. Defaults to \code{0.7}. +#' @param n_neighbor An integer specifying the number of nearest neighbors for clustering when computing multiple +#' conditional means. Defaults to \code{1}. +#' +#' @return A list containing: +#' \item{cond_means}{A matrix of the conditional means computed for the reference data.} +#' \item{rotation_mat}{The rotation matrix obtained from the SVD of the conditional means.} +#' \item{sir_projections}{A \code{data.frame} containing the SIR projections for both the reference and query datasets.} +#' \item{percent_var}{The percentage of variance explained by each component of the SIR projection.} +#' +#' @export +#' +#' @author +#' Anthony Christidis, \email{anthony-alexander_christidis@hms.harvard.edu} +#' +#' @examples +#' # Load data +#' data("reference_data") +#' data("query_data") +#' +#' # Project the query data onto SIR space of reference +#' sir_output <- projectSIR(query_data = query_data, +#' reference_data = reference_data, +#' query_cell_type_col = "SingleR_annotation", +#' ref_cell_type_col = "expert_annotation") +#' +# Function to project query data onto PCA space of reference data +projectSIR <- function(query_data, + reference_data, + query_cell_type_col, + ref_cell_type_col, + cell_types = NULL, + multiple_cond_means = TRUE, + assay_name = "logcounts", + cumulative_variance_threshold = 0.7, + n_neighbor = 1){ + + # Check standard input arguments + argumentCheck(query_data = query_data, + reference_data = reference_data, + query_cell_type_col = query_cell_type_col, + ref_cell_type_col = ref_cell_type_col, + cell_types = cell_types, + assay_name = assay_name) + + # Check if cumulative_variance_threshold is between 0 and 1 + if (!is.numeric(cumulative_variance_threshold) || + cumulative_variance_threshold < 0 || cumulative_variance_threshold > 1) { + stop("cumulative_variance_threshold must be a numeric value between 0 and 1.") + } + + # Check if n_neighbor is a positive integer + if (!is.numeric(n_neighbor) || n_neighbor <= 0 || n_neighbor != as.integer(n_neighbor)) { + stop("n_neighbor must be a positive integer.") + } + + # Get common cell types if they are not specified by user + if(is.null(cell_types)){ + cell_types <- na.omit(unique(c(reference_data[[ref_cell_type_col]], + query_data[[query_cell_type_col]]))) + } + + # Compute conditional means for each cell type of reference data + cond_means <- conditionalMeans(reference_data = reference_data, + ref_cell_type_col = ref_cell_type_col, + cell_types = cell_types, + multiple_cond_means = multiple_cond_means, + assay_name = assay_name, + cumulative_variance_threshold = cumulative_variance_threshold, + n_neighbor = n_neighbor) + + # Extract reference SVD components and rotation matrix + SVD_genes <- rownames(reference_data) + centering_vec <- apply(cond_means, 2, mean)[SVD_genes] + svd_ref <- svd(scale(cond_means, center = centering_vec, scale = FALSE)) + rotation_mat <- svd_ref[["v"]] + rownames(rotation_mat) <- SVD_genes + + # Check if genes used for PCA are available in query data + if (!all(SVD_genes %in% rownames(assay(query_data, assay_name)))) { + stop("Genes in reference SVD are not found in query data.") + } + + # Compute percentage of variance explained + percent_var <- svd_ref[["d"]]^2 / + sum(svd_ref[["d"]]^2) * 100 + + # Center and scale reference query data based on reference for projection + ref_mat <- scale(t(as.matrix( + assay(reference_data, assay_name)))[, SVD_genes, drop = FALSE], + center = centering_vec, scale = FALSE) %*% rotation_mat + query_mat <- scale(t(as.matrix( + assay(query_data, assay_name)))[, SVD_genes, drop = FALSE], + center = centering_vec, scale = FALSE) %*% rotation_mat + sir_projections <- data.frame( + rbind(ref_mat, query_mat), + dataset = c(rep("Reference", nrow(ref_mat)), + rep("Query", nrow(query_mat))), + cell_type = c(ifelse(rep(is.null(ref_cell_type_col), + nrow(ref_mat)), + rep(NA, nrow(ref_mat)), + reference_data[[ref_cell_type_col]]), + ifelse(rep(is.null(query_cell_type_col), + nrow(query_mat)), + rep(NA, nrow(query_mat)), + query_data[[query_cell_type_col]]))) + colnames(sir_projections)[seq_len(ncol(ref_mat))] <- + paste0("SIR", seq_len(ncol(ref_mat)), " (", + sprintf("%.1f%%", percent_var[seq_len(ncol(ref_mat))]), ")") + + # Returning output as a data frame + return(list(cond_means = cond_means, + rotation_mat = rotation_mat, + sir_projections = sir_projections, + percent_var = percent_var)) +} + +#' @title Compute Conditional Means for Cell Types +#' +#' @description +#' This function computes conditional means for each cell type in the reference data. It can compute +#' either a single conditional mean per cell type or multiple conditional means, depending on the +#' specified settings. Principal component analysis (PCA) is used for dimensionality reduction before +#' clustering when computing multiple conditional means. +#' +#' @details +#' The function offers two modes of operation: +#' - **Single conditional mean per cell type**: For each cell type, it computes the mean expression +#' across all observations. +#' - **Multiple conditional means per cell type**: For each cell type, the function performs PCA to +#' reduce dimensionality, followed by clustering to compute multiple conditional means. +#' +#' @param reference_data A \code{SingleCellExperiment} object containing the reference data, where rows +#' represent genes and columns represent cells. +#' @param ref_cell_type_col A character string specifying the column in \code{colData(reference_data)} that contains the cell type labels. +#' @param cell_types A character vector of cell types for which to compute conditional means. +#' @param multiple_cond_means A logical value indicating whether to compute multiple conditional means per cell type. +#' Defaults to \code{FALSE}. +#' @param assay_name A character string specifying the name of the assay to use for the computation. Defaults to \code{"logcounts"}. +#' @param cumulative_variance_threshold A numeric value between 0 and 1 specifying the variance threshold +#' for PCA when computing multiple conditional means. Defaults to \code{0.7}. +#' @param n_neighbor An integer specifying the number of nearest neighbors for clustering when computing multiple conditional means. +#' Defaults to \code{1}. +#' +#' @return A numeric matrix with the conditional means for each cell type. If \code{multiple_cond_means = TRUE}, the matrix +#' will contain multiple rows for each cell type, representing the different conditional means computed via clustering. +#' +#' @keywords internal +#' +#' @author +#' Anthony Christidis, \email{anthony-alexander_christidis@hms.harvard.edu} +#' +# Function to compute Ledoit-Wolf covariance matrix +conditionalMeans <- function(reference_data, + ref_cell_type_col, + cell_types, + multiple_cond_means = FALSE, + assay_name = "logcounts", + cumulative_variance_threshold = 0.7, + n_neighbor = 1) { + + # Compute conditional means for each cell type of reference data + if (multiple_cond_means) { + + # Matrix to store results + cond_means <- matrix(nrow = 0, ncol = nrow(reference_data)) + colnames(cond_means) <- rownames(reference_data) + + for(cell_type in cell_types){ + + # Compute multiple conditional means per cell type + assay_mat <- scale(t(as.matrix(assay( + reference_data[, which(reference_data[[ref_cell_type_col]] == cell_type)], + assay_name))), center = TRUE, scale = FALSE) + assay_svd <- svd(assay_mat) + cumulative_variance <- cumsum(assay_svd$d^2) / sum(assay_svd$d^2) + n_components <- min(which(cumulative_variance >= cumulative_variance_threshold)) + projections <- assay_mat %*% + assay_svd$v[, seq_len(n_components)] + clusters <- bluster::clusterRows( + projections, + BLUSPARAM = bluster::TwoStepParam( + second = bluster::NNGraphParam( + k = n_neighbor))) + cluster_means <- do.call( + rbind, lapply(unique(clusters), + function(cl) colMeans(assay_mat[clusters == cl, ]))) + rownames(cluster_means) <- rep(cell_type, nrow(cluster_means)) + cond_means <- rbind(cond_means, cluster_means) + } + + } else { + + # Compute a single conditional mean per cell type + cond_means <- lapply(cell_types, function(x) + apply(as.matrix(assay( + reference_data[, which(reference_data[[ref_cell_type_col]] == x)], + assay_name)), 1, mean)) + cond_means <- do.call(rbind, cond_means) + rownames(cond_means) <- cell_types + + } + + # Return conditional means + return(cond_means) +} + diff --git a/inst/extdata/compressed/VisualizationTools/calculateSIRSpace1.png b/inst/extdata/compressed/VisualizationTools/calculateSIRSpace1.png new file mode 100644 index 0000000000000000000000000000000000000000..904a7bd67c9979eb288b63708e09b19fe5a02608 GIT binary patch literal 21868 zcmaHSc{o+y+czOeB4fx@h9qN#OdT1MkS2sqi70Uzh{)M8WR3`xd1yo@Q;M|7JkN8S zWID#<7|%KO+56qTzu)h;1!ZE^Dv7*IM_w*S+q~eczwR^HyeDdqws#F)?wO zpE16`#KeYYV%n9)%7V5OnuGS}kL1lW7c7~W{AHP#f?=x9?AiHXTag^6kYIun!5 zD<&qPdubImdguei^XD#^AP6!=S|~0qR#H;3v9ZPn0mH+?Uhlg)7h7;(#Jeabc=F}) z7fYzVc-F&huj)&6^hIW83o~P;2sW`Z=sz6y&e%U>ViJ_u`C~S}AiaV%vU{6bny`~N zcm>#nS6E+sWn#MY+T8f`CBKoSDooshRf_DDJn4o-AnLKh{8xQDmI}*{%pfPrykLJ& zW3^$q%oydQbGi36-{PNkF3};sM!s_uYhJ4*sx|(H+KoBQwsk}bY9j|;gfRAfzIgfN znMfCbX6F|s{SPaZ9qJE1w)6M)_M7k1-ZTBnQ(s5UzxP5?Qf2HH#@Cx88r=bw(e?Eg zdkgp*|IFI9s&-5f-sNn)5sZ!=d}lAe-^o0SV%E80d)Qm%X@i7*EjF?%iB%-pi>h$I z6geg|_J2PC*7pH2hPJ#UNqYDMZb^l1`1+Dgz&XzXU@oQkt;#(o|_v^6AJq7bGL{IT-ae z?b|JcC22bgpkyZOHkO0WmH=k(qUtT=*dK!W_SU6!Le`b;x5J?EatLEv<4ucBktvik z8|rV@)<>7P!|R{ow}xsj_z$g zfLnrO9$eh^rR-r*@@UieOHfB@W1y~ZKg2=~&WKR5W(+I=k`FGz%hk_@Ky#fEbbeYK zm9I>>DRnP$xdCfqiHSxjbGOd{u%Buy#QmxX6fACIPp9U4@WiymWjJ?RLeGs!XLrS7 z{LxL!R^U(<+}ZIM__9*P1Yzfi*eU;83=@XCJ_p*!9ZSvsU;6y(ZB;ujhv$|vW`vC# zLN9H!)XouIa^+mn-&~Iz?l;L%(V3Pd#Y{koCAT2A$O6&3I~GoACSLC4UHL+OL#ldL z|MwR<^gbY|_J8~s0r1t&UEQ{ZRjk&@-GxU12W9`=Ks{@Mf^S3f`YgFj&U?M;yz*@7 z`SMPcvhPO$&(qyw?~FVDkl(*Dh%?;f6mB@CHjZJ?79z%w6ttt!htV!ePo);(4ZfWB z*b*K9?7<e;726a6>b#+0ai)6{9an{g9ipkrD~h zdQa9j2+Zi}Znxa!e!Su}7{A9!?uLfoSet}J=x5Wsh8)rDe;u|QQTS}=XX!3G@Be=2 zDiRECc5|UeY%ts})LddVgWl3I39}wTk7nW`eu#kI~*;^yBS^hH< zTrD)JVg~dKrPHv8X3%^gRKuLRQi6Ymj`TL%M8r3Vd!utS^pb3K7L3n_v z7VD}46*4=`;`rc_ee4N~5V!w!=O+$tgF%LySO$0J_v|#f~)`SW(Llu>ECdsay z5`>5U5fl1uIpDzh?*T+8TK7zhwRT3tMy?ng*jdECH`FK$VUA31JB$0=j9izJ;^(d6 z1>WwhTNtoMSk{!U1tI0TyA>$q%m|Ksa**A+RGVLie(^(!B!^mq<0JTJ1NP&B@+;C_ zsz)RxMQcwCIMymrPaQqncL&iStC?tQ5*|371eRW(v*G%`mFIt|(vc8QWKj9|OXO)I z(K*NVFQJDXmL7|ZgLd{yX9jCsM|emZ)kFxd!83cCRbffZ!!n8eLa z4w}@~eNU(}bp8&uLnarS;;^$6-=p5K&Sw|-!zKr7p;i_6 zLbuEuHazNDBq&^tHP1R%Fd1@hH-jen&)^MkAsea$xuffEk~>$*1rNfYX<(?+=|+k&i5*_EG5A7l-xe|h&B7hr*|f~`SH zcHu$g#-I8Ztc;U4-J&tu4Cd&zRL@+bJQ9*=K=SuL3L+PRQP)g>EO-$P3TFjMFa?qsvGI?&&;R14vi{tY%QzSWfg0PnfXe4l-U*vbF!0=d3 z-WqdeY=@PRg$6*-YC~E2L&&Q=38$WWyKx-_rTlFs{BFP_m@qHmijYqbs5C;TZq*Ib z63eJ3Y*`!K-Tq6Ai<=-^WN>m1tz0)22WsZ6)vI&a^;D~SpZmBsE#lqJc*i{l-MDOg zZ@_!siK{~<{%gJx3GLUPcik#sg$)mq?jhbAVH;$X%ZTB4+gk>u*k2Mv(&)J?omrjc;AN6kJDZDNHC+4qkB9= zh=UQ`Y~iCaL-<8yZy+AOU-F1_6Yc+gqvanMx5RZ zs+c9Bimx_(o9e*HOW%s>62kR0V7$}8DFyDs?;5C``T3%@Ga9Yg(U9T+M1ftXhUgHt zGT0EZP(LG8)17pJTX@DYV7$`4H-mV7yj0H14r$Z68N93G({@z}=vm68{}?rYguG{zsssH$ZwheE1ep z;q{DheE|FNv<6&p!R6bY?HWC>f7go*qXt+2)r(l9TriCa2;+qjBjKrG@cC;XpJf&% zX0L7-Y5A;Zo!FHx#fa8Sqq4sxNdvd39sTxDUNT+lBIzB~12;ar)bhi{+&ckOf3iAx z9Dl)LyZka_IsUoLBH37uW>GsBO~Pu^fawP7VhGfO&$c6@71RyRHpF5>`Mlk4*pmXW z#eo;r$gk+BmQo)9sYkPC>GW*Vzk`nh4Yyj=>H)7K#iQ$&noyH?D1A!VN#A%wxzf@)Y~Ls9w)srMAi>3m+_s{{ z^=7qPj4rf#wO%R?T=-_G{(e>kvV8guirtnN^7KFEa%OeVH!+~|%*vO2*2eGutZAKe zR{gI&5cmgJA$j`K%mycrY0BWWQ$I2ytZFup&qz4nQ%y>8Y@hAK`(E}D$}A@o7IR~qOc58$W-K0k~Vnte;$|UNR zRe1H=gwn^nmbLwO_#y^lW35aH58?4YMN8Jh(GehZQs0E@NC%O?f-f6_EzJW>YRZ90 zD7|-7YZJG(gU0Exxxaky0xUbgUGyJOve6c(E69*hhqFk=iV!qc*-sh?pCE=LEgpv? zK6x2q;J-CtWhY0M1b54yRk|CPRdYP!wx1Z$;-XOIijX14-Sd*obq!wG&+0h;;C^JE zfb>5++qkQ)JAuOmr;c_14OBH_LB_DzVqcS-Otcm3g8C}|*cs7z-;^GZ+l;|k^V!M~G2EkqMRLwDJi=bjYYqgR z;ojn1$Hh^E68Xl?B|!DEuhxvMNH@JExm8yM#?K5G{jk{CL;6{g$P3NRU6W*ir7Qh# ztTARd`YcVtFUM@M+L(LI#H~E_{{CzQd*?CM$-xs=$<5d=0{%qn4B`{@S6K(V?7mey z5}ouSW8dt%C(HZ;%3WvzuO)$F8KAo?;_rg_tv5!D%dnsr#1N=XYl;S<@;VmJHP&Z~)9w z71wxx)V(A(WLbL+_ktuB@Pd*#9cQ&q7%u7FyBK}~!QqNc^PY}yuS0w$S~(H`dm1(T z%z9!~(~1`0?c9_GW|8|i0QSxk^K6@Q0d9(qrHXj7C4lZ)`*~^}sWf7Z&OTgceTy1! zDC+X~_(OO(L#pQB{lPT4M7%~lnMl`-!T9xaMB`+myo`}=Z;T+pFOq1}cF;v087uy2x{0tI1hn+-B*cf0)?`hPjLmf@qSAaq`pn=}O zx`?Uuz-idDeC}g>%@Ek~;-FA1@lN~%ko6d=8XmTxtoTy*E!)HC-w_jp2i_qd?0n!) z=A3tt-JVrn;W&F;>kt=Di+?=~1?wl}Qt-z`SF`~O*G)C(^7Gbd@;ZT=s!&U0RikbH z(NLty>vRx}=PJ&@aM=WI7r`*RY?tC{ris;2t?#x+@eDYKa_%l8{uS(*j(n_O|?`HteAbLPX2xH7y6^4K2E~_G&0ZyD2_$0(jyGNL@y*)j>r}z!H0K z{erH6{iCC&R+q3E@e@ca0gV3uL}yoSk|j^Q!?PS}HZ-HAmcfTWQoLEn`$eycKD>q79#N55lsNW(dA3Qw>ri0LtN@ZIM|_#F6Wl zjROa1!kzS=n2nCq!XcX~fn(2~uWUtyWePJCp7BDAgQggEhVCoL=*xS*Xmf1S`I195 zINr;`nty>TLd!h%LVz9?`6)40=b;UU^L3+G#-S8{Z3qck3P19c>PX zQ#Yy}^EtyE@f$HKUuwQ;ndz3#W*_i;D|O7fS;~J`{*AqezR9hJxL|*(9yKUEfFXa~BGpPrF_i*{oyaVsVn|Twb9uqEk>2 zwDBD;TiWp3WqtAKMiBB)9QwjuR_>Pf?ttc}fliuM4;S14N3)2(>;Rqi0XGj^*qVF= zi-oy|hChBzvlw7*`(hJa9t!FTb=yXafX5q&?we<5(LS1MnK$Z3Gbg#9BKzN7UrBAC zT2##Bh(eZi*0S^Lb=)fSt&{xuL!33U<(IYx@umhB?;;m2>*OBW_y7zuF#|X~YODXQ zw{7d>7ho1}pS=z1t$!1G>-4aj1KvM$Z0{;2pk=xrr%K;qvO(-aLy;Q?cG2_-EHbN0 z)1``KAeoutw_WLjVIg0dJuboxi$v`rfFk#lL<*t|O8xVd;_l;@fF=@=#dw1p&#<6M z`17SwMXHgm@bHtxZ#%s^2NJ&W!dLH}aXZiTY{l1m1rx-{qW)zvTq|RZf0#SZLrIO= zx(2GK+kh$;*PS;wQA1L=mf3$EJl_RuQclmM_JJMyVi^+E*np>*KW;_bK7&TyfG=zj9fi%$F92V1j6E@$tmkvZnibq9!o`vdp zAl*9i$36=B4OxACqNjHQsjo1C1;liD|04RBGrjyeDe&mwiYWTn>TFyu#N8vA4xYWS zh@Ezq#WQC5a)ry1uqo)>asy@;X{?0JHI>9pGQ&y{xW@NN8G?<7=Fb@|DHI6>`duluf67zRC zavLu6Tu?r)3f}^Olp5??l}&=-rD9-|y_RUPhB*kg)Os{mKI13+Iaan#z5?9P zCdUSdO&Kvd3~QG}ylLm$MPzKL!l%CH6+WeIq(Qp%hGc8g!dS`z2i$Y`x!9DlGCq9{ z`v%ZY1Rdq>7$T>_(4I8TSO{$Iz`ECB8{0846yGzQH4}ulSCF2`vR=?gh{NDKa;uRb zN44Pl=o=5}NiSl>fs&n_IO>M%4>nP8VufGcU5`U7zl69Uw~l0^DC0ES*;U&VI&-iJ zJ>{=6RhPlQLU1SicWY(a4e41)sUe|nw-b--1CAVIxOG-n;pts@E50%0 z_G_okOO*^@Lus=SREy44=tIK~^uDs<$c*a22>tyzC8`#_j%y|mjXi#8j!06_KV_}eGpPr(eX8TCBlt^4F$46H%7{RXqW_CJE@HX!TqQ?x#b2MJdgp)+vbW5R5A z>XB)u@Y15=F6yC$krgc!S=vS|%eHH`#u+{E+JT)QukD;pY3&{n{ z6{EdM3Uct*a8vr_ddi*BeLq!`;K8_GWh zzYJ>6p>craBRN2(U<^)}hN0a@nMZJfJMtuL+@!ch|^RAUcwP z^vA0Oj@9igw>AlZ>zJSt=iU8WnMJ(gk6*xHpMMAN6#*u*%I~KZI<*_S3_} zI1MsD6+kb&f(2RH&h}@zTvUYW+)<`{IUXADz7MfzVOeAtedp{+3cwwP%6>r)Uq)hQ zX%cT7FqiC+5g2E7{Ku2-BoeVLF)13`8A<$p>O22nNc6F1r_zXz?bl0DspDs;HMgXnZ~r#n8(RYRG0K5Kv8Qn#S!NXkEgZAOjZ|++(x2Uj{9? z#BTpzxevxul6gK9?BE1;aTgk}(27k6_JQu0B7zfyHG^a2$FC772pJ<<#qvkl-_C!^ zq?T?g4@gA#UN6eKwksFVO+%it+iwysz!$#6jaQS1oHK;ujDRP{?Ude3;#?UE2Vu?7 z1lr&z8JUiR|496@5saRPiSv+W8ImGZ|FKu4cwXWmtE+U@O{H&e@B+eRo~i~nWN`y* zYXq!6Bj|3y?UrkqcTw+oDFVZ1?1I@tNeSmHkq|MOdE$2;mQ6p;WID69VRdH`V6$Nd zl&5Kh@{5o;kWumeBtVP{OB-xGEbaapPWMU~ZGXf$UK)$KtDcjDrR~Iw%%2{YR~+9b zegR+~5JeEX}S`IHH*)k3V38X%{E5Ln@_ zI{1X_Y+Cf~8}Kb%8$O6DO1fxbjwqZSu{I|gL+76+d_ZjP!0DV6M6(W?0|X?NIwhgm zB*sF-9MS%;DOu9MlD?I4juD{v6uhv2ZMuf3>*4_hzi7%cMih&I%&6neM~~=)6#=bv z))&I(f_1Fxf3i$vYmHiOzJDYWrUySTuPy`lr1#{WYCt3f(Ac%F^F*ZF(GDWAJTZYq zv>Emm1E3{kkVKUJR3Pz*&G#qf;l-^>>6l~^vh)oWV`rWYONP?LKp&(Z=`nzPv~r!q z0gGX7p*%oqRI%{lh!`Av`M`z~VD??|;|P0qvHuH`7S-1!zWsTC=c<0m^`9yu1B#3Q z%M2IO6U8+7YT{9`&oIQu_VTU1Z_id12`m|qIu*39U^4e`)11?9M6@JSa2T>2 z+75=X+B>l{%j1D>&|A`AuNaES*U79Vy?Zx9k8JUIhp~9AP`Bo1aGQ+GHtgLa(Bs1q ztr_CEeQ=$S1CX5zI>W?6u$2~j0UOn$$^>ZsLk>DY$axM@eaVEbU)Z93i!K&xh2)|l zgRAoGig98=nnzv)=u1{L>^?n^a&P9kEW^#7lDN5^`TTlD<)dNHS73>FsGB=iOla(} z*omtyD_Rfhdb=j=Q(p(b#fu(hU=kY!v_-Ut-+Nce4G+LXWb)!%5Gj?CdS?9>zVrxr~Q;Gh~UP zrYWH{q33c?tc3C#fc+{_(iMkuwa(cQ(aAu*!}PuBD_REiREwqa%AU5hP}Ly@=g{?x zC-@hKVr(Pu&54E-n_2;nFFCunO3l3Y6iIYcWNM^2)jY6XS6-CaG=iw=jBOaUo^ECZ#?vws?fwLprn|YGJ!}AxRkZ$mHPVdWyYLJ`Www7qZ5_bFb zKxsAPDm>J@`2l#HHAP#BJOyrO%lBPSdE(vc+vEUH22GJ`l0f&ao=;nP1s8MoEq;kz zEhIbk?EM|Ke2=NSj%|b7wAGu!?Ydj4!v6iST)?R3-5KOobVoz5xatDKr8oN{t_Ueysyg2P8hxfs*&v_V21Qk@z7crv*QqM9Mx0CD4EqG~FzVV|O+X*%@M* zu*fzh2(;P#XsGzf0vz=<_)jEl(is{%oj$VOvP+{-0(HBNBNftcas;;XHPE_GNVx^e z{~Je2vS=V&go~*g=ip8G@@KC@&KJWbxfETGZ>eL)JZeKSPxN8 zi@=K*{W;gnIXgN`S$yyqF8NT|2FajirU5&H=;L4R?@|Rg3W#a%xUI{~gNj};z6kb! zIgn4cI5zOKMOEXl%|k8`DwRHN64-IxalZiLyG=k>Q4!N)cBU20Iqn>6k%$)|ek8kb&_!So@)`MflDd zMjCP}UL@)}aoz!&1gc6lmE*v>M@KraqmAXt@&Jbk8*Jq~7Bx?NL(LCCvR@+QN1(Cx z>s&8VzvjO?%&Rentxg!>M%~GrM5wS6aR>$z?aI|Q$+DsIlhY=Hmyl~}fCuV&e>;`Q zO|U^Wu@cF=##(e+30N==XnGF3yDpt|jjNIem~?)*iitrN&Kafj3eaD3K~r&R)Kp>l zPOoTM=Sj*!40z6K7!HpOe`9c;p+G{aiC?BE;I^ zH5Hwa%LB+2dHEE8@hHZT!#imL3lV7Q(KZoU&?@0xuoZ*ZAMJm;iN3Ee3KlQ1i;rC2 zxk%sNgIKuPXsN~_i<~b*@84s{<>xjw{*3?FP(kg1qlHXK0jYNZJ@bj*1X-0km7zp#+5UwYPRbp`Ky zRbC#>uxuDC9VfBm2dnb-pm?a=NOXxPo`#>|g^TBzZLOs%lHP60Kw{tBef8yiI?p+$ zoUy|~uhWy_WDY-lg?eJ)j-^eA=mxQ7&R1#zcS&DTvf$m^00-bJ<^T7xTnl*riEU>B zt#u7u>gtwbOam=H1dgRH#7jS@oGo#Dz<+jDd5?uWkUP(8H@3{v+jaohUi$>GEsHQx zd7Nw+l@>>O4pd-9u~G z$7$ZF*c?V-ZY0N~5NjhUjP zWPgOf;iKlk_n|M#3&fYy%a0*O{r5PX4V*2iHkqpTAh$$gLpEVaKhx9F zO~rxX@Q`-wlk_-eS)gQte4O$_K;pw$DDY2Aq5sRYP}id;06p!60V+f>pZ_~_Vf4b= z#;PaXsQY^PNfU3b-7na3buii|-vv|3Efd=`rLE(N$fhi#!JaGIT9gSieK`~9R3}zc zMPzJ_g0v;;rcQ~Bx^;pTR5o+Fux!<~m0EHM&WKC*MorQz?h0h@4BOp5)v;ks4#`4%}szCbcsbnnB5#ITn5}tkse3ZC-7c+=bZ_%TSOnr zKu28huVcl*?)uWf*8clHBUV-_0#ZKf8!cknI<%idX-* zY_D|_$>6sCG52fI_np1REc02a#fuE4J!L`Q+ePBks=`?fRh@C{_A8M$k;>}ii{W4MARf1Wt}sKlu!)5jTB%m(TQB%;@RH?1ipDI(P3S68p#mYwsK)FH)u zU~;EosDK(yWTX<|P^Q1W?xJW&PaB3CBfx)h4))?oyA$FA4MJ4Y3`O)t?}y!aj)z~GLA*-Kimf_6JtJfmY4^@4evvYd!tjn(P>F=JA7Drth{41* zXd&>dQA?AI*eNA|FAcMPU;s=sviH39jQ35v*nU#s`SVA!ITSA>&G?x_{T-ZZ_LKfH z`T{A4fyK^tf!dS>3bh1BS;)w>gHKR4*6is*_$?~)0gQAp3KT{IRTB+{nKS# z+}K&4_$cn31iGr>*dsC6twZ&z*2$CPNVqKXc+UBMP*Rm#sexZ(%vV#6%M(oJobF;= zxZG`C&6W>rOBB=v_o(2<9`a8S+~M$bdYD1ui|~(+`9>0rP>*?Y#(j2EPd%s5#_-X} z>c`I-zg@sN?L`BJuZxU?V05a$$>0?Frw_2}hg`ftJ&n^nFk}~VRg!+B+LuJdiFL~e zk^;YsiyPz=Ud~6+I(}Av>h03Y8sfRqHy}kTuL8Do`h3I;QB!Xv&TSIKSHgF)3VZbY zAGTpLZ(R@Xxp?sxG2MdGrbe-Ti9LF1ub`v+1w1fkc_L4sL}ytmdAL6|D)tE_6tp^mqu4Rr^w#ydp#%`hIRu zhaDMi|C9x{n>0(_d$tXYS@FTx8Z>;|lkEzDl$D?mCOsI>Tj|9(tLXqdtkVVxb>?w< z8W(KKMGj2w|a$x?T@goKXG=3V50vm(tZ7c_;p)pI=4d0$KxF zJ~gg10K3lqp$rP}R^=ggjeE>+f57R#@-@V@+N1A^whkSm=AVT3simx24_I#Iik_fvIlBAMx29Qc?)f%ye;uc| z5my$0m7orN&4}eSdLUYiiGlTX3%kLkzvOnlk7wiUPJ^J^>Go%!+k4OVuxB!8`5AB^ zsAyP?hfo%bf*M2m)q+&69ZlW+v_kMlkX_#gW)m(RI}a1C+7A<@0Y7S9H&li&HBQz} zQmWi~7@=VI^i>qdG&CPuAWC(E5}t$LzSsiHwQqU9Y@NAQhfKBw&b3=+hl=H8?~da< zSrHL3!`I9C*WH}O@2nJM@5w8X+mx=ArRYZ(%(%?hlKgP=ZLYDWVE81VnNk`+jywb` zDMR(auvMUPg`*`}gO$DYk&?I;lXy8A%g%_ttW3jcS>(xU>qZTO17=Xg>~+Z1;$i0r`Wp!BB#fXD4Hn}T*_a>&OJDPCaC7ahWiliRv1x~O|vAxq! z=V1i+8~Swz$TrF?;GT6)-+2SNhap?j$9~;QIj6mp0mkNi%%*WZYgi>e!ORlS+A-^t zroCZS`uy{*k@7_oM^ryQskon8r9&bnP!@8hRgd0JpCZ$4gsl^d!mYnDy3WBTQb=*V zTaqcDPmoz~?l%taVc1fy550u@q@uI0ZI>S?(9hZ8s{qWG(1gn^1~&(6l@cCCY5htA zQtL>SA>=W0>GNZu?$*ogD_e`DmUCDxh&wimn^M@JGZ9@C7?>3;yASBJ%0G`NrL=N= zaDcFtlM*8Jzp&XE80*s5l=Xps9zTMYLeNy>3ZfaxL$rN%ZbuzZJD@-iINKuBy=Uq-|(EOEkEDaL34;@c+du+enZnQyOHpmJY+64ZV? zFiQooGDJE;7m3c0YyQ$?7Y;rw)PO&!>-lj5*T0Utv8}}yce+N%)ZlX7a5|~4`(vOWE-%{b$~WJa;bnVeFbk;_;i9id|RLq zxxn)G=9V-+a$cBTPni9+BaTmKE^AY{T~Rg}%+qU`AxeEL!f#(JR)s7LSl|c7DXM^; zl;vw$F6$Cjsby8`N+${N!~mSr*xA$S_#T|&K|zc zyavQSLNX@I7fx>#d@L;L`~+VWQ(Gmhb)7|UMKsRmBpDhf9xY<*9_7*Fxg;qLyhc*O z?;?TCkjxw;#)z0_Pm5K1$!XxB=DmN4`7;L!i8Ka4_4S|jjOf3UhB#~j+%C@g8Vbw1 zwfS)5vU|meKg#dCfaj)$7gB8r%M53&d0 zjC4&`F>B6s!k?1eLTU^!?`XA>G!*-Zoi((-)RAd|yg#*O?)!V!EG}pXtizTifQd1Q z@_k@}WC&=}DfU}_8pSP4`J-Yw6XZX`SNnW)J^$c%3#V|&q`eyi={^^XaX-sF(c_Ee zL)=>83V;pEj`oS>H?O63gZyc;T5Z@X{|$3I@~U&AxrZPrr;Cu3M_%f13Wu?EhCOgX zJK1v;(Xl4Y94OemjiwudPNYg+^?GOkxqT~qe+|)b9-N%Q##y1dqB(OJ#niA5zl4w~ zeAV->vdCjiyt+eWB{_Vaxe42?-$>)Ui+m2tQ+j%Jp0~@M&uakcf8HI*@bkr~7R25;hnG3T2v? zQ}@;qOWsKX?~8sJZmigHGeiSEw8v(aX1q9yN?AjOwypf!e-`8P1r~UYB`EhD7eG&( zr3tn5-8CHo;CVIKt4P}229K7dm%KX+xsrasu3fhnvaJ%WA!!iU-AvM0c%FGKLMalI zLGhppqcuA7T#O9&Fq!DOVr#{)+g%1@Pj{mkZx1Q>(KynY#Ra_mEb(CE1Hd-O9Z4y_ z0O|CC>XSdD8L9Uk$f5ZYQSP#R4Y0n4+9obKqtg#Piyjx+v0@N+Vb}!OKCv5o2HsYI z-l9YnlzY(m22I;RmLU_$Vcj6Fjao$V#KjM`m4AUF-Oa8sNI4G>cUfjMSg;*?^eXUz z$0-{b9Eb7%5}NTiol4_89TE1u|4Oi)<*VX?`cESw@ZKYF;3SX5k0*XfvK5B#-vL$m(Bl2}2Zckc zxr|&+*QK#3SL?IAxY==oU@WWala@_a(Qa3eDR#b!T`qylP-iT=gg>XAY z{o$6rcA{tHzKKW^B{Om|cQ*;<;x?Q|OCHbwv?_ZunRqwBzP)KfUkN}yGoZpF$bN1| zgitmd8p;}I9ZZ^1o1amB1^Xd5LB=Ayf!9|JSe`^RkLSDNa&A!1oXTPkn>TJW z8=Jd!ZXwRNHrvb)QJSl`Q_csm@WVCd}YU>L;l<6ZEmY;YqmWXRY1BNkH%nU zpNfK~{3`3673WHRcu^LHH*wPO(X5hMct-5=g*nMRYrerbbY{))kU&2Np~Ra;i*KGk zfi)xU!Xuye0MTQ3&J1GT;To(?unBG3R~7C#g!k){rKO;8sp4BjRzgkKyT;F4?k#Cw zvBiQ0-qDN%j`ww(*}DV$ziOp`ALWdf<$>YxZzRcw6^N~(^orJJU;@oj#FGhjY|O<= z`c&o8WuDvhP++R$_6;>d_t~E3@%N+dLW3563H8*)4DNlAA7P-BuJ!0**xi4gY zr@^K=UT+RBHc7^0tZvo@5Ti zcrKstGYZR^$a8<+L{V;l1&}PkXgi;{`Y83t6*P>@35h+f5~$!< zYR(R3x*PDYz%4*rr2EPrvOzu@)P8ax3HIA1F17$+4nc{5p))th*XwL6@lDV zl8!N{p*uCLrZnMUjLX4?ajTInd8--BrA;urj!6^c!r{;b~-{U}BV-5GKfk+rs({ z=_52I{$(+Xp)iarv~>+3%+SBaPU+*4EUd;(@dZrkM(K>e^)%pZbVng_i)vBX6h&li zX$_qq%crsLLei^OQXWMLgBKI;d|&W@wEOyg-mYVKr=KJ}_^LBh_!7{U;^Ds1FQq&STt{vEeKDGa zA6jALBYR?+P$Kx3po0*8L%L{`bhjrpRUb$U`c3umZ-2K&C@C?9+6hQQmKWUqNe{Vl z44-~l6g`m1!)W@Jhb8r}=Yf4d?+t1``~UCf>eNn+X~80*cv#05&IzC_T!!41jG||h z`Dd|mDHGd(t`K<56)l!FQbb$h?8}qeC0bj6N5*ATpA$eD>W=iHrJ(8u;nf>d2L+Tl z^=txHh}!4W&+o4uZ=P|>^^NDh+h5B19;HzUv^NkeP{M-FK7jAwIb=N|ivcCg13XS^ zmRdS=373LQp;nWwAv$e{3&t+Z;{V5MIQ7%qK&RwQ z-stcv9iWZ)nZLN%S)v+MmwsRxi>a}%@M*w4;M_|cl|LLfd)HHzjedSRj)JU+lf0`IwD1tO_*xUDO!sKsjKj>{1PMfPj1NkEE=m zxMD&7!^58F$2qUPfEKHg7|nBp=)K7$03s$gK%Y_Q;+a zm=I|06z>I(Q5`<)QQ$62163#s{j4Af){xbK~NvCfhw75YqhfU!vZ zpsOv8SXX8q$|RlI7w$DriESV_a7*79^EcJR6k7>~jb}Pm=u#FGL;ktp?s5L%WQ)sY z#A+YJVDx!^S?QJknPOv?JGYN#i*p$*c>HO%CV!u4Rr^i zSIf_n^al|Or@_j{aC0AuAJaj1u7dyu13a~QI&8{I*$47S22OvQ3KP&s+sIP>skWIQ zxnRUHARI0}CwGL)C`J8hq(Yw{cq(u9&m8vC@}op3QG5MO5>=N!!Qnij^IFbgB}LY- zKORkfpn;8HkYA6{nia+njb3*wIFmSmOLu1jMJSp8rhs8(LtRDf{#UloWAg#rLpx{5wcL1f@Hd27?c? zL;hhP|7*C`4dj8#&nznO1gR6-$fr{gOB8uKvh z`A^v(#J`BT5gAXapaHsi zz=6GDV>GIdnG?J*^LMZ_3xl$9Zz2OhpGQGOLj|g?kWSNzt@e)D|9>Rg1?5tMfs`)_-lcsiTQQ0C=3U%Hb*Fj!%A&=Gxz}CLV_@~} zFy3i#`5ff$!6R35(%AjWZ?CEm6Rk~xcvlqoh{sv%kU2`iPoda0AYE5Tby$xQ4pl&Z zBd4_soRrM^u|rW#r|SAuavcc9_V2W&E+x<0R{N*fqq2Hzr!=6~5)t3Wi?S@5u7)n} z6)V)ii&1s)fqU+ChG7DY`W9F_4&U=WKqSeX;`-JVk;?%Y{TU73t#Bu6j9|A#rEC>U z_~{9XJ!fZF#9=BA-TvbE*PR76O$nJqQ-bk^CXf!UsZtyqzkn!D5TZBpS%6NCjk)e; z7h#?I@W^;ofaw|j1FJ(~T4s<8S7tZG9L7JJqCbYp#&f_J#P;xd6qb0gX!pYcozI{9 zBRkqr;?HY0%t^vcQ%*im=|II0XvNdU-RZR%?cB{a*oTtT_AhJM-JON;{AwqHpo+y- zI7FlRZPPgS*R3udGT!~(w?>ECjRYW1yi>^y<(_!{AxJ);Z5uvmVIgi+J-PWt_=NUX z>$6)oLmKtP$xcX{j(7#iCa?R*sLb0m5;d7^$KL(YHM*qZa*OFnK&_DeL3zO+Onkp0 zwf5Q976tUJ)a!~*i(6AU$HFPg-4vAifBz2RWmvj3BJVE<%%!kDUN>$4WCfOv)W}1r4w`~;hXrO6u+gt)eZP+u$sMmF(0!>zrr;)xxWh_^2 z74|?y(%&mj0lyEcm*7kI&q0I-rE#q9yCT5%GcEzNRNUPi*1RbsumgVE!_h;FQjY)g z14C%8nQn{=L8SgqALkv^^wRZl=~YUUB1jXEPCyW(gCrm#s9c&911Mf85`iEflF$SK zh=3a4(j@d=1qqTulO`CEqI4b3c4zj?Idi_}Lw`u! zcCZHV#oAL*MHn6xndhR%1n!(&`qmO$_~$kHWf8)Jahn1<`lv_q`%B@fgz*Hq4iWZ~Mu zS2V#yQE=Fp8Gi>%@=HwGt42@P$TU0Fct>x1? zb*w0D;)-TI$}`2d38)$-G|knpLso=vm)cfun#4n;t|Gq)ayh(}!nqljxE&YHO43Xn%}8U{pHtRl@oVL)6a#NAyn->raJ9^@M6k4{p32Q zn_7);8LZiP*S)eO!|l+1iRykAkNHGjzPr>ifqvC9@rY&LL*Jo1f(43#X9Cu>AixtW z!5Hc1YER@fV$ppCoz^-!J5?Q1tzf}mDbQx*Ad71dNa^8(2*-U6JD1=xynLFW;DSug z&>UxYo?KJbVT!gyyIWVVek4$uuN-J zlBlG#nW6i*OWA2M&*@ni!1X-V3)qQ%C4O%|D1PhI2N#V*(zHvmIEC@sEYr;-Z+ai? zJ@5UhFM30h-01oS8F{xOyE{A1c;LXxYgI{X zYMEuqe?)}$m>Cyv+gtN}%(k|VAQWaFQ#k&%m^QvJw8>;breNL!g0LMtGx>A=w&Hv4 z{ke!g=SVOSup@k7KePj8S6ZMNf6|BiaOgCUE6q?o&Kr2Vw~*&+X(9R43J)Zw0ALR+ zp%fH;anGUIGCeFzoO5aN4Xyo(^Wi@Nyk|mtTGz0Q;xWxiO7y^!4e?8qFLN1jdwV&h zkjA04SUpgBY6Y0X_G3<55J9n9L-7}{`)n-$#n;Z>J4|S*Qyl`+qizK$|cFc?f56yFNplv z13QDCY%)hA-nz~!zC9Sp((bw>3f@(2Fi`toyp59HJ!1SHrIy;_-uTh z$Z;q9kJi(iPno$Bi({sI|z8jnM+m!Mwx&;9+0?&n6`~zCa=xP z>1$97?@X%@kSt1JXr0x?SHzXEkl)jIekN_9d2OWwuP2k70nMOj*%QH;7yVa{_TU?D z`TXt&3_T=iYa-#Bruy2l-K@DRhkJRGEz+1y^VNg$SA`Z zS#Wn;E`9!j=UtU!h1@oCyph8vI3 z^CnLyt+C@k^o3$*ULdu1F*U-w?pQsTZ`1ZCgGhW=WgowoOI{izm@$~RXvtH#$fw$r zN7=`6W>C1L9$opo@0NP+8DhULoN2e?spZx4lx{#SJ?!XVf)WhPkoukp>GG@sy?#R= zN8ZiyB>)TDYV6TyEMr{i-UYE_(?frLEQ~%Ezcet)K6|v72~|TusSB)`?Z*OSk3&xW z^Q@mVqSV{pLHC&|?|0MUsI zRpZxaC3c?Q=?(XdJXCu9SpAqUkF8RTXm>A%(98m9?t)NWlxUETMj&i#8gf)NuEtCI z=_{yOi|*C4+Q~!R-QdINQkiFZCbYel;f;PinICS9dHWNH2acthoL3dvJ_jI+Qf5abZa!uo}A7Mr@BMyh1s?hLq|}OB1=&HE=>L z55lFr<(lhx8>&`kHPSOyxt<^Q9o=l6)fun27?%Jgyi=l{1PyJ2E%>H5IMUa*m^~y? zCOrLWEq1c&4oQJ=ZtntGxlutwVND!-@Tjf#PqxNBI{6lgonR#^)w&~Rqkec}%~K!F zPi9sF`Bt}f!mpB_#P`mHpBYCg8F{ECP(RLtIYkI}M9Urse=^)kzhqPKd+RJy8e=^)#BK)im^j6rq(KSky(lu3!J9b&iL?p}5zI@eV; z=`x`Fb*11c*JK%xe!Xp|n~cXJ2?UAlS>Nf*p?Z;SN!{BFe~~|{fW#2AuZ!|-eK0N* zX<=G?5l&}urNFjDQ0_?;*D9S8Wuu-7B>Ut!#p%-tR(32it}Wa$1D_&+L00D=MW}=& z-xN|x$7U;9IMcm95m9)&7;B|V!qD$j$zcW`wcPX1{`q}{xP1ZoV6X*iT4}_=>4*iC zXVaCaHc>!y^NL$+sMRO&S<9~F1 z(l;urDj+*}1exIf!Nc$(J8==mQaCgJuSjzA# zeb#_;8hZ|dTU5KovKRh66f_Eywwoj&4Z(l|SxF@L z#^uWM@bA`+St@wi8WUZtE;_G6eTJfQAKAKe?GD?iy#zajE@Lx*<9%ymMA&^-gdVY+ zh!o;`ETMmYd+NjJ((b$v#pi*zsF{3l##e9fMoOTgAK~6dOV7W2 zBP3>F4-*IdZ%owKZ`j=3p(Wg^;o_80b2b|jr~G<$mFTKy=G(}ldN@sz{>?6Hh)40V zJlC1D;5h!$f@&*o0>+DO%V4o6&V<37w`a&xTUh6KoDz(D2O??di&>yXYAY~!FPBJT za%@D$WzthMBkU@66R^Uc|y^lZ$Jr=()_QTlqiuS?a{=r}$Cc&^; zRTneXK@qO#m#ByGoZw5v|DL)YAsGE|zVEtj2eO1JYAZU><}R=@ZNIh>U`l>rp3KWExn$dS0A z8sFwko+*sC9?dD`W}2WorKJhw;w%;0LoT#`Q+t=Gnd1PGyxZa%Q;cmGpYyq$KN&5c zKd+nB((Y5pXw zQx(TVuK3$fL6I7}z*8ZqtAAvX(aO)J@PLwnWo%lOOvjmQkH#p1511OxplF=Xq{UpX z0!IV{La0)tUVZ8yi&}Gw=lBG*FWtbFKRqhHCN#SCO$7XB2{g!@nfBIRQVL5iSEfEH z8t1JItz$xx?i#qW7=^kY)q~?z>WY#qowD-sylEx@kmhK3c4KUsJ6p9FqJH!;BbVc5 z;eblh3urlw6OI3NcgJec$r1`ah|-zt9nUfKgMSceT+lZ!6%^~`cMV-e_Bb+zn(DTJg%F(y3VXRsCJN`Om*b-+nHzo+}iILulSCaaZb~3ZVPH1bP!Y)cGTDXH+eO zA1hK1W~wWQkgaPb(79OBv#~fzhZ7p3GXT~$b>wr$fp{_S7%a- zOxs?lFE{IEo7?nytYAsyUibRYc(T=w&+ASazQx;gU5~?u>38jo$(nj@m-^bs>qMD# z2PN{_8O%jgG)22=!B3X8hPQCi1i9fkZMvf|76;q`mlGQJ`eM2$^d6} zYS4OH@gl&u-W?Cq!sYyXNt&Kb?m}t9_?xTS6~$fha=3~jV7AN_zl^PrY^c%oRu6|t zq__x}G4wN_UB@)eVXCojR$re+j&1~L^I^@}z1GtW1myBT4*AsV&dTCLZpb}W#GSps zsj2zX8BJQLA=+giY|DYX$YjotF9+Msx;d7G5&JO6u&B4}KiwMv6r<#rP*za*5p4VA)9QK6~xCoW(y!E@p7bAVu(1tD)0 zkzHr82f_mr?~uzv!c zATag2h{!X((n7uHOk7^+d?x6a_!}2BF>1Th>bc+dNIR|$Z^^#@Ji*#V?81$x0?JR8 yh^jIm`1Ht5(_zN9pKFR)JDMWoL}2<)_@zGZs*5NcDBMsQ&{pPX9sx literal 0 HcmV?d00001 diff --git a/inst/extdata/compressed/VisualizationTools/calculateSIRSpace2.png b/inst/extdata/compressed/VisualizationTools/calculateSIRSpace2.png new file mode 100644 index 0000000000000000000000000000000000000000..c94156e196d4e5637bf7ff548eb760d4ed9fe2f7 GIT binary patch literal 27308 zcmZU)2UJr{*fmO#-V|O0q>F$G1W;+xq<0XcLqI@5iUA>^iYUETC3F;mMoK^sfkP2N ziqa$mLJNY_B=joPdwk#T|NeW|U0JLwCNpzp&dhnv-uu~;aMx6ih4B(21qB6*fxfl{ z1qBV9f`Y2-95q;yD-f~=ehInhTbNK#go;s6JddWJI0B2FuTW4RWGE;|_bDiD=TJ~^ z`4s@>D&P+Sca1G|{`~o~zP{ex-fn1Un3Iz;F)&iD%HRn>aU%a~cID=v<)5=+lBcYs){{4gH+qw4!^q>UQYEl_ul|%X= zXBc9Z^dJcF%A0>G%Hk86w6r0Sve)0a6^HWN#L9u6HOy!> zrs`F*E}#p%$VE?S{GVbT;w55HIot&Ex5!GC1EKurZyjAgZnZOibZ_MNvwHa1(O#d! z(HX{0yjhIU*u3$xG7r%pfz?rnXgk~7BjdX`5COxL`GoZKYO{@q@K99v;W_D|zzE;5 zZD0=)vM_-S?prL&Jgl%jGM|qzXqZ6f z`yJSa=={Qrbq+SI-Ory-un|Pv>2T5qcyTZ;=1=G~!a>d~ zt;-7W{ZFh()A4xaXZ{RGz|DBTvdfnjQOAJ|KiZNoS?P3a2*5s8#G3%k-hl`IR_R3) zyE|VZgz2G<;d@IqPrikp;#7(9P(GtZvO@C#MxTrDfRk2|$7SqC0YV!vGxNpBGPmlD z?nA*zso!#BT1|$3F)`+T*47;~@OHN_ecc)D5fRdtslHnajs=+V`>~b?&&|IIy_jB< zSLGa}2{fDeiGI)g%r^TFkDBtqOfZZaAsEs8%(2ikMtBZ0 zRasXEPG(2X+~v{9ppiK0cz139*--N1aQaytk+8*Ec3d z{Nw}8fc3efDqP)LU~~ieHeJTJwby@F`Vuh^pA1idn^^&3_`gNvISBy&EfFaAA*|eYvzEfM#solb3szjaxQIE-b3 zd(Aoxuy3-HBdwFm=hCg6fgQKzyWQo>!0~_k&q}rx)m#O<$$%}denokgVER$m)=mB4 z$v;SsKYkCU)T|xC=!u4=Yi56{6HhxORK%1U$P1ZA<_fJt@r#9<}owAxD+pV_4iN zsf_XM+>YC@#+}YY(7ucPDB>}@n0OAhuL;zRMm38y9}{u!ifahR5lF zDCOsaPBpBk&^`Xf6f^fTdjlX_MA=nOg861zAc4IsOnw<#|MaJuzyMPNzN~}dSy7c| zaq51wk%HgOEGKP6hdd2uEG5|o0$nwm6slu05+xcM7NF`7(0Cv!8mHO0d2CVt$K~|ZfZ1Jx#O)9d!p4GBq2)l3Q=GQ^$kwK@^3W3<4hTdOU3rh zSg8<&$NiVwYCa0h>kw&S`@dk5gQs0C*zo!Vh6_(!JfEGzwmuDw{*rldPFV0v*~W>& z*I#+n;miIan;sWhjJZx#A77$geAngqTb{5LNNVqpy5pB`mRai!t{Vlb^G2nZHOw8i zk7nVQ{h>=sBO}MHi5?;1j&B}%Ap)G9tac;p=&J? z-yGCi5$#;*_rG|5cpkxbaUsaIv6TbnE84)u8kBL{v3ozUOBRoKt+<3u(~6lB?{Jm< z{Y-p9SvmK$tkb(orJ*<}SpTMJlvvmwG81cnL4bgr%vSP?v{XwgVv zbl@5>%Q)APymZp+na#$p2#BcQ%ZQPB@>yzJJ@if)9Sv6cHmdknwsJ&k--M24&+4Vu z_j=jg>Y{DzVTxW&&3_etgrS_+YioUsRI zUXh;EO3ecUh`rVej*m|J5M^VfjyCS$y$GfO{f20Md%!(Ax<46i1P0csoD57NQ54F7 zc<7FD=rG1T?wU-;oB$W_OM&AX*v`O^8)0YVUf9#t&GZmGm=LMRD@NlD#sN`<8^v7K zc6E2I2*-4@^1l9Q+~=#Zq7QA@R5{3e#kiZF=)H&f`E`{7`#ASmR**nk$<+@%^W@d=Jve(ix*R4`;Eqyr|%>T8~Be> zS|?!*>cZl_yoL%1eOWlfj#(pzKU;E>h;;B7%K?LKf)zBM2Milh(ro%>1*&4rR+4s_i0?@^#mG@Eb|MoDtK%FJS>HQ_u~W12&oIr z<9WT6deN6*C?!;ckgFGAG|qJS?S_NTT1EIZM<-4XWp>0!qe=q80k}mvIvGtBo;vm8 zf3Ul}%*AMGxSE4~DDZ@M^g-zZFB|IEC4rMxe`Dfl=gv-Fh{T|P2q{sVrT@@7NDw$Y zZ3|u|WMg8m{8YBgHBZd*^ZVc3 zIB8dU{()@BU29(~8(JS4MdkPLNkn5gt-1Li(BL&Xm~M{SgwUN;_gr}0P~{yl&gR4V zDsnx+9yF?CkXf@oIC4qo(trLS>Dlqx40^yEq2`HBSpX>Lzt*RVzA2IFE6Tb&QijRD z28EA2iHz1*1|WfcSD{<-jfa*2=nyNSyo}LhF=!+8hur6I5l&jXoA=ytdFpv$I@}14 z7@2)bZP2h~$9{`cWIl8mT-CVgtu{;UJB2F_A2H2SJ81GPa3@?Pkfh9IN&Bdd=M6(@ zl5RmGD-POLFZQ&7^kr7Nlf9LQ0809=j$GFAl-E=n;bcAN-!)Z?b&(5y{_eztp0~4o z2H4Gnvj`451Ad58`fv}oBS^Q5@p*95&e_F@k=|BLVY@QiY}R+K-0)}q=_>5UTKvNq znEg(-D?t3>y1-poxJLJjp-zk>RyB{6k`FI)__mw5ek~0zB`-n{D^S(My~@jgy^G7U zxzS3iD&Ow%1A%UeBr4-tGWXtLBsru;5v>7>^rxr9 zXIZf2t81zc#j|B{z^E#QxB4JraPTnz%Kbrk(?&g^M=gu#50BH;41uI7W$bCuq0W8C z>VOxOlgO`(0YOpfF||+w1qsrP%5fio(ByH!0`K@uJ$R^dPx1F# zeJEwW4Fx7EtlI~KE{@fDd{LzP0>_8$w{nhaIPP-eE-uV1i z%Ga{5qr%f6+600gRgQ`nv|ZXqED_b?N=}CN2o6i%Xx8GJ?SA7tDLlC9GP3!Uzbu~z zVH0r2Vt_hLPl!MP{}$TBhK)fpakvtl;fzYW9y@C30`s&pg;n1wIE`^}$WEoc?2`+@ zinc3P3H|6Wx|=%1%N=evq9%S&!973hqQTZmUw(&A-yr@vnfRh-TOtHap;g^Bp85P% z@XLni)v>s?uaZG3R@tzz_u!qbdrJ6{$r$>vkSVQ7xfemRk{N|hq*Vo1hBr!}Ye;Oa zP?=f%?VV7cYRuGcsRUVvy{y)+_mvB@?o8e>a}x?R|7=VSh6dFXnEQS%W=hiwJ6}Zs z{pEu#n$EP1m%f~}KA~bujnYGI&2K9{LYm!hoLP)XR@xsrU$=;Iu{2BLV9;ALpg|&6 zf(q!250L4Qy`pZX+h?n982GKT#KeHYosE7>EEI}DExYYDgmP{9+f z4ADb2e%;Pevi=9br18vPX}k;>^)?H>la;L0-MAE>%@@i*yn&rZ{>7aEE6g}7nSS;rHBI~~!Q>nmP3VydkMsajVM~lYuxoV4 zoU0P*L(q6){2XHE&;4bwtL^1n(W1DqzhB5OJ(6;va(V!Z4y)hGh3n8Sq4Wmv^_Tz5 zr9i}wIt(zPbM^`FdsDymt#=o0q@{=4czTm4m2hx(%UBnBccVz927oM+m(~jHHpTn{ z>ruxhB(F+Y_A`%VptA{~To91v>2;?};to<$aN#Lb9KioBNxArriHY)z^%_gX>B}P(;*At z!O96nToKl~4dFtn-7Jat%`Dd*GVohse$w)D62qsS7VHycP7sH!-rdOQNwcgwpj^pR zXMD?r)lKf*>vxd6XK&Q> zktr^}7r%d|KKE;1iU&AQz?sWoOLY;1H^(iM6cs6M+!#OR&e#Z@x$K}c60qPyiRJf8 zdk?3hf*)SCcN=nEf#S<%VAcmXEw18-?HT9ba#YX0yTL7W+)0>VY?2EuC7CbE*~&k# z;ypZ|1d>vn&f4lBaF3XJZpFV#{XuMApI)?DYDWhTqxRqxE^izO|Cj z*vC)fsw@Tq9{xdA#n**%kFUH5H?h z-UlX|Ro%}5*0`HskN;xxPyPfd;@IPb+W$d!Rg6&}kH;hcCjTG!`9O^H>G1?MMmA01 z_w(i1!?CAcL%p5qsxL8!x=>z~CHwx#ay49wgE?2X$W6K^5BoaLSSERM>x5^eQE(U3 zAxxGqt%K_8CjP=vTGU}Ayd&q|yzxkq2x1W-hvP*$*Y+vh72nTqXFeRC#6E_%K#U;| zSZT33LPF_HuMy+%gh&AP0^V$nho0xv#U7gqL0PSNFOzGR?(UGVkGHbVQHaF`Jy@1~ zDPC+?p4vxFhu4)tf@IPm)KnJ4$mlQ;s8GC+zc;@=^2Mm+PY=2*L%e^yPt`|;XX4wB z9|jgRXzBoNo4q>3%VIPI(ouaz0xh555z;DG8Z5m%^q}kag`nrCG;uUny)}Ue6FVVj z+)H_o(SCWqr~`euD)SOTU6Hc>+S@=d09v3HI{!LCk>}EVT~WUfg9<6|EuX=UXo}JV z_~`Y^j$LxIVQ)G08%)`(KD3~jb_%{k!+C5Uuk2a_(=My`gHxc6R7m7U5kmzLk_FM9 zE(XqX{(Q-LvsD_(=yQ@wPu31dnq7UjCoq4;1YRV3G5rAm*v!Q zKFFq6pZ2EY_)hVnc%`4a62Y=6f2Rb)hJ0Z8x_Nmmu$i<I}Kjv)woSG#x)p9b#rC zefx?0HDPEnu*r5ZT_EdN^hfOFx_72tsE-X`3L%;s0D)HtoNGS@il`~gfn9(W#rtLP z;ru#dnM2;LZE3!9U&;FYL$!w3I#D!9J*t?g5_CwPPC z4)g5_%${^1`I(+PutfDmxJROyK3DT0gAR9tB__gtSsq4CVo#xyIknxrslc>GhK!7m zkSp192xD8l-Zgkb^_cCDVR#5bL?ITN;g7MQXY0{f|0SDOm!rB8A@q=FA>jXP6X8XpC0b>>3YzxS1)16 zxV)smj2V;T^ZGW=M_-n}GQQ=_I;5N~<&g*E&CVuM1uVsUTT8C$*LB2qx2L8u# z0hMXxn~)9OHwOf7Q{`X%5=4;KNPm^hJoHq2y)bkF;btDk$t)-SYD}Q91MI4TP0sgy z$K8eOUx(d;tsU}J*4z(9ff80Sh(URXff$U#g0*0gwVe@=`k;p>9 zAJ9p+sN%x7?wsFvckfSj@!g4Yc zT%mSLlE!)HcDwcQQ~9i=|LRAfjl!~qHrMOnN+zi~c9&6okF4?x!0^S$E!4t8ixpMG zr$_y7?{4P!ImZ-mILc}liyvYi@>l1=Z_~-I`dqRh$!}?G+r!wkp@6CwbJ~am zNTvRtFTdx%SS=OPSGq4Cc>DQA#sw{cj^@V#${^x7 zro4=G^|Cl7uDcr?#l5|V)8LKu^>=5{RTpPf$`fmBhJHTeEAN%k<3w@=a)iGc9Z;I8 zQ&@IA?sw#RRLKqF$$43D7sIsg1-YF+!hQ6Le1i{9A8}X6fgaZKsEF)uEb#9AMCHiD z#y18_ZnfYaYs}@mA8Ch8fDq1DHm83gS2`{qHFq>$`ngJOv}!U#zw$~dnMkoC z%t|n1Wj_UXZhfwzV#>Ej_znqZ>kQ-K!dlMplwc^9FA&Q5WJ3zxUQ7U*IU4h6oz3#r z2ilda9cCr=GGT*+j~aZDijK44^>Q0tzhBN^7jslh?{SzP{L?KXlQ-QSK6B^K4Px6S z0;F^zB>#<)^V)rg?xYIN3$rREt#BT(v=N*(PD;VwSpInj;Q7j3w+7V8L<4<@c|2ck+{C^b|4{Z< zKVxou-|3#?29~VV@GRwkb(DR?NZ%{5B7ZFX+e=lq&Ab)`4)Mkq3Eza?rZRyc{2izQ zU}hBrZ1sNDtbiM7ma#CO>Q;$rdbeUKNTz}$;8RwsV;KucW#;IHZqf({x(x<*@j%Y#8gK;2co+!)ceIwbA!4$ zvF+b{MYK$MdLCYch($AvKy`gwH{GRjyzjq?&8ztM#XL5TF8R*aFa7sl9*#@%AgcVI zjXh3?FSA!xAi4e&Cec2(j2h+RRt(S>=;QUmkW#e_9UQaD(clue0!=T z{;?-Iq@5o1XC{-!xT|^Mbs)(Yd&CCQ3g zxRy4!10X7V)gSjXtzYej&Cn$y8#xlhMR)Te)a_+^lu<*K9sB%VdNseKuuZ&$Jowzi zrfR6&3M=3sf2xEV9Zk3!p>+w1$${Gkvv~tzb9-0ui}9v@Hwr4v1`f19W-C<4q9Q;) z`}z0epSwhvaRK#bYrTqI-)Mprqs>Ez^*UIVHzqP^;rlQH+8vf=p*W-5l+Qrcs%kx$ zk~`Cy$;q5Il3+HyLP{l?_4^IJ5BX^0toCs2@bk+uhqCqjDgN`ny5r& zmcX3MZ)5MH5>M4#2<9RcqJtt=2`c>JRj+&&km{3ax$86riPXOXMEfEc}^WJ;?Q-2N2CnR_| z-M2zbDG5b@$y&dV*xjo)5k!I2U)rv8*)Wm!FpulZ-*dwVx)5)W{enRsl+ORU1?apP z{>^es74)SvME_L+v5yBAyB#_}R{WkVRC2;)(qcfwR6&G5iW0~#%abr_P}Ujo|%iTWCPGmweV&Rgv<4`q|*L4quSCQ?x8OP28Sl`MC;nlwpRybe#cy+ZS6`Geg*`{`vDZ+^x)jgt)|cxt zJ~0DVA~W(sk{!!);Q^(XO^)V_0q38J7i5fmQB3XL1z!3AA*6Det4cTaSJc_0rC-xF zN=9l=J?rVaHcXRaI->AUxj%NTxMQvEpy)HjCo_!H0(r2j;>-2KjGc)~a~Ndn0WQ)V zb!@h7eAwOlYTa&Dd@tREyU>IxbGnb|*LG$YM3|6Cu&m2}h>bZqZ_T~GN~DG}_483w z(uVVNEiOh5d^c=h79PC!b??1CSIzSi3(NZp-2dcbjXet)zL~vMBH9mkH!k#uZ*+{6d3&abU7ELB60g(8u?WT6*suv6 z7D>Nb9k20tNNh7Lva?pud7#Og`Anx4P&(bxS1{vmYt;ryMc0h?Hw8>v#Z*Z9YuVB5 zRJE7QpJym@sLSOMGOzX{KJezIdN^5QDh{O>Lc4047eD9xA|>k)ZUl@N297HbJ#BLP zv@Io=ErGivW?6sGEbz}n+I1|V13937;4mk=W9;NoBP_u!ew~>)GI+D^W8&#INoQc6 z5n^wPV=0^JW2Bbw9#rWrW0UbepJp36p2bN0Ry~;fu!jVEjt3x7Li=SMGHo%<>^Wpi zN3@10soR)|=r1Y6E~u5v2fvkM5^7-cOo1iB?lmtliS0~{QDZd;KG!z!`{H^3Gbv&S zZqc@#rwEHd z<5N5Klfo?>nL|*rp+)4|K$U30wkg>H2#E1^X(AsX9H)_zK!(bVWi&Q_Je&o7YzC`59*I8ps=3f`TbCIlMv>xC83_xZD@QSxVnG91pb9-r{v82mnJ>) zf)R=H+wnIbZ2c*xF4%fsb3mb=(AVGk$0i?}657MMoN+~KFqF&VN9(@OAu-7c{A5AVFG8(%E4EF(uOT5qanazwClb2{RS-lTEX$Wo$v~vzfm>XB zRSI|ck~Gf~Zb7P(gRVEK81W-*#s?C9afjD8b^aQ6(9}(;(8-aBUS;sF*KsVGZLlL& z->Z535zaBC7fA2x_o99DAs@u;1r@|>*Xwtj!NH4B#QWyHX>fWD*wc5oW}_!3$=z2| zftlFkLKc7V*XnlFu!i`)u*B7-mq9!p?(c){9l}fvkWRj}eYcLh>Kb$6usUt_cG3Et z#Ovz~FRl6!o;M$b7?mjE_X1~RtDix+KZ;olO!$ zCC&(7G!F_6o2o;5R~Jxh#J_knA00H(;XSPO_nG-xERgCVlmRA7mQ<5~2cISZK^@Go z*OLG`H*Wcpj0o0iX|a&0x*lGy!S_#A*&NZZ!cxPA0RW>nT`x)UvcfN7>jh^<<0p z7)VOyJ7&m!uMJfNjSyaZ$5Ew~Dkf#%VBc(BO8{;gG*o_z%Tt~s*Q!bg-Vb+rT$MeK zhzR<}3L{d>xuOymMFO=kCpXH%~MGD@>w(bpGVfS$X!HP+FV`hOH^B0{~E z-QNkk8UAOqWr4va9mZLHSl{XVF^_D$^3&HUAL?#z*br)T>lHa$W$ypQDXsqFlrx+}KRE1PdO+RzM{LwD#NGfjNziEt;u~Qz&%HsMXnNAFm?6>R6;11sG2yq*GO`9UM533XVPv9KlbRiJvs(mQATJ3^T^z zIn+*?~8L7~>AsHt%T&eX75t%D< z!VaVX)5?)K7Qt@Rj&2HWg?TNa!{)@>CpsAT#`EVp$x!(=sklBS7=!hS^C1QKqVGA6 zx$XOx)ps5hodh3at@yo)+C1|G5lul0qh8;BQaE;UIS+rmyOSwi8Qk_kprSVA?=VCk zKOer>ocz81ABc1fv5VO#33XenVwteWs%s^Xm^zt;_Rw%8WQ#J<>5WXWdUi&Um?wYM)10n*f|jURSa@ zyFd6^l^I3hQlI%MIJZrU9-J+Yl^v&PyIquvEq&g{gz{!0P6jDHYh@>Ti^~o ziy5~d(prDhn@ren8o`A5Z2@j|NJp0n9B0*zd6f*)@{J_QlZ@z?#||X_*A<`^(S~VPcNJ_eK>d3Fh2KkG4vrIcYP}6ka?*y$7CI5;NKtr% zLx{cTTJvA}M_7v)`rGry+lRB;v*Mf)UbnJf8c=+YEw0R0H*4aO2)j{an7_hfgFs%* zvLw`p#KxBYyQp-w2XZYxDip;&=f|p-o<|UwjvscRzptJ>b53+y@X;tsnr>V}fg3@l z5S9O;Z>Kv+q9uFS7|4`mo;Em37CUm+*o4R&$L{L_ouc|&#c8#Wzq>%p5PWW-V(*?c#+N%0QzRI8bV~eP8eH#wJR4LeSh0!=g$@)zM z6YN|ceyQq*1u*J+PN9oWBBaCVkUB6`X-w%RN>A=4s%LN_!Jt8-4Th|G3+MeqI2nu0 z8Zcf*z|d zjR~x(F}cSX+g%pcLRV#&A3W&HcNfsBBwLe5F+$pm|INb!r`iU0(O!@Fvr03>vZPo~ zlW+6w?dpLPzH1Bb;mzl=Zbskqq}Y>*rtVW{7TJAv`~K5qX`iZfP){hGnfAT-27ck4 zyE##_r}hIrJPX4zRI4@>1=qu^5nO?$D5YA52-b?dOw1J$*p#+nje6eBOBxCq23lZZ zFW@|;DTZ8Erp<_~(@p8??k20RI^X(jCB%b@xH^fpv4y2)k;PgO>Dgr|9iktz0+i*m!ZywIk@Gk- zTwU_`I^hfUAP4UGgXXR1n%+^Q89k_J9GWXIA>*)XdN@Le)~CDr$T?6-Z*$nQ0dCR0 zfNzCL2|bp%Ko%pMpz~w;sEKW^g!j#Gj~O2FRbCIvbsBK|_)-ID{R6M&m^Db@QrWQ? zL&T7?q}z>&$#izjarw&P~Jb5H+7v<^(-xE&XFu2?oO`M(TzA z3x4I0EyWB! z+}j%Yg_@8jppt$Yk`^oQ^9>eSX>DAig=g_^X9?iif2L0~we$;#!tWzg9R)zjSZSj^ zN*u*6J_-340phU`f(iosvV83bsJPr$l+)*)rDElK794W;wM{U^NZ|saRq_0~h2iHG zMwQjaZ#}EZ zjE)oW2)zbJfa+Z=xFt?sLbTr$rZrr7)cstN`{`QEC`fc(WAfVP;-??!MeL&z(bY*5tqAu^pR*xA0+s&Q^jsi=7q>!f+;SJLyA84PD z$RINoT*_LviztOWktNw$P72W zcipAyuOAY>)7d&;>*Zqe=Bt-{hW?!QUsTgluKfCB`J{`otU?7hZ2(D6=Id);1esj! z+We;0K%Z|3`R-1sh&Jn00(A<0pkvXBJUX9-x4Ll?*V2g{^~nAzOH@cPI&c_uN9JZ@ zt7h2$)6A76p*KZrK{bm%8#E{UE0q6>9$Zop&Z z%07u+PJ^`NgnvG`&1Cl83Ta8NK7A7l573IqHgBjA^$!}qgK-ch*>j>QCLjHExQ%_- z@X}Ynf9yvvIb06v6Ilf6o8P1U%Py)~#C*}m>)U(~^>0uqL?_kRTA#Qq&_gt*n$}5j z$2;ZqiVJ)7AQ^7>=Sp5uQmF9=xodxOePGH{c-^|s%)U~)?M@d;0s9bB`z`B)1qAKsQ;IUG z(JNvEOI*xho7?X=H%3Li^KT7*+~)({H^KD+)GBC}wQJsA_{cTk+IOJ3)Y5Mdl2E6G zAV?A~8>EUXIB9#3q4+`#!I^dyWsPzA+2@K$xBTuk1E3>zUapvyX)e3*cO1qw>96a* z=Gk;+Zzd#jK{L7_6>FkHVDm0tj)yDpomYK;E)_9oSiif|2=@@cDjgP_rcm4S&un*T z>FEH7%1w{j@6zTx;58xIy|VG~P7GStNos=$2;MA695p-TZxY+C^qiX&UIe# zJZA3!-^Yc>8`uP}baBY^uPK4yl%@S$lu|gaAm{#?oHq@h+%TsHDUThvnLI@gZvP4ZyWoAfkV1NK3s_VMB4H#-XC3P@8&C@4bOqyFXoi~)WeO`*&}m5hyPsCs}js^^uo`Nx3N;HU*b|+ z{AxNj?@5wwL>Dl*@!rpkgf9l9X(}7yMAvU$gZ4p_WSL*&*MlQS1e^}%ALu|E*^i&o zcm-;7>!!Nk_-T@kUy*@dB{+nyn>p77|HjbNz9p$9hj(CPFa{4aOo_46)=q)K{+s;C zMZl_35NjphiVC|ukEqDWfRV=>aaN!u_cmzLnr3fu$U{qjuGz5xZ35N*b>Sb$;v6S- zG34v{fe6C&u}LY0!^u+W!-}HdN8hnkZJqN!e0>JcOSxhZjXGi8g)N+*o#Ni_giJxp zXWIH+aV!hK_Zdsi7ok(<^5IhW$6a7VtH>8t{l5Ae=^62I{Czf73YO3!)okg4{^3ym zKZY?QPTIc*z&H?D|JTXJAAt_hb#q)=`&j&xpuQbAeg}F`d(v(FZRYyOqYsb2Cb<$j z(0a7?9Osp=F>jW5Qp`VN`tGYg6%$~lbv%vkQw}lumDa-FU9x^~cC(1bkCjdhE*Qpq zt2iz8x7H`rD+hLe;-{xZRJS=aybt&;$AQ*2pC{@6c64mVesW>vL1S4eI|$roW`8V- z#pOw;`@hKrMf~FMf3Brpls`Dj^V+A?NJ~zqOFGQ77W@5=HnjTWiQ(OW31&fVRKKwQ z_~qk{wepD|IaI#)!%~Qb^OIF2)15IzJbwggyP*AP5yR7)aSabzx^rWd3&yCEx-!qe zql(1srgNt+F`S0TKyl#v+Z@m_@d!r~F*ueSpjgm}42WWiSD$1$5q8Fg5sUP3dkJmW z6kX^#n?I{YrP}gs`=l>b?SgL&(q2`Pu@BET6&Wfd140fj*SYyYS~rYVc5f3|KCelR zeWogW!kRs%ImJ6%@j zP3v1gbDsYtW-_;Z7B6V}-9pnHk1qBveOlO4QID+sALsYe9UKjkU28&kb^e5x3!YubilrEUPP|?F!F^S zq^PLT%(V)^|LkPz+6ST9R-1CgZZ<9-b{ByDsk=xl#Y_dRWNUB3CM6V) zwJ{O`HZKdw`(TMSP)T)F&}+SSUtCc0REvdJw7m8aj_Z=R^?Nz_ak&s-r^5_=%k*p8 zy~i(bti>GFZ?Y-oIIlBH+nSAcu4Pv!f69@1!4`+~ET7TluMmPRcK3R%Y@2hXjqNZ} z2g+OXaC_VSWTb}U)8p{^P%f?}pRM;e1YM={@il@F>DPvMtW);}s96o8n~(Lw{LExp zbl2a`)rriC$!HVKW=z6j#pkvdsqZYbA6)|v)%830rHvwnok1YM-C6S`QY;pK%#!TV zXfM+!q9Vys>`!gf9x4u%ypRB_*LB*Xva{HgLvcF3vRB{0C#sX3uqjsDGsix*+|hWy z1@5=u9u0gEIhT~IzX9R|WAc0W4oD^Q?AUK)8FH-$T0)P(8!}ccvyf{T^X1wJ{nRT` zSJ)Gyg|!+`cP|&YS(<8_~gP_X*>T+^AU}MX=q( zsSwz~ByHXN-2u8q<;pOu6t@@T-$@A^=+=-hlbu^ zUT?uHtUk{w!_Jjmc`m+C5V zZTO@{xElikJrx#CPY0+}huc7Vf}ya)q7o zn2><@G5;7ijc$t0>pAu_$4b9Xb9%cMS@*dkN7LPw(z$z9e?jHXr*76K(8#j>MNAD5 z{oR;rbK4zgqi}aAcEwZ0LibcwLmfxkumvF0IVI@x+wkW z2~6^tIyx@CA(V?%4319avWaCC(S1q! zLG)RL#I)M+6!~0%jeg+!_FjM*8hHX1fK3f z60@l29<%>!dkHVB6f3aoSOia8UcuA6E5#P}79=RdW({+TegXlGsoRo4P))ckmG)PA zn);daGCC~FlMw%*+UBOaz_+!k%p2FzA+&avyq1R@k$aMv)mHw5+t#|kcWYJrV&jH@ zCASQ|OnPD@fr>~gsVY+?!OUDri(Rytin31a|3mRW5vI<88R@Wu@T8g^b^BEjEZ(5J&4UN?obO;-Ru*im~ zVyi6hESGM-grGK@o||)&qP20HLwr^V0GPv}Fm~$hmE{-wn3oLp5bfAB#8E`i1a|PU zOmDYY)cjj$kccfTltBO2!!Mu^DPgza4{ZXindSMgCYUIgC?u2{9XxP1u)OjALte>n z2L%oU)9Pl{1QUr36;p3->R#Xy>Cg0#M;anlAi8^9U)0VSo z#an#oLD$>ks{)KPpWQQT8`DQlJXulr-oPMYP0WrCe&dS)870XRlaD{gKY+;Z*5duw zv%2=^Useeo8UPR9|2;(YpX!geY2MI&1I48E=YGx4)Uyzy<~7Q*vV8{U>MaE5xGG$N z-&w!K%%vRx8t875bdS?jFXX5~a@2#1{Ctwnuu$W73Ybn^d$V4jz2~P2C@OAg%=q~K z)74joMb))Yi%2QWLrJJ0peRFkDIh7O5&}a=Dd12uG*S{ucXvt*B?1FTN%s&VG4z1L z&?)^L^m*U!{k|V_UGw9d9cS--pS{<;);hS`k!PMU@N*zwPq;7NaZMS;s~>IeaFO(@ zL3=IVvAr&N;PxFl7Jj(k;8^V&{Pnb=#WEtC2KoLX+Gi43w%;wg!={GdLwbDc%-&n7 zuW{<8Ot{ABuvGEzA8=MmVZLD@{BO(lmE9%>vR8}{9Wb>1oe0!drZNHzRe30@uLHfZNxxA0v|04VKh_~%{%AKVvsd0>jaN2a=^&UZzwOkGTHAjA zmMTB+n&yGlfVp=ap2Dm{#4XBnB@}NAiTqMz?cgWG=Wx1d#l|RQdVq$`NO`v$zlsKH z4)05|a`u<@yT1F|V89t@%A4PQb2OQ`d%rT@nK*%ebyT)c>hqO-n$q<0=NM$%`~vOr z)F4{c*n?Ah!YJ=<6|KchsmCLaUj_3#Zl8|b)!r!Iv#)7?XCXB;X45_(CTlRdJH+< zAxMG_h?qp=&)Ns*9Alo{(rvq!=vQD)(!|8gXGuH=i&%?yL$5`}GUGID?_6seJ8!#d zG|2Dh?sHZGz(0S-&p`)K$a2Qn$^6K#@6lV*$gwtzAC%U?{iDJuCBxTt!#H^Ru=_HQ zILvQ4QNi|eFYaD+cVGy}UQ+{=bO3U3gv}A@=#d$oOn~n*@Z%d!Mj0OVuSaI9s;Ba# zK#{7bhcTOh2x#)C7PC36){2nEg~t3Y64?+T@p5dDF+6(4m@1J7Srebfw9&5I)&kT8 zqz}bGyfCX`bSWF!YsR?oUZ_P37kDjCmLv%^U#WS6Fj{u5aL{=)kD8kuz5ERpl;T{t zv7P%%sJTtM8&=k)>EQg2Y$1~RvUkR81^qRP0&oENC57ON$9>~ZmS#u2_b3wIUz~h# zS6gvRx0d|Z|1l?6gDaU0H%uKetoBfi4f|+|0(*{Wc$we^B`=bVXJfK7_KNz_jpKdeT$BMeJQJ>68Jw8Cw~pPn;rAOlSW5U| zHINH3otRHyYOXwMedWy)@qR*T&8bk7wb{K$&+#!WP!~geYUi+fo{~93#SiyK3^vu{ zbeToh=JXSd>&wi{Wb-BI?%QqdyEz7G-Wz_7cx}$mCc_K@1Wx$p83IsNY%&x`Vw)69 zQh>+>s@$+|$?Z=bqFYs!6D7}wu|#FXk<6j{W){l~N~_}M02Sdt1K|Y3ib6?ay$&E? z)pF*-a{x*hKRa50CmMa@if`Jf<$V04ws|57lAgWdWb05=pD&TIP1M1uyj6im9#9}Z zpi}A!xH7ymM(GG|;1H@;qNBJjonwP|{a7a0o{I-Dwm2eXf^9A4ruLDR{D5INB3+Qj z8NKF`XLL0D_k?()OdBMh`d@MTMBrkQN98j$EO(8hHcQNvqbqcT@*^hw$=XOOz z{dt1q)F~xWx>sv&A_rb2N>#!J9`6j%cd)sQg1tO<@qnRaLo{IMuTalt!lFo?acYHk z2_K7Dj1XjVF${GyCtG$J-M~q5^R_i&`-?QYON9m}X?+13T8uZ6OngfUq$!NzX~8sE z8{1^!$c-FYGDqooWy0jYMMZVB+W4B?ZHj}0C3fW*5+=H)C3MdSu)wTK=t)5K6^ra^ zHHn)xKcx;9> z-)Y2^Ch6#Bu#Nb`N`R~(-zCOd4oIlQc%thqMVuesNvx$Zd9jh-kOz#UTp^U`L;^#?&)2ahl=RyM#gtx=lPPkGTxdrgff+$6r21NQ^VL=AVAT2DKF4CyrqO<6_iZmXz>oit+PVU@44X7w7VhcBk2?EeV$Mi}C2~etwbZW15Nf z8NZHF%QivP(xI4!K!jee=5!qtwZU0K)LXh;;%s z(#oprsx+D@jZiYp;=@p#OOPG6Dvq_g<-4K)z^?m6`uN?p+ZGu1Hz`NdW*{hu5}VXbIcX z4itcRvG2fiSW4#2#|$ty)?#>S|B*%o{8qtBuxRwQ%Vxq4zYLD(m`$v-+I96*Yn8nG z>3VCGvH!yP@rGIxEl&w?UDk8LATiGVuJ(1QPNRRysS#j>eVs|^M3Cq_>fv4 znTR{#LNok|&ZSaH&*ac%fDyoHsM14RE%KC{R;Vl8o=e`I+kb|@Pde}l7WA-kPLuS@ zw7g$6I}TgQQi@QVDXVWdnYg-%GTlz@_Mm!hbAcctE?Uwrt}2M)C?76KC~LGXS^Q>= z_7L9j!CHb}Eyoet`84uYlIG1M26P@%ID#4VDMUktEQ#u|chb;ccRKGE##Ws?!+iKO zQp-6^R5UCyxS31redw`1VCUhr0gRI2D*uDmSxRsHcWxX>@y@h52Bie#yDB@g17w!D{J!U`xVnNPwNF^1@d2 zrw`vY+92xLjnJ@zOXoDg{|%v8eSpBJgS)EV{2(eVX6pe1?3~hQSDTnpg3ws){0$ zDiNSX36h|@KO{5Qu~f$9Lsdn$5|;tBq*9!tyM{ATrwtSBWt72}c}|i8G7Wc&_v0en zecsXI@iJ-VRWzi(G)&Z{qmuYnG@lS#;aplUiTuPYhn3}rv1p#hbJSmQIX4TMx>kMe z=YZ7lb|liCp&7w6>I|xNj2X2~!RL=g=a$>r9u-j6gF=qOBvH)>2MW#mb_wER_%EQ zurRXVxFwDU>Va8;-WU6_$Wa7@*kk$a4*~7kTo^bAcb~zS9R@a9uIgB?e`QM4%rWu; zO#j05i%4gD}>_znM|T!pWm+vE)}w#2}x|PIslLD{Nb{rZ+}Ku0>LYS5t{({G?wr`$2C$pg^v9) zTtZbqBx-|ME)U?OA5R%M_M&VaqSqv)sLxkm$L}S$z>2mXHeT!AkujuaI@XGOsKtp&kMQ*cZt#tiXmz;1vSa%M-oZjH3-oQ;YwP z(i`^&xJNTz{v~1Ibir50$miVPV_2@&28D2*=giH3m%pjLJsJ;n#U%WmTTW0a=3`by z01`-YFSSLJH`5=+Vey8GmobJe={Cli2yPf!3}nB8xRupiHtD{}A<26Mp&2S6A~gmz z#`8i=#hLLQSgG~@_`+&l4P zN&q1$>Eo^^G#+ZAqH7s`OiaN?kCj#io;07wfVwa{@8Z|FIR3o~MFcbgjd!b}@v5Me z@>5Qr2Se28IRI5}jGRRdgd=38Ng zTm_22x(2xw{lmjAoxe+}E(MH?XRJb67Nq;4-32498t`Uciv%el8GW1=T2uGBx!P(RI5oWRPH$Yg(PF7T)M4HCEz-7s7T4aN zZAhjC0&CxNUP*LQ^<|c{`K`2HUG1G&ER+Y*->OaTo#oC@UnYwA1NJsNH~#c;$RHb< zIom|x0ZV~deS&DXxoCqC$z)Gsl;-UV$Ynp6!TSISaK=BTH5s$SAB4-(a4x?T9X`Nq zd9z>KcLE3X=W0MWdyMq9m|EW%&GwVD*hn`%H#e=9F)C2$o?uCCseI0v!#Z0}XkvnAWhWy84=miv;1_Ij@CT4Qai+qg+Fr07M zLx}e*gpM}oI!-O{sSP+MkF5G5+ z>(L*hCEZhGP80(f9=0;b|9rrSicjr<;egB`7FaI!RyX?keN!duls+SZMgDGO{>OI* z_X2L1wh06r0UoNk7v)B^PB&+QpcMUP#K2oV+Q$he6nt-dh-?f-(QWqN5jrc|(?N|b zYV*@<%zsJC(g3k%jm zSOk$}Na`J)14OE>JfeZb&;7TTeoc}WVZh}F$Ebu=8A-8yDj6lkb#w0+kSkT4ul)kH zg_J|d)MXvt*+3A>tpH|LjjXQeJS8LDR**eN*ksJaO6kA9V`Q2;qgi;tBe$nDW(RR- z&SR?k{lwGZo5$0*nJ9>m!AymHl`&E$1n~e(b+&vmIPTIwa5rZYiQW;H+78*pi{~%J zjCwt(!9P4$PlP6i@{!m})#{R``iVBQlgQW|mI612Bd$#_%zxxO4Ur!J2nN(+AMLA! z4RKYBHZ*UJWlr_eR|`~zR)cRN%k<4c>cot4=7h$+8;sCXcww|NQi1SVterSKXu%BU zZ&)5FYom+>u@_O%!SHMtC^yhW4R(jrg zI*4&c_nBuW`CqL(r_e=KfdVDY~A?Xv+px90f)w8;VJa=TS0z$P; zYJp9IfO8@|)Q`To;gGh;hCP)R&YYO2=Fu1>s zVe`p_dh&f5LPah1=ps+q11{16zoEXCKi&~?6)gn|@m=DVyo1_(2C2<0#oA_Kk!e~^ zTY(+?4piUUXrSKsa+t^KXx16;#`yInQ@=OndAIk1osi+{TH+{-&Y9l&ZqSj(E03pea;iZ<(vwja+~9Duh|XV{j9dy}*h!17I}^tEYoMoWSIcs1mUPRv07B zT5v%P28wX>Pr2eZ4BydmA(M9-+enADK0C~3E+cWmdp7p0P&gJ24A~83F zsmpPjhmY6VdOs+F7uZ~>H|1MgkB));+_FMHGy%j_L^a_o(wr({;X1l4@rq-Oj#Svy z!uHZY9zhF`$>}E(kp(Er5PpAv3~T+K#k}a4BL1}AXsqUJzW)Gt0~>!t=IcGQ7k?VA z0CE5e37kWOgw@15irG@`lr9Np}$UMsG&R@d|I0;w1KUVQ>|Pzf(k}O)SVaV z+amm!jD~PLHcZ0hE|###t`k?RzvYCs8|6Wj0d&tYL)f3>;a^z&kNpky8)W=nCb(hk z8d4aLDeh*kAk`x)4k?i}pO4!5wWX~vEFzS!cnagzN^Tje_s(uY(DY(^2uAn2uG&ZC z0~Xd#HW!F95vP0$!xG|XiN(kIxbBJDSOFk_c5Vhhar!kzwRa8Iyx~lmWS-GDIGgqYEV9q+ZVE z_=tUl2HHCsuu4%SJ7h6qFS5FFikQr=&9*2HJ{@aKrNvc($b^UDAr?jAk^lpW|LrVw zVJW&r*}1aOR(=4doOnt%yvu63UhK0aUsz@B?h5cHw8t>U_wXG;NL$o@w2rKO4$S*o zMio3Y3UWOq*wzfxpl?hGV|n4eK|z!VI^t36%iXv(y1yRwXyS*#%5(Y7Xf6?(6LBBpH2fGKyy5FaUt zQN|5ddRtO1ruh6}hWupnud(b}4bZ@Qz1I}55y{ZxfJY+nx}WqHsnLg- z`u!#7(*Pt{q=^=@{G4w>u!){xAAS|+oVG^y*XF)RA7e>(@{923vmX?gv~Xq+4r0)K zE$lw&OzX&^%~O`Npe7IF!Eq4>_nDGkRMZ=Qt>(v-5eK2#b|ZSz5YhsoiChFQ;7|Jk z?u-i3g?zezffsbxP#AK9)D8lPA|G(#)4eB%fQmMUZ~yQ~OAR&s_|*P#z*aT!Bq0ye zs!1v}6|7AOB!RfsIvISm%zgr%vwx{9$yT~_cKtVe`s@0$i77>;@(p?+N*g*)atxsr`Q(|ZT?$O zZCe63Z52Esb&11c2xF!pNA?M=Gy;4jTLAcu0W^4{x62IyM}t3crTV_bMCQ{-)dEyj z7hw)u5n!bn_^csR3*5na<6ld02{J=X7Wuj~UUBPVIdBacuU-S&ddQ;!OSmCAA`Zto z^;i7I`lZd`?f)fTK_5f}{&>riZ-DJC!|aL&TeMYTG=c+O?F4T{G+Ia7ja>dfS%n<5$&J zwpuHPVVX}7TeYbjE}?Kii_buYN(nu=_AIE6h9YY#6iUPg+mFZ=u2^+AjP$vQ32jF$p=!0;s=Lq7YK1>l&zuFuX3{!gX|gr1FS0xX2j5onq^`ZshU}nE0T)=#s>jan>H@oPHP4!?w(V9L zTEaS2Cc>hw*l19@H?Y|L1&x*O5u)=_fb`=8%=>W!S)KDC)T1V;r8hQuQ&}DhAjW{$ zW-AA@25W^6ue8<2sxesoD*BV`yl%ug0AnO+9!EZ^vTS}nN%8u@QL3n5YA=w~0;-M2YTW0X7N{E5hu-eG9E2#l0V;N)J%N;}l zVC_|Hu+(J}j(Q)^uB8K@D8MSMr%g%b7_r1I?g{LoR9(*k;7DmbbYkxB`Ymp?TEF^t zy|W5afX`)_Ijuw~6LRBfSqy}k0=@@-b{1EMpg*titEFi5A7vaMKf5yB`LnX@>LK|Zuae&kKx(;&8%_V|_|c@n(YCGZc*3Syc8?|xQVf1| z^>$)0B&ey_BSwcg6q5VPycGN+G!o`1+%9suW*bc&)$GD%*k^6Z(zVRRI`uZYu04E_ zJh)i$u^E2tnTWM3rq(PVwfA*N(@X{w-7}KLUDcUif{;(oe}t0z_Kq9vU-Sq(U3I{{ zHpLQwuy0&BC+k-Le-dVe#k=aaX^5=#Hmh;26(@f}OimO-GxH>a-<@~qPgi|mozdk< z7Wv*zt&6RHM>{Bg@v+J7C5wHkre=1rW)EN#(Ox)hYzCBvjs1XH2vHsnR=b@hHCqd?Q=Pj(J-H&@`cX zz@ZAiXKGptb2@+TxKDqs0bT=N#F%rxj|LJ^c$-G;x9RrVec}eXtdsI0&Du?7?Qaok zFdEpBud7Ar`b52L&kD$Adi$le{Nl6CZ*O?e!;tLnUsj+Ct@w9>c9GSM9Ng;_VLESn z=cU>;Ud(a6d`PevqrHd3___}E4_-B(5oo}hptS6n_ZOJ zU|5CYThFV{XJc=_XPuB~*A_KnX>LxiQOVf9cS$8+wSq7mT+SoIGb;Ay7auFRhHaUA zZI?SJ+Wfb6wQeoITl{79X37vuoFlbqb`}1N&1-h*myDbt@yzR$Z>$$sjf!;KQVoxvEFx9#tGT3z8x;(bBRXkf>8c0WF~Dg!zm z5dKpd_RSqD(Bhdew^O655xBsyD+=Mon?m(|;al&kk-w$UM^0W$>>Wj<;;v>PqAI8@ z>JqbdHG{5(KaN=B%@LE^T?Lv8MV_JwvMmc+LH5EoVBkrZZdf6|s&S-&B;0-lNGT`4 zaUw}L-FLkPF!b~MZp{rx3G|F8ro3XFwa1RlFqKhUXRF}#WTIQ0iWmGy|`~2g1@?4bD=5Cpa`shcbBN^lb|A zNZw>k0_^7Cm=O5%#T2!Ui+)?N5Jv_VI!IQi;(0FSL4FHbBrkNyuhWyMD%yQ^YQ3AW zG#=2!py?)F6uP{(cv`GMeNXZl2niT^$!cu_E-m%>A}yMz*L=I4_2ao-DUx^%K5Qk~L^Jf}s5*WY3t$(&4T zxC`Exk^W@N;(dNU?Z@}8b7wq@%gvXnOJg6;FwIzhl5^SP{ycE`|H?F zbWIWDa`&k2@4ey1=V!ooYTBAJ^_3h=9}az%*}~Rt+b^#Q52Q`xn($j+=ejg(lpSsF z8CEoD21&I!q2n5I=kDI`KoelsR)L4Yt-i%(8CTvcVfWCP`g`SB*J_va@4kcpSb_t6IYeKTs_qtYd_oXW~OcB0mGZ&KfNmd(~SOv~;xb`Oel;-5)VIMX|Zfqpv<)tnd|5 zYZ#n}F;4jwI7^|z>-+Lf-Km`c+^B43!?XMYVB^w0=WjdWV-2H{JUl1#U5-JH|> zN_C-e(MS(NCpECGb7O7~YXn(4~eT+X?1}z`4u|JO!?W^P#I@65W5}2?%95 zsGRAga&EoZBt(w)phtg=50$Rbd&y3YKPVxD6PYwOW@}uYjh~;mlOMhQu8@}3gG}W$ zm=Os3NqmMN>jG-bXGQGKB{ld#Yn*@yHh$QAAwsg%pmlUu@lSeMl%C1L=IH!DV&+?N z3-1dj78vs445LlK#eTRjpf~?Bq=?BCbwXO3{bn=hs{WU_@tXy!N$g9TpBlswnD~v+ z3L3|w`!gS>1NY?)Uac$4NRx@a`^MWznc<%v(!jN_5;dSDMJq|=imN-m#Mfdt>F~bP T^94Qwbxq}&h5}s9)c^kg&(`5F literal 0 HcmV?d00001 diff --git a/man/calculateSIRSpace.Rd b/man/calculateSIRSpace.Rd new file mode 100644 index 0000000..12437b0 --- /dev/null +++ b/man/calculateSIRSpace.Rd @@ -0,0 +1,104 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calculateSIRSpace.R, R/plot.calculateSIRSpace.R +\name{calculateSIRSpace} +\alias{calculateSIRSpace} +\alias{plot.calculateSIRSpace} +\title{Calculate SIR Space for Different Cell Types} +\usage{ +calculateSIRSpace( + query_data, + reference_data, + query_cell_type_col, + ref_cell_type_col, + cell_types = NULL, + multiple_cond_means = TRUE, + assay_name = "logcounts", + cumulative_variance_threshold = 0.7, + n_neighbor = 1 +) + +\method{plot}{calculateSIRSpace}( + x, + plot_type = c("scatterplot", "boxplot", "varplot"), + sir_subset = NULL, + n_top_vars = 5, + ... +) +} +\arguments{ +\item{query_data}{A \code{\linkS4class{SingleCellExperiment}} object containing the numeric expression matrix for the query cells.} + +\item{reference_data}{A \code{\linkS4class{SingleCellExperiment}} object containing the numeric expression matrix for the reference cells.} + +\item{query_cell_type_col}{A character string specifying the column name in the \code{colData} of \code{query_data} that identifies the cell types.} + +\item{ref_cell_type_col}{A character string specifying the column name in the \code{colData} of \code{reference_data} that identifies the cell types.} + +\item{cell_types}{A character vector specifying the cell types to include in the analysis. If NULL, all common cell types between the query and reference data will be used.} + +\item{multiple_cond_means}{Logical. Whether to compute conditional means for multiple conditions in the reference dataset. Default is TRUE.} + +\item{assay_name}{A character string specifying the name of the assay on which to perform computations. Default is "logcounts".} + +\item{cumulative_variance_threshold}{A numeric value specifying the cumulative variance threshold for selecting principal components. Default is 0.7.} + +\item{n_neighbor}{A numeric value specifying the number of neighbors for computing the SIR space. Default is 1.} + +\item{x}{An object of class \code{calculateSIRSpace}, which contains projected data on the discriminant space. +Each element should include 'ref_proj' and 'query_proj' data frames representing reference and query projections.} + +\item{plot_type}{A character string indicating the type of plot to generate. Options are "scatterplot", "boxplot", or "varplot". Default is "scatterplot".} + +\item{sir_subset}{A numeric vector specifying which discriminant axes (SIR components) to include in the plot. Default is the first 5 axes.} + +\item{n_top_vars}{Number of top contributing variables to display in varplot. Default is 5.} + +\item{...}{Additional arguments to be passed to the plotting functions.} +} +\value{ +A list containing the SIR projections, rotation matrix, and percentage of variance explained for the given cell types. + +A \code{ggplot} object representing the chosen visualization (scatterplot, boxplot, or varplot) of the projected data. +} +\description{ +This function calculates the SIR space projections for different cell types in the query and reference datasets. + +The S3 plot method visualizes the projected reference and query data on discriminant spaces using either a scatterplot, boxplot, or varplot. +} +\details{ +The function projects the query dataset onto the SIR space of the reference dataset based on shared cell types. +It computes conditional means for the reference dataset, extracts the SVD components, and performs the projection +of both the query and reference data. It uses the `projectSIR` function to perform the actual projection and +allows the user to specify particular cell types for analysis. + +- **Scatterplot**: Displays projected data points, with colors used to differentiate between cell types and datasets. +- **Boxplot**: Shows the distribution of projected data values for each cell type, separated by datasets. +- **Varplot**: Highlights the top contributing variables for each discriminant axis, differentiating between positive and negative loadings. +} +\examples{ +# Load data +data("reference_data") +data("query_data") + +# Compute important variables for all pairwise cell comparisons +sir_output <- calculateSIRSpace(reference_data = reference_data, + query_data = query_data, + query_cell_type_col = "SingleR_annotation", + ref_cell_type_col = "expert_annotation", + multiple_cond_means = TRUE, + cumulative_variance_threshold = 0.9, + n_neighbor = 1) + +# Generate boxplot of SIR projections +plot(sir_output, plot_type = "boxplot", sir_subset = 1:6) + +} +\seealso{ +\code{\link{plot.calculateSIRSpace}} + +\code{\link{calculateSIRSpace}} +} +\author{ +Anthony Christidis, \email{anthony-alexander_christidis@hms.harvard.edu} +} +\keyword{internal} diff --git a/man/conditionalMeans.Rd b/man/conditionalMeans.Rd new file mode 100644 index 0000000..77782f6 --- /dev/null +++ b/man/conditionalMeans.Rd @@ -0,0 +1,56 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/projectSIR.R +\name{conditionalMeans} +\alias{conditionalMeans} +\title{Compute Conditional Means for Cell Types} +\usage{ +conditionalMeans( + reference_data, + ref_cell_type_col, + cell_types, + multiple_cond_means = FALSE, + assay_name = "logcounts", + cumulative_variance_threshold = 0.7, + n_neighbor = 1 +) +} +\arguments{ +\item{reference_data}{A \code{SingleCellExperiment} object containing the reference data, where rows +represent genes and columns represent cells.} + +\item{ref_cell_type_col}{A character string specifying the column in \code{colData(reference_data)} that contains the cell type labels.} + +\item{cell_types}{A character vector of cell types for which to compute conditional means.} + +\item{multiple_cond_means}{A logical value indicating whether to compute multiple conditional means per cell type. +Defaults to \code{FALSE}.} + +\item{assay_name}{A character string specifying the name of the assay to use for the computation. Defaults to \code{"logcounts"}.} + +\item{cumulative_variance_threshold}{A numeric value between 0 and 1 specifying the variance threshold +for PCA when computing multiple conditional means. Defaults to \code{0.7}.} + +\item{n_neighbor}{An integer specifying the number of nearest neighbors for clustering when computing multiple conditional means. +Defaults to \code{1}.} +} +\value{ +A numeric matrix with the conditional means for each cell type. If \code{multiple_cond_means = TRUE}, the matrix +will contain multiple rows for each cell type, representing the different conditional means computed via clustering. +} +\description{ +This function computes conditional means for each cell type in the reference data. It can compute +either a single conditional mean per cell type or multiple conditional means, depending on the +specified settings. Principal component analysis (PCA) is used for dimensionality reduction before +clustering when computing multiple conditional means. +} +\details{ +The function offers two modes of operation: +- **Single conditional mean per cell type**: For each cell type, it computes the mean expression + across all observations. +- **Multiple conditional means per cell type**: For each cell type, the function performs PCA to + reduce dimensionality, followed by clustering to compute multiple conditional means. +} +\author{ +Anthony Christidis, \email{anthony-alexander_christidis@hms.harvard.edu} +} +\keyword{internal} diff --git a/man/projectPCA.Rd b/man/projectPCA.Rd index 7775c88..de0cc12 100644 --- a/man/projectPCA.Rd +++ b/man/projectPCA.Rd @@ -14,16 +14,16 @@ projectPCA( ) } \arguments{ -\item{query_data}{A \code{\linkS4class{SingleCellExperiment}} object containing numeric expression matrix +\item{query_data}{A \code{\linkS4class{SingleCellExperiment}} object containing numeric expression matrix for the query cells.} -\item{reference_data}{A \code{\linkS4class{SingleCellExperiment}} object containing numeric expression matrix +\item{reference_data}{A \code{\linkS4class{SingleCellExperiment}} object containing numeric expression matrix for the reference cells.} -\item{query_cell_type_col}{character. The column name in the \code{colData} of \code{query_data} +\item{query_cell_type_col}{character. The column name in the \code{colData} of \code{query_data} that identifies the cell types.} -\item{ref_cell_type_col}{character. The column name in the \code{colData} of \code{reference_data} +\item{ref_cell_type_col}{character. The column name in the \code{colData} of \code{reference_data} that identifies the cell types.} \item{pc_subset}{A numeric vector specifying the subset of principal components (PCs) to compare. Default is 1:10.} @@ -34,14 +34,14 @@ that identifies the cell types.} A \code{data.frame} containing the projected data in rows (reference and query data combined). } \description{ -This function projects a query singleCellExperiment object onto the PCA space of a reference -singleCellExperiment object. The PCA analysis on the reference data is assumed to be pre-computed +This function projects a query singleCellExperiment object onto the PCA space of a reference +singleCellExperiment object. The PCA analysis on the reference data is assumed to be pre-computed and stored within the object. } \details{ -This function assumes that the "PCA" element exists within the \code{reducedDims} of the reference data -(obtained using \code{reducedDim(reference_data)}) and that the genes used for PCA are present in both -the reference and query data. It performs centering and scaling of the query data based on the reference +This function assumes that the "PCA" element exists within the \code{reducedDims} of the reference data +(obtained using \code{reducedDim(reference_data)}) and that the genes used for PCA are present in both +the reference and query data. It performs centering and scaling of the query data based on the reference data before projection. } \examples{ @@ -50,10 +50,10 @@ data("reference_data") data("query_data") # Project the query data onto PCA space of reference -pca_output <- projectPCA(query_data = query_data, +pca_output <- projectPCA(query_data = query_data, reference_data = reference_data, - query_cell_type_col = "SingleR_annotation", - ref_cell_type_col = "expert_annotation", + query_cell_type_col = "SingleR_annotation", + ref_cell_type_col = "expert_annotation", pc_subset = 1:10) } diff --git a/man/projectSIR.Rd b/man/projectSIR.Rd new file mode 100644 index 0000000..81ec597 --- /dev/null +++ b/man/projectSIR.Rd @@ -0,0 +1,76 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/projectSIR.R +\name{projectSIR} +\alias{projectSIR} +\title{Project Query Data Onto SIR Space of Reference Data} +\usage{ +projectSIR( + query_data, + reference_data, + query_cell_type_col, + ref_cell_type_col, + cell_types = NULL, + multiple_cond_means = TRUE, + assay_name = "logcounts", + cumulative_variance_threshold = 0.7, + n_neighbor = 1 +) +} +\arguments{ +\item{query_data}{A \code{\linkS4class{SingleCellExperiment}} object containing numeric expression matrix for the query cells.} + +\item{reference_data}{A \code{\linkS4class{SingleCellExperiment}} object containing numeric expression matrix for the reference cells.} + +\item{query_cell_type_col}{A character string specifying the column in the \code{colData} of \code{query_data} +that identifies the cell types.} + +\item{ref_cell_type_col}{A character string specifying the column in the \code{colData} of \code{reference_data} +that identifies the cell types.} + +\item{cell_types}{A character vector of cell types for which to compute conditional means in the reference data.} + +\item{multiple_cond_means}{A logical value indicating whether to compute multiple conditional means per cell type +(through PCA and clustering). Defaults to \code{TRUE}.} + +\item{assay_name}{A character string specifying the assay name on which to perform computations. Defaults to \code{"logcounts"}.} + +\item{cumulative_variance_threshold}{A numeric value between 0 and 1 specifying the variance threshold for PCA +when computing multiple conditional means. Defaults to \code{0.7}.} + +\item{n_neighbor}{An integer specifying the number of nearest neighbors for clustering when computing multiple +conditional means. Defaults to \code{1}.} +} +\value{ +A list containing: +\item{cond_means}{A matrix of the conditional means computed for the reference data.} +\item{rotation_mat}{The rotation matrix obtained from the SVD of the conditional means.} +\item{sir_projections}{A \code{data.frame} containing the SIR projections for both the reference and query datasets.} +\item{percent_var}{The percentage of variance explained by each component of the SIR projection.} +} +\description{ +This function projects a query \code{SingleCellExperiment} object onto the SIR (supervised independent +component) space of a reference \code{SingleCellExperiment} object. The SVD of the reference data is +computed on conditional means per cell type, and the query data is projected based on these reference +components. +} +\details{ +The genes used for the projection (SVD) must be present in both the reference and query datasets. +The function first computes conditional means for each cell type in the reference data, then performs +SVD on these conditional means to obtain the rotation matrix used for projecting both the reference +and query datasets. The query data is centered and scaled based on the reference data. +} +\examples{ +# Load data +data("reference_data") +data("query_data") + +# Project the query data onto SIR space of reference +sir_output <- projectSIR(query_data = query_data, + reference_data = reference_data, + query_cell_type_col = "SingleR_annotation", + ref_cell_type_col = "expert_annotation") + +} +\author{ +Anthony Christidis, \email{anthony-alexander_christidis@hms.harvard.edu} +} diff --git a/tests/testthat/test-calculateSIRSpace.R b/tests/testthat/test-calculateSIRSpace.R new file mode 100644 index 0000000..390ad92 --- /dev/null +++ b/tests/testthat/test-calculateSIRSpace.R @@ -0,0 +1,89 @@ +# Load necessary libraries +library(testthat) +library(scDiagnostics) + +# Load example datasets +data("reference_data") +data("query_data") + +test_that("calculateSIRSpace returns correct structure", { + # Run the calculateSIRSpace function + result <- calculateSIRSpace( + query_data = query_data, + reference_data = reference_data, + query_cell_type_col = "SingleR_annotation", + ref_cell_type_col = "expert_annotation" + ) + + # Check that result is a list and contains expected components + expect_type(result, "list") + expect_named(result, c("cond_means", "rotation_mat", "sir_projections", "percent_var")) + + # Check that the class includes "calculateSIRSpace" + expect_s3_class(result, "calculateSIRSpace") +}) + +test_that("calculateSIRSpace uses the correct cumulative variance threshold", { + # Run the function with a specific cumulative variance threshold + result <- calculateSIRSpace( + query_data = query_data, + reference_data = reference_data, + query_cell_type_col = "SingleR_annotation", + ref_cell_type_col = "expert_annotation", + cumulative_variance_threshold = 0.9 + ) + + # Ensure the cumulative variance threshold is respected + expect_true(sum(result$percent_var) >= 0.9) +}) + +test_that("calculateSIRSpace works with specific cell types", { + # Test with a specified subset of cell types + cell_types <- c("CD4", "CD8") + + result <- calculateSIRSpace( + query_data = query_data, + reference_data = reference_data, + query_cell_type_col = "SingleR_annotation", + ref_cell_type_col = "expert_annotation", + cell_types = cell_types + ) + + # Check that the conditional means are computed only for the specified cell types + expect_true(all(rownames(result$cond_means) %in% cell_types)) +}) + +test_that("calculateSIRSpace handles multiple conditional means", { + # Test with multiple_cond_means set to FALSE + result <- calculateSIRSpace( + query_data = query_data, + reference_data = reference_data, + query_cell_type_col = "SingleR_annotation", + ref_cell_type_col = "expert_annotation", + multiple_cond_means = FALSE + ) + + # Ensure that the conditional means are calculated based on single condition + expect_true(is.matrix(result$cond_means)) # Expect conditional means to be a matrix +}) + +test_that("calculateSIRSpace handles n_neighbor parameter", { + # Test with n_neighbor = 3 + result <- calculateSIRSpace( + query_data = query_data, + reference_data = reference_data, + query_cell_type_col = "SingleR_annotation", + ref_cell_type_col = "expert_annotation", + n_neighbor = 3 + ) + + # No direct way to test n_neighbor impact, but ensure result structure is correct + expect_type(result$sir_projections, "list") +}) + +test_that("calculateSIRSpace throws error for missing arguments", { + # Test if function throws an error for missing required arguments + expect_error(calculateSIRSpace(query_data = query_data), "argument \"reference_data\" is missing") + expect_error(calculateSIRSpace(reference_data = reference_data), "argument \"query_data\" is missing") + expect_error(calculateSIRSpace(query_data = query_data, reference_data = reference_data), "argument \"query_cell_type_col\" is missing") +}) diff --git a/tests/testthat/test-plot.calculateSIRSpace.R b/tests/testthat/test-plot.calculateSIRSpace.R new file mode 100644 index 0000000..3717aa9 --- /dev/null +++ b/tests/testthat/test-plot.calculateSIRSpace.R @@ -0,0 +1,29 @@ +# Load necessary libraries +library(testthat) +library(scDiagnostics) + +# Load example datasets +data("reference_data") +data("query_data") + +# Compute important variables for all pairwise cell comparisons +sir_output <- calculateSIRSpace(reference_data = reference_data, + query_data = query_data, + query_cell_type_col = "SingleR_annotation", + ref_cell_type_col = "expert_annotation", + multiple_cond_means = TRUE, + cumulative_variance_threshold = 0.9, + n_neighbor = 1) + +test_that("plot.calculateSIRSpace generates plots correctly", { + # Generate plot using the function + p1 <- plot(sir_output, plot_type = "scatterplot", sir_subset = 1:6) + p2 <- plot(sir_output, plot_type = "boxplot", sir_subset = 1:6) + p3 <- plot(sir_output, plot_type = "varplot", sir_subset = 1:6) + + # Check if output is a ggplot object + expect_s3_class(p1, "ggplot") + expect_s3_class(p2, "ggplot") + expect_s3_class(p3, "ggplot") +}) + diff --git a/tests/testthat/test-projectSIR.R b/tests/testthat/test-projectSIR.R new file mode 100644 index 0000000..eb6195b --- /dev/null +++ b/tests/testthat/test-projectSIR.R @@ -0,0 +1,47 @@ +# Load necessary libraries +library(testthat) +library(scDiagnostics) + +# Load example datasets +data("reference_data") +data("query_data") + +test_that("projectSIR returns expected output structure", { + # Test that the function returns a list with the correct components + result <- projectSIR( + query_data = query_data, + reference_data = reference_data, + query_cell_type_col = "SingleR_annotation", + ref_cell_type_col = "expert_annotation" + ) + + expect_type(result, "list") + expect_named(result, c("cond_means", "rotation_mat", "sir_projections", "percent_var")) +}) + +test_that("projectSIR handles genes not found in query data", { + # Modify query_data to exclude certain genes + modified_query_data <- query_data[-1, ] # Remove the first gene + + expect_error( + projectSIR( + query_data = modified_query_data, + reference_data = reference_data, + query_cell_type_col = "SingleR_annotation", + ref_cell_type_col = "expert_annotation" + ), + "Genes in reference SVD are not found in query data." + ) +}) + +test_that("projectSIR returns projections with correct number of components", { + result <- projectSIR( + query_data = query_data, + reference_data = reference_data, + query_cell_type_col = "SingleR_annotation", + ref_cell_type_col = "expert_annotation" + ) + + expect_equal(ncol(result$sir_projections), length(result$percent_var) + 2) +}) + diff --git a/tests/testthat/test-regressPC.R b/tests/testthat/test-regressPC.R index 296544a..1416079 100644 --- a/tests/testthat/test-regressPC.R +++ b/tests/testthat/test-regressPC.R @@ -14,20 +14,20 @@ test_that("regressPC works correctly with only reference data", { cell_types = c("CD4", "CD8", "B_and_plasma", "Myeloid"), pc_subset = 1:5 ) - + # Check if the output is a list expect_true(is.list(regress_res)) - + # Check if the list contains the correct elements expect_true("regression_summaries" %in% names(regress_res)) expect_true("r_squared" %in% names(regress_res)) expect_true("var_contributions" %in% names(regress_res)) expect_true("total_variance_explained" %in% names(regress_res)) - + # Check if R-squared values are numeric and of correct length expect_true(is.numeric(regress_res$r_squared)) expect_equal(length(regress_res$r_squared), 5) - + # Check if total variance explained is a numeric value expect_true(is.numeric(regress_res$total_variance_explained)) }) @@ -42,18 +42,18 @@ test_that("regressPC works correctly with reference and query data", { cell_types = c("CD4", "CD8", "B_and_plasma", "Myeloid"), pc_subset = 1:5 ) - + # Check if the output is a list expect_true(is.list(regress_res)) - + # Check if the list contains the correct elements for (cell_type in c("CD4", "CD8", "B_and_plasma", "Myeloid")) { expect_true(cell_type %in% names(regress_res)) } - + expect_true("indep_var" %in% names(regress_res)) expect_equal(regress_res$indep_var, "dataset") - + # Check if each cell type has regression summaries for (cell_type in c("CD4", "CD8", "B_and_plasma", "Myeloid")) { expect_true(is.list(regress_res[[cell_type]])) @@ -69,7 +69,7 @@ test_that("regressPC handles incorrect parameters", { cell_types = c("CD4", "CD8", "B_and_plasma", "Myeloid"), pc_subset = 1:5 )) - + # Test for invalid query data column expect_error(regressPC( reference_data = reference_data, @@ -79,7 +79,7 @@ test_that("regressPC handles incorrect parameters", { cell_types = c("CD4", "CD8", "B_and_plasma", "Myeloid"), pc_subset = 1:5 )) - + # Test for invalid cell types expect_error(regressPC( reference_data = reference_data, @@ -95,20 +95,20 @@ test_that("regressPC works with all cell types and default PC subset", { reference_data = reference_data, ref_cell_type_col = "expert_annotation" ) - + # Check if the output is a list expect_true(is.list(regress_res)) - + # Check if the list contains the correct elements expect_true("regression_summaries" %in% names(regress_res)) expect_true("r_squared" %in% names(regress_res)) expect_true("var_contributions" %in% names(regress_res)) expect_true("total_variance_explained" %in% names(regress_res)) - + # Check if R-squared values are numeric and of correct length expect_true(is.numeric(regress_res$r_squared)) expect_equal(length(regress_res$r_squared), 10) - + # Check if total variance explained is a numeric value expect_true(is.numeric(regress_res$total_variance_explained)) }) diff --git a/vignettes/VisualizationTools.Rmd b/vignettes/VisualizationTools.Rmd index 3f22f87..81ba8a6 100644 --- a/vignettes/VisualizationTools.Rmd +++ b/vignettes/VisualizationTools.Rmd @@ -44,7 +44,7 @@ dpi = 72, fig.retina = 2, fig.align = "center", out.width = "100%", -pngquant = "--speed=1 --quality=1-5" +pngquant = "--speed=1 --quality=5-10" ) ``` @@ -207,6 +207,43 @@ ggplot(distances_df, aes(x = Distance, fill = Anomaly)) + - **High Mahalanobis Distances**: Cells with high distances may represent outliers or novel, previously uncharacterized cell types. Further investigation may be needed to understand these anomalies. +## Project Data onto Sliced Inverse Regression (SIR) Space of Reference Data + +The `calculateSIRSpace()` function leverages Sliced Inverse Regression (SIR) to analyze high-dimensional cell type data by focusing on the directions that best separate different cell types and their subtypes. The function projects data onto a lower-dimensional SIR space, aiming to highlight the most informative features for distinguishing between cell types. + +The process starts by handling each cell type as a distinct slice of the data, allowing for a more detailed analysis. When `multiple_cond_means` is set to `TRUE`, the function refines this by computing conditional means for each cell type, using clustering to capture finer distinctions within each type. It then performs Singular Value Decomposition (SVD) on these refined means to uncover the principal directions that most effectively differentiate the cell types. + +Ultimately, the function provides projections that reveal how cell types and their subtypes relate to each other in a reduced space, facilitating a deeper understanding of cellular variability and distinctions in complex datasets. + +```{r, fig.height=5, fig.width=10, fig.show='hide'} +# Calculate the SIR space projections +sir_output <- calculateSIRSpace( + reference_data = reference_data, + query_data = query_data, + query_cell_type_col = "SingleR_annotation", + ref_cell_type_col = "expert_annotation", + multiple_cond_means = TRUE +) + +# Plot the SIR space projections +plot(sir_output, + plot_type = "boxplot", + sir_subset = 1:3) +``` +![](https://raw.githubusercontent.com/ccb-hms/scDiagnostics/main/inst/extdata/compressed/VisualizationTools/calculateSIRSpace1.png) + +In this example, we compute the SIR space projections using the `calculateSIRSpace()` function, and then we visualize the projections with a boxplot to analyze the distribution of the data in the SIR space. This example helps you understand how different cell types and their subtypes are projected into a common space, facilitating their comparison and interpretation. + +The plot below visualizes the contribution of different markers to the variance along the SIR dimensions, helping identify which markers are most influential in differentiating e.g. between B cells (including plasma cells) from other cell types. + +```{r, fig.height=5, fig.width=10, fig.show='hide'} +# Plot the top contributing markers +plot(sir_output, + plot_type = "varplot", + sir_subset = 1:3, + n_top_vars = 10) +``` +![](https://raw.githubusercontent.com/ccb-hms/scDiagnostics/main/inst/extdata/compressed/VisualizationTools/calculateSIRSpace2.png) # Visualization of Marker Expressions