diff --git a/R/application_visualization.R b/R/application_visualization.R index 20fea81..cdc9e01 100644 --- a/R/application_visualization.R +++ b/R/application_visualization.R @@ -552,19 +552,28 @@ make_heatmap_bidir_lt_ggplot = function(matrix, y_name, x_name, y_axis = TRUE, x #' #' @description \code{make_mushroom_plot} Make a plot in which each glyph consists of two semicircles corresponding to ligand- and receptor- information. The size of the semicircle is the percentage of cells that express the protein, while the saturation corresponds to the scaled average expression value. #' -#' @param prioritization_table A prioritization table as generated by \code{\link{generate_prioritization_tables}} +#' @param prioritization_table A prioritization table as generated by \code{\link{generate_prioritization_tables}}. #' @param top_n An integer indicating how many ligand-receptor pairs to show #' @param show_rankings A logical indicating whether to show the ranking of the ligand-receptor pairs (default: FALSE) #' @param show_all_datapoints A logical indicating whether to show all ligand-receptor pairs (default: FALSE, if true they will be grayed out) -#' @param true_color_range A logical indicating whether to use the true color range for the ligand-receptor pairs (default: FALSE; range 0-1 is used) +#' @param true_color_range A logical indicating whether to use the default color range as determined by ggplot (TRUE, default) or set the limits to a range of 0-1 (FALSE) +#' @param use_absolute_rank A logical indicating to whether use the absolute prioritization rank to filter the top_n ligand-receptor pairs (default: FALSE) #' @param size A string indicating which column to use for the size of the semicircles (default: "scaled_avg_exprs"; use column name without "_ligand" or "_receptor" suffix) -#' @param color A string indicating which column to use for the color of the semicircles (default: "scaled_lfc"; use column name without "_ligand" or "_receptor" suffix) +#' @param color A string indicating which column to use for the color of the semicircles (default: "scaled_p_val_adapted"; use column name without "_ligand" or "_receptor" suffix) #' @param ligand_fill_colors A vector of the low and high colors to use for the ligand semicircle fill gradient (default: c("#DEEBF7", "#08306B")) #' @param receptor_fill_colors A vector of the low and high colors to use for the receptor semicircle fill gradient (default: c("#FEE0D2", "#A50F15")) #' @param unranked_ligand_fill_colors A vector of the low and high colors to use for the unranked ligands when show_all_datapoints is TRUE (default: c(alpha("#FFFFFF", alpha=0.2), alpha("#252525", alpha=0.2))) -#' @param unranked_receptor_fill_colors A vector of the low and high colors to use for the unkraed receptors when show_all_datapoints is TRUE (default: c(alpha("#FFFFFF", alpha=0.2), alpha("#252525", alpha=0.2))) +#' @param unranked_receptor_fill_colors A vector of the low and high colors to use for the unranked receptors when show_all_datapoints is TRUE (default: c(alpha("#FFFFFF", alpha=0.2), alpha("#252525", alpha=0.2))) #' @param ... Additional arguments passed to \code{\link{ggplot2::theme}}. As there are often issues with the scales legend, it is recommended to change legend sizes and positions using this argument, i.e., \code{legend.key.height}, \code{legend.key.width}, \code{legend.title}, and \code{legend.text}. #' +#' @details +#' If the values range of the column used as the "size" parameter is not between 0 and 1.001, an error will be thrown. +#' +#' The sender cell types can be ordered by encoding the "sender" column as a factor. If the "sender" column is not a factor, the sender cell types will be ordered alphabetically. +#' +#' By default, the top_n ligand-receptor pairs are shown despite their absolute ranking. So, if a receiver cell type has LR pairs that are only ranked from 31-40 and the top_n is set to 20, the LR pairs will be shown. If use_absolute_rank is set to TRUE, only LR pairs with absolute ranking from 1-20 will be shown. +#' +#' #' @return A ggplot object #' #' @import ggplot2 @@ -589,14 +598,18 @@ make_heatmap_bidir_lt_ggplot = function(matrix, y_name, x_name, y_axis = TRUE, x #' #' # Change the size and color columns #' make_mushroom_plot(prior_table, size = "pct_expressed", color = "scaled_avg_exprs") -#' } #' #' +#' # For a prioritization table with multiple receiver cell types +#' make_mushroom_plot(prior_table_combined %>% filter(receiver == celltype_oi)) +#'} +#' #' @export #' make_mushroom_plot <- function(prioritization_table, top_n = 30, show_rankings = FALSE, - show_all_datapoints = FALSE, true_color_range = FALSE, - size = "scaled_avg_exprs", color = "scaled_lfc", + show_all_datapoints = FALSE, true_color_range = TRUE, + use_absolute_rank = FALSE, + size = "scaled_avg_exprs", color = "scaled_p_val_adapted", ligand_fill_colors = c("#DEEBF7", "#08306B"), receptor_fill_colors = c("#FEE0D2", "#A50F15"), unranked_ligand_fill_colors = c(alpha("#FFFFFF", alpha=0.2), alpha("#252525", alpha=0.2)), @@ -617,6 +630,8 @@ make_mushroom_plot <- function(prioritization_table, top_n = 30, show_rankings = stop("show_all_datapoints should be a TRUE or FALSE") if(!is.logical(true_color_range) | length(true_color_range) != 1) stop("true_color_range should be a TRUE or FALSE") + if(!is.logical(use_absolute_rank) | length(use_absolute_rank) != 1) + stop("use_absolute_rank should be a TRUE or FALSE") if(!is.numeric(top_n) | length(top_n) != 1) stop("top_n should be a numeric vector of length 1") if(length(ligand_fill_colors) != 2) @@ -635,27 +650,53 @@ make_mushroom_plot <- function(prioritization_table, top_n = 30, show_rankings = requireNamespace("shadowtext") requireNamespace("cowplot") - # Filter to top_n, create a new column of ligand-receptor interactions - filtered_table <- prioritization_table %>% dplyr::mutate(prioritization_rank = rank(desc(prioritization_score))) %>% - dplyr::mutate(lr_interaction = paste(ligand, receptor, sep = " - ")) - order_interactions <- unique(filtered_table %>% filter(prioritization_rank <= top_n) %>% pull(lr_interaction)) + if (!"prioritization_rank" %in% colnames(prioritization_table)){ + prioritization_table <- prioritization_table %>% dplyr::mutate(prioritization_rank = rank(desc(prioritization_score))) + } + # Add 'relative rank' column which is basically 1:n + prioritization_table <- prioritization_table %>% dplyr::mutate(relative_rank = rank(desc(prioritization_score))) + + # If use_absolute_rank, use 'prioritization_rank' column to filter top_n + rank_filter_col <- ifelse(use_absolute_rank, "prioritization_rank", "relative_rank") + + # Create a new column of ligand-receptor interactions, and filter table to + # only include LR interactions that appear in the top_n + filtered_table <- prioritization_table %>% dplyr::mutate(lr_interaction = paste(ligand, receptor, sep = " - ")) + order_interactions <- unique(filtered_table %>% filter(.data[[rank_filter_col]] <= top_n) %>% pull(lr_interaction)) filtered_table <- filtered_table %>% filter(lr_interaction %in% order_interactions) %>% mutate(lr_interaction = factor(lr_interaction, levels = rev(order_interactions))) - celltypes_vec <- 1:length(unique(filtered_table$sender)) %>% setNames(sort(unique(filtered_table$sender))) + # Check if filtered_table is empty + if (nrow(filtered_table) == 0){ + stop("No ligand-receptor interactions found in the top_n. Please try use_absolute_rank = FALSE or increase top_n.") + } + + # Keep order of senders, if present (if not, sort alphabetically) + if (!is.factor(filtered_table$sender)){ + filtered_table$sender <- as.factor(filtered_table$sender) + } else { + # Drop levels that are not present in the filtered table + filtered_table$sender <- droplevels(filtered_table$sender) + } + lr_interaction_vec <- 1:length(order_interactions) %>% setNames(order_interactions) # Make each ligand and receptor into separate rows (to draw 1 semicircle per row) - filtered_table <- filtered_table %>% select(c("lr_interaction", all_of(cols_to_use), "prioritization_rank")) %>% + filtered_table <- filtered_table %>% select(c("lr_interaction", all_of(cols_to_use), "prioritization_rank", "relative_rank")) %>% pivot_longer(c(ligand, receptor), names_to = "type", values_to = "protein") %>% mutate(size = ifelse(type == "ligand", get(paste0(size, "_", size_ext[1])), get(paste0(size, "_", size_ext[2]))), color = ifelse(type == "ligand", get(paste0(color, "_", color_ext[1])), get(paste0(color, "_", color_ext[2])))) %>% select(-contains(c("_ligand", "_receptor", "_sender", "_receiver"))) %>% mutate(start = rep(c(-pi, 0), nrow(filtered_table))) %>% - mutate(x = celltypes_vec[sender], y = lr_interaction_vec[lr_interaction]) + mutate(x = as.numeric(sender), y = lr_interaction_vec[lr_interaction]) + + # Warning if size column is not scaled between 0 and 1.001 + if (any(filtered_table$size < 0) | any(filtered_table$size > 1.001)){ + stop("Size column is not scaled between 0 and 1. Please use this column as the color instead.") + } # Rename size and color columns to be more human-readable - keywords_adj <- c("LFC", "p-val", "product", "mean", "adjusted", "expression") %>% setNames(c("lfc", "pval", "prod", "avg", "adj", "exprs")) + keywords_adj <- c("LFC", "pval", "", "product", "mean", "adjusted", "expression") %>% setNames(c("lfc", "p", "val", "prod", "avg", "adj", "exprs")) size_title <- sapply(stringr::str_split(size, "_")[[1]], function(k) ifelse(is.na(keywords_adj[k]), k, keywords_adj[k])) %>% paste0(., collapse = " ") %>% stringr::str_replace("^\\w{1}", toupper) color_title <- sapply(stringr::str_split(color, "_")[[1]], function(k) ifelse(is.na(keywords_adj[k]), k, keywords_adj[k])) %>% @@ -666,7 +707,7 @@ make_mushroom_plot <- function(prioritization_table, top_n = 30, show_rankings = scale <- 0.5 - ncelltypes <- length(celltypes_vec) + ncelltypes <- length(unique(filtered_table$sender)) n_interactions <- length(lr_interaction_vec) legend2_df <- data.frame(values = c(0.25, 0.5, 0.75, 1), x=(ncelltypes+2.5):(ncelltypes+5.5), y=rep(floor(n_interactions/3), 4), start=-pi) axis_rect <- data.frame(xmin=0, xmax=ncelltypes+1, ymin=0, ymax=n_interactions+1) @@ -699,7 +740,7 @@ make_mushroom_plot <- function(prioritization_table, top_n = 30, show_rankings = p1 <- ggplot() + # Draw ligand semicircle - geom_arc_bar(data = filtered_table %>% filter(type=="ligand", prioritization_rank <= top_n), + geom_arc_bar(data = filtered_table %>% filter(type=="ligand", .data[[rank_filter_col]] <= top_n), aes(x0 = x, y0 = y, r0 = 0, r = sqrt(size)*scale, start = start, end = start + pi, fill=color), color = "white") + @@ -707,10 +748,10 @@ make_mushroom_plot <- function(prioritization_table, top_n = 30, show_rankings = limits=color_lims, oob=scales::squish, n.breaks = 3, guide = guide_colorbar(order = 1), - name=paste0(color_title, " (", color_ext[1], ")") %>% str_wrap(width=15)) + + name=paste0(color_title, " (", color_ext[1], ")") %>% stringr::str_wrap(width=15)) + # Create new fill scale for receptor semicircles new_scale_fill() + - geom_arc_bar(data = filtered_table %>% filter(type=="receptor", prioritization_rank <= top_n), + geom_arc_bar(data = filtered_table %>% filter(type=="receptor", .data[[rank_filter_col]] <= top_n), aes(x0 = x, y0 = y, r0 = 0, r = sqrt(size)*scale, start = start, end = start + pi, fill=color), color = "white") + @@ -719,7 +760,7 @@ make_mushroom_plot <- function(prioritization_table, top_n = 30, show_rankings = geom_rect(data = legend2_df, aes(xmin=x-0.5, xmax=x+0.5, ymin=y-0.5, ymax=y+0.5), color="gray90", fill=NA) + geom_text(data = legend2_df, aes(label=values, x=x, y=y-0.6), vjust=1, size = scale_legend_text_size) + geom_text(data = data.frame(x = (ncelltypes+4), y = floor(n_interactions/3)+1, - label = size_title %>% str_wrap(width=15)), + label = size_title %>% stringr::str_wrap(width=15)), aes(x=x, y=y, label=label), size = scale_legend_title_size, vjust=0, lineheight = .75) + # Panel grid geom_line(data = panel_grid_y, aes(x=x, y=y, group=group), color = "gray90") + @@ -728,10 +769,10 @@ make_mushroom_plot <- function(prioritization_table, top_n = 30, show_rankings = geom_rect(data = axis_rect, aes(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax), color = "black", fill = "transparent") + # Other plot information scale_fill_gradient(low = receptor_fill_colors[1], high=receptor_fill_colors[2] , limits=color_lims, oob=scales::squish, n.breaks = 3, - name=paste0(color_title, " (", color_ext[2], ")") %>% str_wrap(width=15), + name=paste0(color_title, " (", color_ext[2], ")") %>% stringr::str_wrap(width=15), guide = guide_colorbar(order = 2)) + scale_y_continuous(breaks=n_interactions:1, labels=names(lr_interaction_vec), expand = expansion(add=c(0,0))) + - scale_x_continuous(breaks=1:ncelltypes, labels=names(celltypes_vec), position="top", expand = expansion(add=c(0,0))) + + scale_x_continuous(breaks=1:ncelltypes, labels=levels(filtered_table$sender), position="top", expand = expansion(add=c(0,0))) + xlab("Sender cell types") + ylab("Ligand-receptor interaction") + coord_fixed() + do.call(theme, theme_args) @@ -743,14 +784,14 @@ make_mushroom_plot <- function(prioritization_table, top_n = 30, show_rankings = unranked_ligand_lims <- c(0,1); unranked_receptor_lims <- c(0,1) if (true_color_range){ # Follow limits of the top_n lr pairs - unranked_ligand_lims <- filtered_table %>% filter(type=="ligand", prioritization_rank <= top_n) %>% + unranked_ligand_lims <- filtered_table %>% filter(type=="ligand", .data[[rank_filter_col]] <= top_n) %>% select(color) %>% range - unranked_receptor_lims <- filtered_table %>% filter(type=="receptor", prioritization_rank <= top_n) %>% + unranked_receptor_lims <- filtered_table %>% filter(type=="receptor", .data[[rank_filter_col]] <= top_n) %>% select(color) %>% range } p1 <- p1 + new_scale_fill() + - geom_arc_bar(data = filtered_table %>% filter(type=="ligand", prioritization_rank > top_n), + geom_arc_bar(data = filtered_table %>% filter(type=="ligand", .data[[rank_filter_col]] > top_n), aes(x0 = x, y0 = y, r0 = 0, r = sqrt(size)*scale, start = start, end = start + pi, fill=color), color = "white") + @@ -758,7 +799,7 @@ make_mushroom_plot <- function(prioritization_table, top_n = 30, show_rankings = limits=unranked_ligand_lims, oob = scales::oob_squish, guide = "none") + new_scale_fill() + - geom_arc_bar(data = filtered_table %>% filter(type=="receptor", prioritization_rank > top_n), + geom_arc_bar(data = filtered_table %>% filter(type=="receptor", .data[[rank_filter_col]] > top_n), aes(x0 = x, y0 = y, r0 = 0, r = sqrt(size)*scale, start = start, end = start + pi, fill=color), color = "white") + @@ -769,7 +810,7 @@ make_mushroom_plot <- function(prioritization_table, top_n = 30, show_rankings = # Add ranking numbers if requested if (show_rankings){ - p1 <- p1 + geom_shadowtext(data = filtered_table %>% filter(prioritization_rank <= top_n), + p1 <- p1 + geom_shadowtext(data = filtered_table %>% filter(.data[[rank_filter_col]] <= top_n), aes(x=x, y=y, label=prioritization_rank)) } diff --git a/R/prioritization.R b/R/prioritization.R index 3beca2c..79b2933 100644 --- a/R/prioritization.R +++ b/R/prioritization.R @@ -401,11 +401,11 @@ generate_info_tables <- function(seuratObj, #' Additionally, the following columns are added: #' \itemize{ #' \item \code{lfc_pval_*}: product of -log10(pval) and the LFC of the ligand/receptor -#' \item \code{p_val_*_adapted}: p-value adapted to the sign of the LFC to only consider interactions where the ligand/receptor is upregulated in the sender/receiver +#' \item \code{p_val_adapted_*}: p-value adapted to the sign of the LFC to only consider interactions where the ligand/receptor is upregulated in the sender/receiver #' \item \code{activity_zscore}: z-score of the ligand activity #' \item \code{prioritization_score}: The prioritization score for each interaction, calculated as a weighted sum of the prioritization criteria. #' } -#' Moreover, \code{scaled_*} columns are scaled using the corresponding column's ranking or the \code{scale_quantile_adapted} function. The columns used for prioritization are scaled_p_val_ligand_adapted, scaled_p_val_receptor_adapted, scaled_activity, scaled_avg_exprs_ligand, scaled_avg_exprs_receptor, scaled_p_val_ligand_adapted_group, scaled_p_val_receptor_adapted_group +#' Moreover, \code{scaled_*} columns are scaled using the corresponding column's ranking or the \code{scale_quantile_adapted} function. The columns used for prioritization are scaled_p_val_adapted_ligand, scaled_p_val_adapted_receptor, scaled_activity, scaled_avg_exprs_ligand, scaled_avg_exprs_receptor, scaled_p_val_adapted_ligand_group, scaled_p_val_adapted_receptor_group #' #' @import dplyr #' @@ -511,21 +511,21 @@ generate_prioritization_tables = function(sender_receiver_info, sender_receiver_ sender_ligand_prioritization = sender_receiver_de %>% dplyr::ungroup() %>% dplyr::select(sender, ligand, lfc_ligand, p_val_ligand) %>% dplyr::distinct() %>% dplyr::mutate(lfc_pval_ligand = -log10(p_val_ligand)*lfc_ligand, - p_val_ligand_adapted = -log10(p_val_ligand)*sign(lfc_ligand)) + p_val_adapted_ligand = -log10(p_val_ligand)*sign(lfc_ligand)) sender_ligand_prioritization = sender_ligand_prioritization %>% dplyr::mutate(scaled_lfc_ligand = rank(lfc_ligand, ties.method = "average", na.last = FALSE)/max(rank(lfc_ligand, ties.method = "average", na.last = FALSE)), scaled_p_val_ligand = rank(desc(p_val_ligand), ties.method = "average", na.last = FALSE)/max(rank(desc(p_val_ligand), ties.method = "average", na.last = FALSE)), scaled_lfc_pval_ligand = rank(lfc_pval_ligand, ties.method = "average", na.last = FALSE)/max(rank(lfc_pval_ligand, ties.method = "average", na.last = FALSE)), - scaled_p_val_ligand_adapted = rank(p_val_ligand_adapted, ties.method = "average", na.last = FALSE)/max(rank(p_val_ligand_adapted, ties.method = "average", na.last = FALSE))) %>% + scaled_p_val_adapted_ligand = rank(p_val_adapted_ligand, ties.method = "average", na.last = FALSE)/max(rank(p_val_adapted_ligand, ties.method = "average", na.last = FALSE))) %>% dplyr::arrange(-lfc_pval_ligand) # Receptor DE prioritization receiver_receptor_prioritization = sender_receiver_de %>% dplyr::ungroup() %>% dplyr::select(receiver, receptor, lfc_receptor, p_val_receptor) %>% dplyr::distinct() %>% dplyr::mutate(lfc_pval_receptor = -log10(p_val_receptor)*lfc_receptor, - p_val_receptor_adapted = -log10(p_val_receptor)*sign(lfc_receptor) ) + p_val_adapted_receptor = -log10(p_val_receptor)*sign(lfc_receptor) ) receiver_receptor_prioritization = receiver_receptor_prioritization %>% dplyr::mutate(scaled_lfc_receptor = rank(lfc_receptor, ties.method = "average", na.last = FALSE)/max(rank(lfc_receptor, ties.method = "average", na.last = FALSE)), scaled_p_val_receptor = rank(desc(p_val_receptor), ties.method = "average", na.last = FALSE)/max(rank(desc(p_val_receptor), ties.method = "average", na.last = FALSE)), scaled_lfc_pval_receptor = rank(lfc_pval_receptor, ties.method = "average", na.last = FALSE)/max(rank(lfc_pval_receptor, ties.method = "average", na.last = FALSE)), - scaled_p_val_receptor_adapted = rank(p_val_receptor_adapted, ties.method = "average", na.last = FALSE)/max(rank(p_val_receptor_adapted, ties.method = "average", na.last = FALSE))) %>% dplyr::arrange(-lfc_pval_receptor) + scaled_p_val_adapted_receptor = rank(p_val_adapted_receptor, ties.method = "average", na.last = FALSE)/max(rank(p_val_adapted_receptor, ties.method = "average", na.last = FALSE))) %>% dplyr::arrange(-lfc_pval_receptor) # Ligand activity prioritization ligand_activity_prioritization = ligand_activities %>% @@ -547,21 +547,21 @@ generate_prioritization_tables = function(sender_receiver_info, sender_receiver_ # Condition specificity of ligand (upregulation) ligand_condition_prioritization = lr_condition_de %>% dplyr::ungroup() %>% dplyr::select(ligand, lfc_ligand, p_val_ligand) %>% dplyr::distinct() %>% dplyr::mutate(lfc_pval_ligand = -log10(p_val_ligand)*lfc_ligand, - p_val_ligand_adapted = -log10(p_val_ligand)*sign(lfc_ligand)) + p_val_adapted_ligand = -log10(p_val_ligand)*sign(lfc_ligand)) ligand_condition_prioritization = ligand_condition_prioritization %>% dplyr::mutate(scaled_lfc_ligand = rank(lfc_ligand, ties.method = "average", na.last = FALSE)/max(rank(lfc_ligand, ties.method = "average", na.last = FALSE)), scaled_p_val_ligand = rank(desc(p_val_ligand), ties.method = "average", na.last = FALSE)/max(rank(desc(p_val_ligand), ties.method = "average", na.last = FALSE)), scaled_lfc_pval_ligand = rank(lfc_pval_ligand, ties.method = "average", na.last = FALSE)/max(rank(lfc_pval_ligand, ties.method = "average", na.last = FALSE)), - scaled_p_val_ligand_adapted = rank(p_val_ligand_adapted, ties.method = "average", na.last = FALSE)/max(rank(p_val_ligand_adapted, ties.method = "average", na.last = FALSE))) %>% + scaled_p_val_adapted_ligand = rank(p_val_adapted_ligand, ties.method = "average", na.last = FALSE)/max(rank(p_val_adapted_ligand, ties.method = "average", na.last = FALSE))) %>% dplyr::arrange(-lfc_pval_ligand) %>% rename_with(.fn = function(column_name) paste0(column_name, "_group"), .cols = -ligand) # Condition specificity of receptor (upregulation) receptor_condition_prioritization = lr_condition_de %>% dplyr::ungroup() %>% dplyr::select(receptor, lfc_receptor, p_val_receptor) %>% dplyr::distinct() %>% dplyr::mutate(lfc_pval_receptor = -log10(p_val_receptor)*lfc_receptor, - p_val_receptor_adapted = -log10(p_val_receptor)*sign(lfc_receptor)) + p_val_adapted_receptor = -log10(p_val_receptor)*sign(lfc_receptor)) receptor_condition_prioritization = receptor_condition_prioritization %>% dplyr::mutate(scaled_lfc_receptor = rank(lfc_receptor, ties.method = "average", na.last = FALSE)/max(rank(lfc_receptor, ties.method = "average", na.last = FALSE)), scaled_p_val_receptor = rank(desc(p_val_receptor), ties.method = "average", na.last = FALSE)/max(rank(desc(p_val_receptor), ties.method = "average", na.last = FALSE)), scaled_lfc_pval_receptor = rank(lfc_pval_receptor, ties.method = "average", na.last = FALSE)/max(rank(lfc_pval_receptor, ties.method = "average", na.last = FALSE)), - scaled_p_val_receptor_adapted = rank(p_val_receptor_adapted, ties.method = "average", na.last = FALSE)/max(rank(p_val_receptor_adapted, ties.method = "average", na.last = FALSE))) %>% + scaled_p_val_adapted_receptor = rank(p_val_adapted_receptor, ties.method = "average", na.last = FALSE)/max(rank(p_val_adapted_receptor, ties.method = "average", na.last = FALSE))) %>% dplyr::arrange(-lfc_pval_receptor) %>% rename_with(.fn = function(column_name) paste0(column_name, "_group"), .cols = -receptor) } else { @@ -584,20 +584,23 @@ generate_prioritization_tables = function(sender_receiver_info, sender_receiver_ # have a weighted average the final score (no product!!) - sum_prioritization_weights = weights["de_ligand"] + weights["de_receptor"] + weights["activity_scaled"] + weights["exprs_ligand"] + weights["exprs_receptor"] + weights["ligand_condition_specificity"] + weights["receptor_condition_specificity"] + sum_prioritization_weights = 0.5*weights["de_ligand"] + 0.5*weights["de_receptor"] + weights["activity_scaled"] + 0.5*weights["exprs_ligand"] + 0.5*weights["exprs_receptor"] + weights["ligand_condition_specificity"] + weights["receptor_condition_specificity"] group_prioritization_tbl = group_prioritization_tbl %>% rowwise() %>% dplyr::mutate(prioritization_score = ( - (prioritizing_weights["de_ligand"] * ifelse("scaled_p_val_ligand_adapted" %in% names(group_prioritization_tbl), scaled_p_val_ligand_adapted, 0)) + - (prioritizing_weights["de_receptor"] * ifelse("scaled_p_val_receptor_adapted" %in% names(group_prioritization_tbl), scaled_p_val_receptor_adapted, 0)) + + (0.5*prioritizing_weights["de_ligand"] * ifelse("scaled_p_val_adapted_ligand" %in% names(group_prioritization_tbl), scaled_p_val_adapted_ligand, 0)) + + (0.5*prioritizing_weights["de_receptor"] * ifelse("scaled_p_val_adapted_receptor" %in% names(group_prioritization_tbl), scaled_p_val_adapted_receptor, 0)) + (prioritizing_weights["activity_scaled"] * ifelse("scaled_activity" %in% names(group_prioritization_tbl), scaled_activity, 0)) + - (prioritizing_weights["exprs_ligand"] * ifelse("scaled_avg_exprs_ligand" %in% names(group_prioritization_tbl), scaled_avg_exprs_ligand, 0)) + - (prioritizing_weights["exprs_receptor"] * ifelse("scaled_avg_exprs_receptor" %in% names(group_prioritization_tbl), scaled_avg_exprs_receptor, 0)) + - (prioritizing_weights["ligand_condition_specificity"] * ifelse("scaled_p_val_ligand_adapted_group" %in% names(group_prioritization_tbl), scaled_p_val_ligand_adapted_group, 0)) + - (prioritizing_weights["receptor_condition_specificity"] * ifelse("scaled_p_val_receptor_adapted_group" %in% names(group_prioritization_tbl), scaled_p_val_receptor_adapted_group, 0)) + (0.5*prioritizing_weights["exprs_ligand"] * ifelse("scaled_avg_exprs_ligand" %in% names(group_prioritization_tbl), scaled_avg_exprs_ligand, 0)) + + (0.5*prioritizing_weights["exprs_receptor"] * ifelse("scaled_avg_exprs_receptor" %in% names(group_prioritization_tbl), scaled_avg_exprs_receptor, 0)) + + (prioritizing_weights["ligand_condition_specificity"] * ifelse("scaled_p_val_adapted_ligand_group" %in% names(group_prioritization_tbl), scaled_p_val_adapted_ligand_group, 0)) + + (prioritizing_weights["receptor_condition_specificity"] * ifelse("scaled_p_val_adapted_receptor_group" %in% names(group_prioritization_tbl), scaled_p_val_adapted_receptor_group, 0)) )* (1/sum_prioritization_weights)) %>% dplyr::arrange(-prioritization_score) %>% ungroup() + # Add rank + group_prioritization_tbl = group_prioritization_tbl %>% dplyr::mutate(prioritization_rank = rank(desc(prioritization_score))) + return (group_prioritization_tbl) } diff --git a/man/generate_prioritization_tables.Rd b/man/generate_prioritization_tables.Rd index 33aa3ee..719a002 100644 --- a/man/generate_prioritization_tables.Rd +++ b/man/generate_prioritization_tables.Rd @@ -32,11 +32,11 @@ The resulting dataframe contains columns from the input dataframes, but columns Additionally, the following columns are added: \itemize{ \item \code{lfc_pval_*}: product of -log10(pval) and the LFC of the ligand/receptor -\item \code{p_val_*_adapted}: p-value adapted to the sign of the LFC to only consider interactions where the ligand/receptor is upregulated in the sender/receiver +\item \code{p_val_adapted_*}: p-value adapted to the sign of the LFC to only consider interactions where the ligand/receptor is upregulated in the sender/receiver \item \code{activity_zscore}: z-score of the ligand activity \item \code{prioritization_score}: The prioritization score for each interaction, calculated as a weighted sum of the prioritization criteria. } -Moreover, \code{scaled_*} columns are scaled using the corresponding column's ranking or the \code{scale_quantile_adapted} function. The columns used for prioritization are scaled_p_val_ligand_adapted, scaled_p_val_receptor_adapted, scaled_activity, scaled_avg_exprs_ligand, scaled_avg_exprs_receptor, scaled_p_val_ligand_adapted_group, scaled_p_val_receptor_adapted_group +Moreover, \code{scaled_*} columns are scaled using the corresponding column's ranking or the \code{scale_quantile_adapted} function. The columns used for prioritization are scaled_p_val_adapted_ligand, scaled_p_val_adapted_receptor, scaled_activity, scaled_avg_exprs_ligand, scaled_avg_exprs_receptor, scaled_p_val_adapted_ligand_group, scaled_p_val_adapted_receptor_group } \description{ User can choose the importance attached to each of the following prioritization criteria: differential expression of ligand and receptor, cell-type specificity of expression of ligand and receptor, NicheNet ligand activity diff --git a/man/make_mushroom_plot.Rd b/man/make_mushroom_plot.Rd index 2d02d29..e3d5a82 100644 --- a/man/make_mushroom_plot.Rd +++ b/man/make_mushroom_plot.Rd @@ -9,20 +9,21 @@ make_mushroom_plot( top_n = 30, show_rankings = FALSE, show_all_datapoints = FALSE, - true_color_range = FALSE, + true_color_range = TRUE, + use_absolute_rank = FALSE, size = "scaled_avg_exprs", - color = "scaled_lfc", + color = "scaled_p_val_adapted", ligand_fill_colors = c("#DEEBF7", "#08306B"), receptor_fill_colors = c("#FEE0D2", "#A50F15"), - unranked_ligand_fill_colors = c(alpha("#FFFFFF", alpha = 0.2), alpha("#252525", alpha + unranked_ligand_fill_colors = c(alpha("#FFFFFF", alpha = 0.2), alpha("#252525", alpha = + 0.2)), + unranked_receptor_fill_colors = c(alpha("#FFFFFF", alpha = 0.2), alpha("#252525", alpha = 0.2)), - unranked_receptor_fill_colors = c(alpha("#FFFFFF", alpha = 0.2), alpha("#252525", - alpha = 0.2)), ... ) } \arguments{ -\item{prioritization_table}{A prioritization table as generated by \code{\link{generate_prioritization_tables}}} +\item{prioritization_table}{A prioritization table as generated by \code{\link{generate_prioritization_tables}}.} \item{top_n}{An integer indicating how many ligand-receptor pairs to show} @@ -30,11 +31,13 @@ make_mushroom_plot( \item{show_all_datapoints}{A logical indicating whether to show all ligand-receptor pairs (default: FALSE, if true they will be grayed out)} -\item{true_color_range}{A logical indicating whether to use the true color range for the ligand-receptor pairs (default: FALSE; range 0-1 is used)} +\item{true_color_range}{A logical indicating whether to use the default color range as determined by ggplot (TRUE, default) or set the limits to a range of 0-1 (FALSE)} + +\item{use_absolute_rank}{A logical indicating to whether use the absolute prioritization rank to filter the top_n ligand-receptor pairs (default: FALSE)} \item{size}{A string indicating which column to use for the size of the semicircles (default: "scaled_avg_exprs"; use column name without "_ligand" or "_receptor" suffix)} -\item{color}{A string indicating which column to use for the color of the semicircles (default: "scaled_lfc"; use column name without "_ligand" or "_receptor" suffix)} +\item{color}{A string indicating which column to use for the color of the semicircles (default: "scaled_p_val_adapted"; use column name without "_ligand" or "_receptor" suffix)} \item{ligand_fill_colors}{A vector of the low and high colors to use for the ligand semicircle fill gradient (default: c("#DEEBF7", "#08306B"))} @@ -42,7 +45,7 @@ make_mushroom_plot( \item{unranked_ligand_fill_colors}{A vector of the low and high colors to use for the unranked ligands when show_all_datapoints is TRUE (default: c(alpha("#FFFFFF", alpha=0.2), alpha("#252525", alpha=0.2)))} -\item{unranked_receptor_fill_colors}{A vector of the low and high colors to use for the unkraed receptors when show_all_datapoints is TRUE (default: c(alpha("#FFFFFF", alpha=0.2), alpha("#252525", alpha=0.2)))} +\item{unranked_receptor_fill_colors}{A vector of the low and high colors to use for the unranked receptors when show_all_datapoints is TRUE (default: c(alpha("#FFFFFF", alpha=0.2), alpha("#252525", alpha=0.2)))} \item{...}{Additional arguments passed to \code{\link{ggplot2::theme}}. As there are often issues with the scales legend, it is recommended to change legend sizes and positions using this argument, i.e., \code{legend.key.height}, \code{legend.key.width}, \code{legend.title}, and \code{legend.text}.} } @@ -52,6 +55,13 @@ A ggplot object \description{ \code{make_mushroom_plot} Make a plot in which each glyph consists of two semicircles corresponding to ligand- and receptor- information. The size of the semicircle is the percentage of cells that express the protein, while the saturation corresponds to the scaled average expression value. } +\details{ +If the values range of the column used as the "size" parameter is not between 0 and 1.001, an error will be thrown. + +The sender cell types can be ordered by encoding the "sender" column as a factor. If the "sender" column is not a factor, the sender cell types will be ordered alphabetically. + +By default, the top_n ligand-receptor pairs are shown despite their absolute ranking. So, if a receiver cell type has LR pairs that are only ranked from 31-40 and the top_n is set to 20, the LR pairs will be shown. If use_absolute_rank is set to TRUE, only LR pairs with absolute ranking from 1-20 will be shown. +} \examples{ \dontrun{ # Create a prioritization table @@ -68,7 +78,10 @@ make_mushroom_plot(prior_table, show_all_datapoints = TRUE, true_color_range = T # Change the size and color columns make_mushroom_plot(prior_table, size = "pct_expressed", color = "scaled_avg_exprs") -} +# For a prioritization table with multiple receiver cell types +make_mushroom_plot(prior_table_combined \%>\% filter(receiver == celltype_oi)) +} + } diff --git a/tests/testthat/test-prioritization.R b/tests/testthat/test-prioritization.R index 9838b9d..afd0c0f 100644 --- a/tests/testthat/test-prioritization.R +++ b/tests/testthat/test-prioritization.R @@ -206,41 +206,41 @@ test_that("Prioritization scheme works", { expect_identical(info_tables$sender_receiver_de, processed_DE_table) expect_identical(info_tables$lr_condition_de, processed_condition_markers) - # Check processed tables - expect_type(processed_DE_table,"list") - expect_type(processed_expr_table,"list") - expect_type(processed_condition_markers,"list") - - # Check that processed_DE_table has been processed correctly - expect_equal(length(unique(processed_DE_table$sender)), length(sender_celltypes)) - expect_equal(length(unique(processed_DE_table$receiver)), length(receiver)) - - expect_equal(processed_DE_table %>% filter(ligand == "Il1rn", sender == "Mono") %>% select(p_val_ligand, lfc_ligand, pct_expressed_sender, p_adj_ligand, sender, ligand) %>% distinct(across(everything())), - DE_table %>% filter(gene == "Il1rn", cluster_id == "Mono") %>% select(-pct.2), - check.attributes = FALSE) - expect_equal(processed_DE_table %>% filter(ligand == "Il1rn", sender == "Mono", receptor == "Il1r2") %>% select(p_val_receptor, lfc_receptor, pct_expressed_receiver, p_adj_receptor, receiver, receptor) , - DE_table %>% filter(gene == "Il1r2", cluster_id == "CD8 T") %>% select(-pct.2), - check.attributes = FALSE) - temp_row <- processed_DE_table %>% filter(ligand == "Il1rn", sender == "Mono", receiver == "CD8 T", receptor == "Il1r2") - expect_equal((temp_row$lfc_ligand + temp_row$lfc_receptor)/2, temp_row$ligand_receptor_lfc_avg) - - - # Check that processed_expr_table has been processed correctly - expect_equal(length(unique(processed_expr_table$sender)), length(celltypes)) - expect_equal(length(unique(processed_expr_table$receiver)), length(celltypes)) - - temp_row <- processed_expr_table %>% filter(ligand == "Il1rn", sender == "Mono", receiver == "CD8 T", receptor == "Il1r2") - expect_equal(temp_row$avg_ligand * temp_row$avg_receptor, temp_row$ligand_receptor_prod) - - # Check that processed_condition_markers has been processed correctly - expect_equal(processed_condition_markers %>% filter(ligand == "Il1rn") %>% select(ligand, p_val_ligand, lfc_ligand, p_adj_ligand) %>% distinct(across(everything())), - condition_markers %>% filter(gene == "Il1rn") %>% select(-pct.1, -pct.2), - check.attributes = FALSE) - expect_equal(processed_condition_markers %>% filter(receptor == "Il1r2") %>% select(receptor, p_val_receptor, lfc_receptor, p_adj_receptor) %>% distinct(across(everything())), - condition_markers %>% filter(gene == "Il1r2") %>% select(-pct.1, -pct.2), - check.attributes = FALSE) - temp_row <- processed_condition_markers %>% filter(ligand == "Il1rn", receptor == "Il1r2") - expect_equal((temp_row$lfc_ligand + temp_row$lfc_receptor)/2, temp_row$ligand_receptor_lfc_avg) + # Check processed tables + expect_type(processed_DE_table,"list") + expect_type(processed_expr_table,"list") + expect_type(processed_condition_markers,"list") + + # Check that processed_DE_table has been processed correctly + expect_equal(length(unique(processed_DE_table$sender)), length(sender_celltypes)) + expect_equal(length(unique(processed_DE_table$receiver)), length(receiver)) + + expect_equal(processed_DE_table %>% filter(ligand == "Il1rn", sender == "Mono") %>% select(p_val_ligand, lfc_ligand, pct_expressed_sender, p_adj_ligand, sender, ligand) %>% distinct(across(everything())), + DE_table %>% filter(gene == "Il1rn", cluster_id == "Mono") %>% select(-pct.2), + check.attributes = FALSE) + expect_equal(processed_DE_table %>% filter(ligand == "Il1rn", sender == "Mono", receptor == "Il1r2") %>% select(p_val_receptor, lfc_receptor, pct_expressed_receiver, p_adj_receptor, receiver, receptor) , + DE_table %>% filter(gene == "Il1r2", cluster_id == "CD8 T") %>% select(-pct.2), + check.attributes = FALSE) + temp_row <- processed_DE_table %>% filter(ligand == "Il1rn", sender == "Mono", receiver == "CD8 T", receptor == "Il1r2") + expect_equal((temp_row$lfc_ligand + temp_row$lfc_receptor)/2, temp_row$ligand_receptor_lfc_avg) + + + # Check that processed_expr_table has been processed correctly + expect_equal(length(unique(processed_expr_table$sender)), length(celltypes)) + expect_equal(length(unique(processed_expr_table$receiver)), length(celltypes)) + + temp_row <- processed_expr_table %>% filter(ligand == "Il1rn", sender == "Mono", receiver == "CD8 T", receptor == "Il1r2") + expect_equal(temp_row$avg_ligand * temp_row$avg_receptor, temp_row$ligand_receptor_prod) + + # Check that processed_condition_markers has been processed correctly + expect_equal(processed_condition_markers %>% filter(ligand == "Il1rn") %>% select(ligand, p_val_ligand, lfc_ligand, p_adj_ligand) %>% distinct(across(everything())), + condition_markers %>% filter(gene == "Il1rn") %>% select(-pct.1, -pct.2), + check.attributes = FALSE) + expect_equal(processed_condition_markers %>% filter(receptor == "Il1r2") %>% select(receptor, p_val_receptor, lfc_receptor, p_adj_receptor) %>% distinct(across(everything())), + condition_markers %>% filter(gene == "Il1r2") %>% select(-pct.1, -pct.2), + check.attributes = FALSE) + temp_row <- processed_condition_markers %>% filter(ligand == "Il1rn", receptor == "Il1r2") + expect_equal((temp_row$lfc_ligand + temp_row$lfc_receptor)/2, temp_row$ligand_receptor_lfc_avg) # Check errors and warnings in case of improper usage expect_error(process_table_to_ic(condition_markers, table_type = "group_DE", lr_network = lr_network, receivers_oi = receiver)) @@ -251,109 +251,175 @@ test_that("Prioritization scheme works", { expect_warning(process_table_to_ic(expression_info, lr_network = lr_network, receivers_oi = receiver)) expect_warning(process_table_to_ic(expression_info, lr_network = lr_network, senders_oi = sender_celltypes)) - # Default weights - prioritizing_weights = c("de_ligand" = 1, - "de_receptor" = 1, - "activity_scaled" = 2, - "exprs_ligand" = 1, - "exprs_receptor" = 1, - "ligand_condition_specificity" = 0.5, - "receptor_condition_specificity" = 0.5) - + # Default weights prior_table <- generate_prioritization_tables(processed_expr_table, - processed_DE_table, - ligand_activities, - processed_condition_markers, - prioritizing_weights = prioritizing_weights) - - expect_type(prior_table,"list") - - # Check that columns contain columns from processed_DE_table, processed_expr_table, ligand_activities, and processed_condition_markers - expect_equal(prior_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% select(colnames(processed_DE_table)), - processed_DE_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% mutate(sender = as.character(sender), receiver = as.character(receiver)), - check.attributes = FALSE) - expect_equal(prior_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% select(colnames(processed_expr_table)), - processed_expr_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69", receiver == "CD8 T"), - check.attributes = FALSE) - expect_equal(prior_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% pull(activity), - ligand_activities %>% filter(test_ligand == "Lgals1") %>% pull(aupr_corrected)) - temp_cols <- c("lfc_ligand", "lfc_receptor", "p_val_ligand", "p_val_receptor") - expect_equal(prior_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% select(all_of(paste0(temp_cols, "_group"))), - processed_condition_markers %>% filter(ligand == "Lgals1", receptor == "Cd69") %>% select(all_of(temp_cols)), - check.attributes = FALSE) + processed_DE_table, + ligand_activities, + processed_condition_markers) + default_prior_table <- prior_table + expect_type(prior_table,"list") + + # Check that columns contain columns from processed_DE_table, processed_expr_table, ligand_activities, and processed_condition_markers + expect_equal(prior_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% select(colnames(processed_DE_table)), + processed_DE_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% mutate(sender = as.character(sender), receiver = as.character(receiver)), + check.attributes = FALSE) + expect_equal(prior_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% select(colnames(processed_expr_table)), + processed_expr_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69", receiver == "CD8 T"), + check.attributes = FALSE) + expect_equal(prior_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% pull(activity), + ligand_activities %>% filter(test_ligand == "Lgals1") %>% pull(aupr_corrected)) + temp_cols <- c("lfc_ligand", "lfc_receptor", "p_val_ligand", "p_val_receptor") + expect_equal(prior_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% select(all_of(paste0(temp_cols, "_group"))), + processed_condition_markers %>% filter(ligand == "Lgals1", receptor == "Cd69") %>% select(all_of(temp_cols)), + check.attributes = FALSE) # Check that prioritization score is the same as the sum of the weighted scores - temp_weights <- prioritizing_weights + temp_weights = c("de_ligand" = 0.5, + "de_receptor" = 0.5, + "activity_scaled" = 1, + "exprs_ligand" = 0.5, + "exprs_receptor" = 0.5, + "ligand_condition_specificity" = 1, + "receptor_condition_specificity" = 1) expect_equal(prior_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% - mutate(prioritization_score = rowSums(across(c(scaled_p_val_ligand_adapted, scaled_p_val_receptor_adapted, scaled_activity, scaled_avg_exprs_ligand, scaled_avg_exprs_receptor, scaled_p_val_ligand_adapted_group, scaled_p_val_receptor_adapted_group)) * temp_weights) / sum(temp_weights)) %>% pull(prioritization_score), + mutate(prioritization_score = rowSums(across(c(scaled_p_val_adapted_ligand, scaled_p_val_adapted_receptor, scaled_activity, scaled_avg_exprs_ligand, scaled_avg_exprs_receptor, scaled_p_val_adapted_ligand_group, scaled_p_val_adapted_receptor_group)) * temp_weights) / sum(temp_weights)) %>% pull(prioritization_score), prior_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% pull(prioritization_score), check.attributes=FALSE) - # Do not pass condition markers - expect_error(generate_prioritization_tables(processed_expr_table, + # Custom weights + prioritizing_weights <- c("de_ligand" = 1, + "de_receptor" = 1, + "activity_scaled" = 2, + "exprs_ligand" = 1, + "exprs_receptor" = 1, + "ligand_condition_specificity" = 0.5, + "receptor_condition_specificity" = 0.5) + prior_table <- generate_prioritization_tables(processed_expr_table, + processed_DE_table, + ligand_activities, + processed_condition_markers, + prioritizing_weights = prioritizing_weights) + temp_weights <- prioritizing_weights + temp_weights[c("de_ligand", "de_receptor", "exprs_ligand", "exprs_receptor")] <- 0.5 + expect_equal(prior_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% + mutate(prioritization_score = rowSums(across(c(scaled_p_val_adapted_ligand, scaled_p_val_adapted_receptor, scaled_activity, scaled_avg_exprs_ligand, scaled_avg_exprs_receptor, scaled_p_val_adapted_ligand_group, scaled_p_val_adapted_receptor_group)) * temp_weights) / sum(temp_weights)) %>% pull(prioritization_score), + prior_table %>% filter(ligand == "Lgals1", sender == "Mono", receptor == "Cd69") %>% pull(prioritization_score), + check.attributes=FALSE) + + # Check "one_condition" scenario + prior_table <- generate_prioritization_tables(processed_expr_table, + processed_DE_table, + ligand_activities, + scenario = "one_condition") + # Check that colnames don't have "_group + expect_false(any(grepl("_group", colnames(prior_table)))) + + # Try giving lr_condition_de + expect_warning(generate_prioritization_tables(processed_expr_table, + processed_DE_table, + ligand_activities, + lr_condition_de = processed_condition_markers, + scenario = "one_condition")) + + # Do not pass condition markers + expect_error(generate_prioritization_tables(processed_expr_table, + processed_DE_table, + ligand_activities, + lr_condition_de = NULL, + scenario = "case_control")) + + # Define priorization weights to 0 except for activity scaled + prioritizing_weights = c("de_ligand" = 0, + "de_receptor" = 0, + "activity_scaled" = 1, + "exprs_ligand" = 0, + "exprs_receptor" = 0, + "ligand_condition_specificity" = 0, + "receptor_condition_specificity" = 0) + + prior_table <- generate_prioritization_tables(processed_expr_table, + processed_DE_table, + ligand_activities, + processed_condition_markers, + prioritizing_weights) + + # Check colnames + expect_equal(colnames(prior_table), + unique(c(colnames(processed_DE_table), colnames(processed_expr_table), + "activity", "rank", "activity_zscore", "scaled_activity", "prioritization_score", "prioritization_rank"))) + + + prior_table_top10 <- prior_table %>% distinct(ligand, prioritization_score) %>% mutate(rank = rank(desc(prioritization_score), ties.method = "average")) %>% arrange(rank, ligand) %>% pull(ligand) %>% .[1:10] + ligands_top10 <- ligand_activities %>% arrange(rank, test_ligand) %>% pull(test_ligand) %>% .[1:10] + + # When using only activity scaled, the top 10 ligands should be the same as the top 10 ligands from the ligand activity table + expect_equal(prior_table_top10, ligands_top10) + + # All weights zero + prioritizing_weights = c("de_ligand" = 0, + "de_receptor" = 0, + "activity_scaled" = 0, + "exprs_ligand" = 0, + "exprs_receptor" = 0, + "ligand_condition_specificity" = 0, + "receptor_condition_specificity" = 0) + + prior_table <- generate_prioritization_tables(processed_expr_table, processed_DE_table, ligand_activities, - lr_condition_de = NULL, - prioritizing_weights)) - - # Define priorization weights to 0 except for activity scaled - prioritizing_weights = c("de_ligand" = 0, - "de_receptor" = 0, - "activity_scaled" = 1, - "exprs_ligand" = 0, - "exprs_receptor" = 0, - "ligand_condition_specificity" = 0, - "receptor_condition_specificity" = 0) - - prior_table <- generate_prioritization_tables(processed_expr_table, - processed_DE_table, - ligand_activities, - processed_condition_markers, - prioritizing_weights) - - # Check colnames - expect_equal(colnames(prior_table), - unique(c(colnames(processed_DE_table), colnames(processed_expr_table), "activity", "rank", "activity_zscore", "scaled_activity", "prioritization_score"))) - - - prior_table_top10 <- prior_table %>% distinct(ligand, prioritization_score) %>% mutate(rank = rank(desc(prioritization_score), ties.method = "average")) %>% arrange(rank, ligand) %>% pull(ligand) %>% .[1:10] - ligands_top10 <- ligand_activities %>% arrange(rank, test_ligand) %>% pull(test_ligand) %>% .[1:10] - - # When using only activity scaled, the top 10 ligands should be the same as the top 10 ligands from the ligand activity table - expect_equal(prior_table_top10, ligands_top10) - - # All weights zero - prioritizing_weights = c("de_ligand" = 0, - "de_receptor" = 0, - "activity_scaled" = 0, - "exprs_ligand" = 0, - "exprs_receptor" = 0, - "ligand_condition_specificity" = 0, - "receptor_condition_specificity" = 0) - - prior_table <- generate_prioritization_tables(processed_expr_table, - processed_DE_table, - ligand_activities, - processed_condition_markers, - prioritizing_weights) - - # Check colnames - expect_equal(colnames(prior_table), - unique(c(colnames(processed_DE_table), colnames(processed_expr_table), "prioritization_score"))) - - prioritizing_weights["de_ligand"] = 1 - prior_table <- generate_prioritization_tables(processed_expr_table, - processed_DE_table, - ligand_activities, - processed_condition_markers, - prioritizing_weights) - - expect_equal(colnames(prior_table), - unique(c(colnames(processed_DE_table), colnames(processed_expr_table), - "lfc_pval_ligand", "p_val_ligand_adapted", "scaled_lfc_ligand", "scaled_p_val_ligand", "scaled_lfc_pval_ligand", "scaled_p_val_ligand_adapted", "prioritization_score"))) + processed_condition_markers, + prioritizing_weights) + + # Check colnames + expect_equal(colnames(prior_table), + unique(c(colnames(processed_DE_table), colnames(processed_expr_table), "prioritization_score", "prioritization_rank"))) + + prioritizing_weights["de_ligand"] = 1 + prior_table <- generate_prioritization_tables(processed_expr_table, + processed_DE_table, + ligand_activities, + processed_condition_markers, + prioritizing_weights) + + + expect_equal(colnames(prior_table), + unique(c(colnames(processed_DE_table), colnames(processed_expr_table), + "lfc_pval_ligand", "p_val_adapted_ligand", "scaled_lfc_ligand", "scaled_p_val_ligand", "scaled_lfc_pval_ligand", "scaled_p_val_adapted_ligand", "prioritization_score", "prioritization_rank"))) + + # Check mushroom plot + usable_colnames <- default_prior_table %>% select(ends_with(c("_ligand", "_receptor", "_sender", "_receiver"))) %>% colnames() %>% + stringr::str_remove("_ligand|_receptor|_sender|_receiver") %>% unique + + mushroom_plots <- list() + for (colname in usable_colnames){ + # If range of values is between 0 and 0.001 + ext <- c("_ligand", "_receptor") + if (colname == "pct_expressed") ext <- c("_sender", "_receiver") + values_range <- range(default_prior_table[, c(paste0(colname, ext[1]), paste0(colname, ext[2]))]) + + # Only scaled values can occupy both size and color + if (values_range[1] >= 0 & values_range[2] <= 1.001){ + mushroom_plots[[colname]] <- make_mushroom_plot(default_prior_table, top_n = 30, size = colname, color = colname) + } else { + # If range of values is not between 0 and 0.001, it should throw an error when used as size + expect_error(make_mushroom_plot(default_prior_table, top_n = 30, size = colname, color = colname)) + + # But ok when used as color + mushroom_plots[[colname]] <- make_mushroom_plot(default_prior_table, top_n = 30, size = "scaled_lfc", color = colname, true_color_range = TRUE) + } } + # Expect all to be ggplot objects + expect_true(all(sapply(mushroom_plots, inherits, "gg"))) + + # If a column name that doesn't exist is passed, it should throw an error + expect_error(make_mushroom_plot(default_prior_table, top_n = 30, size = "non_existent_colname", color = usable_colnames[1])) + expect_error(make_mushroom_plot(default_prior_table, top_n = 30, size = usable_colnames[1], color = "non_existent_colname")) + } + + + }) diff --git a/vignettes/seurat_steps_prioritization.Rmd b/vignettes/seurat_steps_prioritization.Rmd index cf091ab..3b2e77b 100644 --- a/vignettes/seurat_steps_prioritization.Rmd +++ b/vignettes/seurat_steps_prioritization.Rmd @@ -38,7 +38,6 @@ library(nichenetr) # Please update to v2.0.6 library(Seurat) library(SeuratObject) library(tidyverse) - ``` ```{r} @@ -170,7 +169,7 @@ As you can see, the resulting table now show the rankings for *ligand-receptor i We included all columns here, but if you just want relevant columns that were used to calculate the ranking: ```{r} -prior_table %>% select(c('sender', 'receiver', 'ligand', 'receptor', 'scaled_p_val_ligand_adapted', 'scaled_p_val_receptor_adapted', 'scaled_avg_exprs_ligand', 'scaled_avg_exprs_receptor', 'scaled_p_val_ligand_adapted_group', 'scaled_p_val_receptor_adapted_group', 'scaled_activity')) +prior_table %>% select(c('sender', 'receiver', 'ligand', 'receptor', 'scaled_p_val_adapted_ligand', 'scaled_p_val_adapted_receptor', 'scaled_avg_exprs_ligand', 'scaled_avg_exprs_receptor', 'scaled_p_val_adapted_ligand_group', 'scaled_p_val_adapted_receptor_group', 'scaled_activity')) ``` Note that we appended the suffix '_group' to columns that refer to differential expression between conditions, e.g., `lfc_ligand_group` and `lfc_receptor_group.` @@ -231,8 +230,8 @@ prior_table %>% head As NicheNet is a receiver-based pipeline, to prioritize ligand-receptor pairs across multiple receivers, we need to perform the NicheNet analysis for each receiver separately. Let's suppose we want to prioritize ligand-receptor pairs across all T cells (CD4, CD8, and Tregs). The CD8 T analysis has already been performed above. We will use the wrapper function to perform a basic NicheNet analysis on the other two: ```{r} -nichenet_output <- lapply(c("CD4 T", "Treg"), function(receiver_ct){ - nichenet_seuratobj_aggregate(receiver = receiver_ct, +nichenet_outputs <- lapply(c("CD8 T", "CD4 T", "Treg"), function(receiver_ct){ + output <- nichenet_seuratobj_aggregate(receiver = receiver_ct, seurat_obj = seuratObj, condition_colname = "aggregate", condition_oi = condition_oi, @@ -243,42 +242,41 @@ nichenet_output <- lapply(c("CD4 T", "Treg"), function(receiver_ct){ weighted_networks = weighted_networks, expression_pct = 0.05) -}) %>% setNames(c("CD4 T", "Treg")) + # Add receiver cell type in ligand activity table + output$ligand_activities$receiver <- receiver_ct + return(output) +}) ``` To generate the dataframes used for prioritization, we will simply change the `lr_network_filtered` argument to only calculate DE and expression values for ligand-receptor pairs of interest. ```{r} -info_tables2 <- lapply(names(nichenet_output), function(receiver_ct) { - generate_info_tables(seuratObj, - celltype_colname = "celltype", - senders_oi = sender_celltypes, - receivers_oi = receiver_ct, - lr_network_filtered = lr_network %>% - filter(from %in% nichenet_output[[receiver_ct]]$ligand_activities$test_ligand & - to %in% nichenet_output[[receiver_ct]]$background_expressed_genes), - condition_colname = "aggregate", - condition_oi = condition_oi, - condition_reference = condition_reference, - scenario = "case_control") +# Calculate prioritization criteria for each receiver cell type +info_tables <- lapply(nichenet_outputs, function(output) { + lr_network_filtered <- lr_network %>% select(from, to) %>% + filter(from %in% output$ligand_activities$test_ligand & + to %in% output$background_expressed_genes) + + generate_info_tables(seuratObj, + celltype_colname = "celltype", + senders_oi = sender_celltypes, + receivers_oi = unique(output$ligand_activities$receiver), + lr_network_filtered = lr_network_filtered, + condition_colname = "aggregate", + condition_oi = condition_oi, + condition_reference = condition_reference, + scenario = "case_control") }) ``` -We can then combine the results from `generate_info_tables` using `bind_rows`, which will concatenate the rows together. For the ligand activities, we will also add an additional column containing the receiver cell type. Note that for the average expression table (`sender_receiver_info`) and condition specificity (`lr_condition_de`), we need to remove duplicate rows. +We can then combine the results from `generate_info_tables` using `bind_rows`, which will concatenate the rows together. Note that for the average expression table (`sender_receiver_info`) and condition specificity (`lr_condition_de`), we need to remove duplicate rows. ```{r} -# Add CD8 T to list -info_tables2[[3]] <- info_tables - # bind rows of each element of info_tables using pmap -info_tables_combined <- purrr::pmap(info_tables2, bind_rows) - -# Combine ligand activities and add receiver information -ligand_activities_combined <- bind_rows(nichenet_output$`CD4 T`$ligand_activities %>% mutate(receiver = "CD4 T"), - nichenet_output$Treg$ligand_activities %>% mutate(receiver = "Treg"), - ligand_activities %>% mutate(receiver = "CD8 T")) +info_tables_combined <- purrr::pmap(info_tables, bind_rows) +ligand_activities_combined <- purrr::map_dfr(nichenet_outputs, "ligand_activities") prior_table_combined <- generate_prioritization_tables( sender_receiver_info = info_tables_combined$sender_receiver_info %>% distinct, @@ -313,7 +311,8 @@ circos_plot <- make_circos_lr(prior_table_oi, circos_plot ``` -Furthermore, we provide the function `make_mushroom_plot` which allows you to display expression of ligand-receptor pairs in a specific receiver. By default, the fill gradient shows the LFC between cell types, while the size of the semicircle corresponds to the scaled mean expression. You can also choose to show the rankings of each ligand-receptor-sender pair with `show_rankings`, as well as show all data points for context (`show_all_datapoints`). `true_color_range = TRUE` will adjust the limits of the color gradient to the min-max of the values, instead of the limit being from 0 to 1. Note that the numbers displayed here are the rankings within the chosen cell type and not across all receiver cell types (in case of multiple receivers). +Furthermore, we provide the function `make_mushroom_plot` which allows you to display expression of ligand-receptor pairs in a specific receiver. By default, the fill gradient shows the LFC between cell types, while the size of the semicircle corresponds to the scaled mean expression. You can also choose to show the rankings of each ligand-receptor-sender pair with `show_rankings`, as well as show all data points for context (`show_all_datapoints`). `true_color_range = TRUE` will adjust the limits of the color gradient to the min-max of the values, instead of the limit being from 0 to 1. Note that the numbers displayed here are the rankings across all receiver cell types (in case of multiple receivers), and by default the `top_n` ligand-receptor pairs are shown despite the absolute ranking. To show only pairs that have an absolute ranking within top_n across all receivers, set `use_absolute_rank = TRUE`. + ```{r mushroom-plot-1, fig.height=8, fig.width=8} receiver_oi <- "CD8 T" @@ -339,6 +338,7 @@ make_mushroom_plot(prior_table, top_n = 30, size = "pct_expressed", color = "sca axis.title.x = element_text(hjust = 0.25)) ``` + ```{r} sessionInfo() ``` diff --git a/vignettes/seurat_steps_prioritization.md b/vignettes/seurat_steps_prioritization.md index 7a51449..516a814 100644 --- a/vignettes/seurat_steps_prioritization.md +++ b/vignettes/seurat_steps_prioritization.md @@ -168,21 +168,20 @@ info_tables <- generate_info_tables(seuratObj, names(info_tables) ## [1] "sender_receiver_de" "sender_receiver_info" "lr_condition_de" +``` + +``` r info_tables$sender_receiver_de %>% head() -## sender receiver ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor pct_expressed_sender -## 1 DC CD8 T Ccl5 Cxcr3 6.432043 0.16714791 3.299595 1.893317e-25 2.563740e-21 7.758812e-05 1.000000e+00 1.000 -## 2 Mono CD8 T Lyz2 Itgal 5.493265 -0.01687003 2.738198 1.728697e-160 2.340828e-156 4.973381e-02 1.000000e+00 0.933 -## 3 DC CD8 T H2-M2 Cd8a 3.416479 1.94059972 2.678539 1.017174e-272 1.377355e-268 5.250531e-206 7.109745e-202 0.429 -## 4 DC CD8 T Cxcl16 Cxcr6 4.182085 0.54826454 2.365175 1.138617e-243 1.541801e-239 5.987787e-21 8.108063e-17 0.929 -## 5 Mono CD8 T Cxcl9 Cxcr3 4.328801 0.16714791 2.247975 3.834954e-124 5.192911e-120 7.758812e-05 1.000000e+00 0.547 -## 6 Mono CD8 T Cxcl9 Dpp4 4.328801 0.16416445 2.246483 3.834954e-124 5.192911e-120 6.628900e-04 1.000000e+00 0.547 -## pct_expressed_receiver -## 1 0.042 -## 2 0.188 -## 3 0.659 -## 4 0.089 -## 5 0.042 -## 6 0.148 +## sender receiver ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor pct_expressed_sender pct_expressed_receiver +## 1 DC CD8 T H2-M2 Cd8a 11.002412 2.3838066 6.693109 1.017174e-272 1.377355e-268 5.250531e-206 7.109745e-202 0.429 0.659 +## 2 DC CD8 T H2-M2 Klrd1 11.002412 0.9199196 5.961166 1.017174e-272 1.377355e-268 6.104465e-17 8.266056e-13 0.429 0.185 +## 3 DC CD8 T Ccl22 Dpp4 9.920608 0.2991720 5.109890 1.590801e-296 2.154103e-292 6.628900e-04 1.000000e+00 0.500 0.148 +## 4 DC CD8 T Vsig10 Il6st 10.070530 0.1411494 5.105840 2.637179e-194 3.571005e-190 1.470347e-02 1.000000e+00 0.286 0.090 +## 5 DC CD8 T Ccl22 Ccr7 9.920608 0.1468652 5.033737 1.590801e-296 2.154103e-292 5.070025e-05 6.865321e-01 0.500 0.320 +## 6 DC CD8 T Cxcl16 Cxcr6 8.101436 1.8384579 4.969947 1.138617e-243 1.541801e-239 5.987787e-21 8.108063e-17 0.929 0.089 +``` + +``` r info_tables$sender_receiver_info %>% head() ## # A tibble: 6 × 7 ## sender receiver ligand receptor avg_ligand avg_receptor ligand_receptor_prod @@ -193,14 +192,17 @@ info_tables$sender_receiver_info %>% head() ## 4 DC Treg B2m Tap1 216. 7.18 1552. ## 5 Mono Mono B2m Tap1 158. 8.59 1353. ## 6 DC DC B2m Tap1 216. 5.91 1277. +``` + +``` r info_tables$lr_condition_de %>% head() -## ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor -## 1 H2-Ab1 Cd4 2.4021254 0.11569357 1.2589095 4.424390e-06 5.991066e-02 5.634068e-02 1.000000e+00 -## 2 Cxcl10 Dpp4 1.6066163 0.35175421 0.9791853 6.700636e-29 9.073332e-25 1.170731e-06 1.585287e-02 -## 3 B2m Tap1 0.7071427 1.13931050 0.9232266 6.936359e-174 9.392524e-170 3.585450e-52 4.855057e-48 -## 4 H2-T22 Klrd1 1.5223370 -0.05659737 0.7328698 1.006291e-111 1.362618e-107 6.202530e-01 1.000000e+00 -## 5 H2-T23 Klrd1 1.4651999 -0.05659737 0.7043013 1.789643e-114 2.423356e-110 6.202530e-01 1.000000e+00 -## 6 Cxcl10 Cxcr3 1.6066163 -0.25400642 0.6763049 6.700636e-29 9.073332e-25 1.918372e-06 2.597667e-02 +## ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor +## 1 Cxcl11 Dpp4 7.197344 0.7345098 3.965927 0.0001621364 1 1.170731e-06 1.585287e-02 +## 2 Sirpb1c Cd47 6.236414 0.7474147 3.491914 0.0006820290 1 8.720485e-23 1.180841e-18 +## 3 Cxcl11 Cxcr3 7.197344 -1.1317386 3.032803 0.0001621364 1 1.918372e-06 2.597667e-02 +## 4 Ccl22 Dpp4 5.075469 0.7345098 2.904989 0.0863610523 1 1.170731e-06 1.585287e-02 +## 5 F13a1 Itga4 5.436884 0.1228459 2.779865 0.0299628836 1 6.837926e-02 1.000000e+00 +## 6 Vcan Sell 5.234169 0.3254999 2.779835 0.0423593686 1 7.148719e-07 9.680080e-03 ``` Next, we generate the prioritization table. This table contains the @@ -219,22 +221,22 @@ prior_table <- generate_prioritization_tables(info_tables$sender_receiver_info, scenario = "case_control") prior_table %>% head -## # A tibble: 6 × 51 -## sender receiver ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor pct_expressed_sender -## -## 1 NK CD8 T Ptprc Dpp4 0.596 0.164 0.380 2.18e- 7 2.96e- 3 0.000663 1 0.894 -## 2 Mono CD8 T Ptprc Dpp4 0.438 0.164 0.301 3.52e- 5 4.77e- 1 0.000663 1 0.867 -## 3 Mono CD8 T Cxcl10 Dpp4 4.27 0.164 2.22 2.53e- 79 3.43e- 75 0.000663 1 0.867 -## 4 Mono CD8 T Cxcl9 Dpp4 4.33 0.164 2.25 3.83e-124 5.19e-120 0.000663 1 0.547 -## 5 Treg CD8 T Ptprc Dpp4 0.282 0.164 0.223 1.44e- 2 1 e+ 0 0.000663 1 0.685 -## 6 Mono CD8 T Cxcl11 Dpp4 2.36 0.164 1.26 9.28e-121 1.26e-116 0.000663 1 0.307 -## # ℹ 39 more variables: pct_expressed_receiver , avg_ligand , avg_receptor , ligand_receptor_prod , lfc_pval_ligand , -## # p_val_ligand_adapted , scaled_lfc_ligand , scaled_p_val_ligand , scaled_lfc_pval_ligand , scaled_p_val_ligand_adapted , activity , -## # rank , activity_zscore , scaled_activity , lfc_pval_receptor , p_val_receptor_adapted , scaled_lfc_receptor , -## # scaled_p_val_receptor , scaled_lfc_pval_receptor , scaled_p_val_receptor_adapted , scaled_avg_exprs_ligand , -## # scaled_avg_exprs_receptor , lfc_ligand_group , p_val_ligand_group , lfc_pval_ligand_group , p_val_ligand_adapted_group , -## # scaled_lfc_ligand_group , scaled_p_val_ligand_group , scaled_lfc_pval_ligand_group , scaled_p_val_ligand_adapted_group , -## # lfc_receptor_group , p_val_receptor_group , lfc_pval_receptor_group , p_val_receptor_adapted_group , scaled_lfc_receptor_group , … +## # A tibble: 6 × 52 +## sender receiver ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor pct_expressed_sender pct_expressed_receiver avg_ligand +## +## 1 NK CD8 T Ptprc Dpp4 0.642 0.299 0.471 2.18e- 7 2.96e- 3 0.000663 1 0.894 0.148 16.6 +## 2 Mono CD8 T Ptprc Dpp4 0.474 0.299 0.386 3.52e- 5 4.77e- 1 0.000663 1 0.867 0.148 14.9 +## 3 Treg CD8 T Ptprc Dpp4 0.307 0.299 0.303 1.44e- 2 1 e+ 0 0.000663 1 0.685 0.148 13.2 +## 4 Mono CD8 T Cxcl10 Dpp4 4.86 0.299 2.58 2.53e-79 3.43e-75 0.000663 1 0.867 0.148 54.8 +## 5 B CD8 T Ptprc Dpp4 0.201 0.299 0.250 2.08e- 2 1 e+ 0 0.000663 1 0.669 0.148 12.3 +## 6 Mono CD8 T Ebi3 Il6st 4.00 0.141 2.07 9.77e-49 1.32e-44 0.0147 1 0.147 0.09 0.546 +## # ℹ 38 more variables: avg_receptor , ligand_receptor_prod , lfc_pval_ligand , p_val_adapted_ligand , scaled_lfc_ligand , scaled_p_val_ligand , +## # scaled_lfc_pval_ligand , scaled_p_val_adapted_ligand , activity , rank , activity_zscore , scaled_activity , lfc_pval_receptor , +## # p_val_adapted_receptor , scaled_lfc_receptor , scaled_p_val_receptor , scaled_lfc_pval_receptor , scaled_p_val_adapted_receptor , scaled_avg_exprs_ligand , +## # scaled_avg_exprs_receptor , lfc_ligand_group , p_val_ligand_group , lfc_pval_ligand_group , p_val_adapted_ligand_group , scaled_lfc_ligand_group , +## # scaled_p_val_ligand_group , scaled_lfc_pval_ligand_group , scaled_p_val_adapted_ligand_group , lfc_receptor_group , p_val_receptor_group , +## # lfc_pval_receptor_group , p_val_adapted_receptor_group , scaled_lfc_receptor_group , scaled_p_val_receptor_group , scaled_lfc_pval_receptor_group , +## # scaled_p_val_adapted_receptor_group , prioritization_score , prioritization_rank ``` As you can see, the resulting table now show the rankings for @@ -248,23 +250,22 @@ We included all columns here, but if you just want relevant columns that were used to calculate the ranking: ``` r -prior_table %>% select(c('sender', 'receiver', 'ligand', 'receptor', 'scaled_p_val_ligand_adapted', 'scaled_p_val_receptor_adapted', 'scaled_avg_exprs_ligand', 'scaled_avg_exprs_receptor', 'scaled_p_val_ligand_adapted_group', 'scaled_p_val_receptor_adapted_group', 'scaled_activity')) -## # A tibble: 1,272 × 11 -## sender receiver ligand receptor scaled_p_val_ligand_adapted scaled_p_val_receptor_adapted scaled_avg_exprs_ligand scaled_avg_exprs_receptor scaled_p_val_ligand_…¹ -## -## 1 NK CD8 T Ptprc Dpp4 0.869 0.829 1.00 1.00 0.850 -## 2 Mono CD8 T Ptprc Dpp4 0.841 0.829 0.867 1.00 0.850 -## 3 Mono CD8 T Cxcl10 Dpp4 0.960 0.829 1.00 1.00 0.929 -## 4 Mono CD8 T Cxcl9 Dpp4 0.975 0.829 1.00 1.00 0.787 -## 5 Treg CD8 T Ptprc Dpp4 0.756 0.829 0.741 1.00 0.850 -## 6 Mono CD8 T Cxcl11 Dpp4 0.973 0.829 1.00 1.00 0.732 -## 7 B CD8 T Ptprc Dpp4 0.748 0.829 0.666 1.00 0.850 -## 8 DC CD8 T Icam1 Il2rg 0.876 0.714 1.00 0.995 0.717 -## 9 DC CD8 T Ccl22 Dpp4 0.997 0.829 1.00 1.00 0.539 -## 10 NK CD8 T Cd320 Jaml 0.889 0.943 0.905 1.00 0.472 -## # ℹ 1,262 more rows -## # ℹ abbreviated name: ¹​scaled_p_val_ligand_adapted_group -## # ℹ 2 more variables: scaled_p_val_receptor_adapted_group , scaled_activity +prior_table %>% select(c('sender', 'receiver', 'ligand', 'receptor', 'scaled_p_val_adapted_ligand', 'scaled_p_val_adapted_receptor', 'scaled_avg_exprs_ligand', 'scaled_avg_exprs_receptor', 'scaled_p_val_adapted_ligand_group', 'scaled_p_val_adapted_receptor_group', 'scaled_activity')) +## # A tibble: 1,212 × 11 +## sender receiver ligand receptor scaled_p_val_adapted_ligand scaled_p_val_adapted_re…¹ scaled_avg_exprs_lig…² scaled_avg_exprs_rec…³ scaled_p_val_adapted…⁴ scaled_p_val_adapted…⁵ scaled_activity +## +## 1 NK CD8 T Ptprc Dpp4 0.871 0.846 1.00 1.00 0.844 0.833 0.660 +## 2 Mono CD8 T Ptprc Dpp4 0.841 0.846 0.867 1.00 0.844 0.833 0.660 +## 3 Treg CD8 T Ptprc Dpp4 0.754 0.846 0.741 1.00 0.844 0.833 0.660 +## 4 Mono CD8 T Cxcl10 Dpp4 0.958 0.846 1.00 1.00 0.926 0.833 0.309 +## 5 B CD8 T Ptprc Dpp4 0.747 0.846 0.666 1.00 0.844 0.833 0.660 +## 6 Mono CD8 T Ebi3 Il6st 0.937 0.692 1.00 0.635 0.680 0.515 1.00 +## 7 Mono CD8 T Cxcl9 Dpp4 0.974 0.846 1.00 1.00 0.779 0.833 0.263 +## 8 DC CD8 T Icam1 Il2rg 0.878 0.723 1.00 0.995 0.713 0.985 0.273 +## 9 DC CD8 T B2m Tap2 0.862 0.631 1.00 0.781 1 0.864 0.254 +## 10 Mono CD8 T Cxcl11 Dpp4 0.972 0.846 1.00 1.00 0.721 0.833 0.273 +## # ℹ 1,202 more rows +## # ℹ abbreviated names: ¹​scaled_p_val_adapted_receptor, ²​scaled_avg_exprs_ligand, ³​scaled_avg_exprs_receptor, ⁴​scaled_p_val_adapted_ligand_group, ⁵​scaled_p_val_adapted_receptor_group ``` Note that we appended the suffix ’\_group’ to columns that refer to @@ -328,22 +329,22 @@ prior_table <- generate_prioritization_tables(processed_expr_table, prioritizing_weights) prior_table %>% head -## # A tibble: 6 × 51 -## sender receiver ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor pct_expressed_sender -## -## 1 NK CD8 T Ptprc Dpp4 0.596 0.164 0.380 2.18e- 7 2.96e- 3 0.000663 1 0.894 -## 2 Mono CD8 T Ptprc Dpp4 0.438 0.164 0.301 3.52e- 5 4.77e- 1 0.000663 1 0.867 -## 3 Mono CD8 T Cxcl10 Dpp4 4.27 0.164 2.22 2.53e- 79 3.43e- 75 0.000663 1 0.867 -## 4 Mono CD8 T Cxcl9 Dpp4 4.33 0.164 2.25 3.83e-124 5.19e-120 0.000663 1 0.547 -## 5 Treg CD8 T Ptprc Dpp4 0.282 0.164 0.223 1.44e- 2 1 e+ 0 0.000663 1 0.685 -## 6 Mono CD8 T Cxcl11 Dpp4 2.36 0.164 1.26 9.28e-121 1.26e-116 0.000663 1 0.307 -## # ℹ 39 more variables: pct_expressed_receiver , avg_ligand , avg_receptor , ligand_receptor_prod , lfc_pval_ligand , -## # p_val_ligand_adapted , scaled_lfc_ligand , scaled_p_val_ligand , scaled_lfc_pval_ligand , scaled_p_val_ligand_adapted , activity , -## # rank , activity_zscore , scaled_activity , lfc_pval_receptor , p_val_receptor_adapted , scaled_lfc_receptor , -## # scaled_p_val_receptor , scaled_lfc_pval_receptor , scaled_p_val_receptor_adapted , scaled_avg_exprs_ligand , -## # scaled_avg_exprs_receptor , lfc_ligand_group , p_val_ligand_group , lfc_pval_ligand_group , p_val_ligand_adapted_group , -## # scaled_lfc_ligand_group , scaled_p_val_ligand_group , scaled_lfc_pval_ligand_group , scaled_p_val_ligand_adapted_group , -## # lfc_receptor_group , p_val_receptor_group , lfc_pval_receptor_group , p_val_receptor_adapted_group , scaled_lfc_receptor_group , … +## # A tibble: 6 × 52 +## sender receiver ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor pct_expressed_sender pct_expressed_receiver avg_ligand +## +## 1 NK CD8 T Ptprc Dpp4 0.642 0.299 0.471 2.18e- 7 2.96e- 3 0.000663 1 0.894 0.148 16.6 +## 2 Mono CD8 T Ptprc Dpp4 0.474 0.299 0.386 3.52e- 5 4.77e- 1 0.000663 1 0.867 0.148 14.9 +## 3 Treg CD8 T Ptprc Dpp4 0.307 0.299 0.303 1.44e- 2 1 e+ 0 0.000663 1 0.685 0.148 13.2 +## 4 Mono CD8 T Cxcl10 Dpp4 4.86 0.299 2.58 2.53e-79 3.43e-75 0.000663 1 0.867 0.148 54.8 +## 5 B CD8 T Ptprc Dpp4 0.201 0.299 0.250 2.08e- 2 1 e+ 0 0.000663 1 0.669 0.148 12.3 +## 6 Mono CD8 T Ebi3 Il6st 4.00 0.141 2.07 9.77e-49 1.32e-44 0.0147 1 0.147 0.09 0.546 +## # ℹ 38 more variables: avg_receptor , ligand_receptor_prod , lfc_pval_ligand , p_val_adapted_ligand , scaled_lfc_ligand , scaled_p_val_ligand , +## # scaled_lfc_pval_ligand , scaled_p_val_adapted_ligand , activity , rank , activity_zscore , scaled_activity , lfc_pval_receptor , +## # p_val_adapted_receptor , scaled_lfc_receptor , scaled_p_val_receptor , scaled_lfc_pval_receptor , scaled_p_val_adapted_receptor , scaled_avg_exprs_ligand , +## # scaled_avg_exprs_receptor , lfc_ligand_group , p_val_ligand_group , lfc_pval_ligand_group , p_val_adapted_ligand_group , scaled_lfc_ligand_group , +## # scaled_p_val_ligand_group , scaled_lfc_pval_ligand_group , scaled_p_val_adapted_ligand_group , lfc_receptor_group , p_val_receptor_group , +## # lfc_pval_receptor_group , p_val_adapted_receptor_group , scaled_lfc_receptor_group , scaled_p_val_receptor_group , scaled_lfc_pval_receptor_group , +## # scaled_p_val_adapted_receptor_group , prioritization_score , prioritization_rank ``` # Prioritizing across multiple receivers @@ -357,8 +358,8 @@ the wrapper function to perform a basic NicheNet analysis on the other two: ``` r -nichenet_output <- lapply(c("CD4 T", "Treg"), function(receiver_ct){ - nichenet_seuratobj_aggregate(receiver = receiver_ct, +nichenet_outputs <- lapply(c("CD8 T", "CD4 T", "Treg"), function(receiver_ct){ + output <- nichenet_seuratobj_aggregate(receiver = receiver_ct, seurat_obj = seuratObj, condition_colname = "aggregate", condition_oi = condition_oi, @@ -369,7 +370,11 @@ nichenet_output <- lapply(c("CD4 T", "Treg"), function(receiver_ct){ weighted_networks = weighted_networks, expression_pct = 0.05) -}) %>% setNames(c("CD4 T", "Treg")) + # Add receiver cell type in ligand activity table + output$ligand_activities$receiver <- receiver_ct + return(output) +}) +## [1] "The RNA assay will be used for the analysis." ## [1] "Read in and process NicheNet's networks" ## [1] "Define expressed ligands and receptors in receiver and sender cells" ## [1] "Perform DE analysis in receiver cell" @@ -377,6 +382,15 @@ nichenet_output <- lapply(c("CD4 T", "Treg"), function(receiver_ct){ ## [1] "Infer active target genes of the prioritized ligands" ## [1] "Infer receptors of the prioritized ligands" ## [1] "Perform DE analysis in sender cells" +## [1] "The RNA assay will be used for the analysis." +## [1] "Read in and process NicheNet's networks" +## [1] "Define expressed ligands and receptors in receiver and sender cells" +## [1] "Perform DE analysis in receiver cell" +## [1] "Perform NicheNet ligand activity analysis" +## [1] "Infer active target genes of the prioritized ligands" +## [1] "Infer receptors of the prioritized ligands" +## [1] "Perform DE analysis in sender cells" +## [1] "The RNA assay will be used for the analysis." ## [1] "Read in and process NicheNet's networks" ## [1] "Define expressed ligands and receptors in receiver and sender cells" ## [1] "Perform DE analysis in receiver cell" @@ -391,39 +405,33 @@ change the `lr_network_filtered` argument to only calculate DE and expression values for ligand-receptor pairs of interest. ``` r -info_tables2 <- lapply(names(nichenet_output), function(receiver_ct) { - generate_info_tables(seuratObj, - celltype_colname = "celltype", - senders_oi = sender_celltypes, - receivers_oi = receiver_ct, - lr_network_filtered = lr_network %>% - filter(from %in% nichenet_output[[receiver_ct]]$ligand_activities$test_ligand & - to %in% nichenet_output[[receiver_ct]]$background_expressed_genes), - condition_colname = "aggregate", - condition_oi = condition_oi, - condition_reference = condition_reference, - scenario = "case_control") +# Calculate prioritization criteria for each receiver cell type +info_tables <- lapply(nichenet_outputs, function(output) { + lr_network_filtered <- lr_network %>% select(from, to) %>% + filter(from %in% output$ligand_activities$test_ligand & + to %in% output$background_expressed_genes) + + generate_info_tables(seuratObj, + celltype_colname = "celltype", + senders_oi = sender_celltypes, + receivers_oi = unique(output$ligand_activities$receiver), + lr_network_filtered = lr_network_filtered, + condition_colname = "aggregate", + condition_oi = condition_oi, + condition_reference = condition_reference, + scenario = "case_control") }) ``` We can then combine the results from `generate_info_tables` using -`bind_rows`, which will concatenate the rows together. For the ligand -activities, we will also add an additional column containing the -receiver cell type. Note that for the average expression table -(`sender_receiver_info`) and condition specificity (`lr_condition_de`), -we need to remove duplicate rows. +`bind_rows`, which will concatenate the rows together. Note that for the +average expression table (`sender_receiver_info`) and condition +specificity (`lr_condition_de`), we need to remove duplicate rows. ``` r -# Add CD8 T to list -info_tables2[[3]] <- info_tables - # bind rows of each element of info_tables using pmap -info_tables_combined <- purrr::pmap(info_tables2, bind_rows) - -# Combine ligand activities and add receiver information -ligand_activities_combined <- bind_rows(nichenet_output$`CD4 T`$ligand_activities %>% mutate(receiver = "CD4 T"), - nichenet_output$Treg$ligand_activities %>% mutate(receiver = "Treg"), - ligand_activities %>% mutate(receiver = "CD8 T")) +info_tables_combined <- purrr::pmap(info_tables, bind_rows) +ligand_activities_combined <- purrr::map_dfr(nichenet_outputs, "ligand_activities") prior_table_combined <- generate_prioritization_tables( sender_receiver_info = info_tables_combined$sender_receiver_info %>% distinct, @@ -433,22 +441,22 @@ prior_table_combined <- generate_prioritization_tables( scenario = "case_control") head(prior_table_combined) -## # A tibble: 6 × 51 -## sender receiver ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor pct_expressed_sender -## -## 1 NK CD8 T Ptprc Dpp4 0.596 0.164 0.380 0.000000218 0.00296 6.63e- 4 1 e+ 0 0.894 -## 2 NK CD4 T Ptprc Cd4 0.596 0.996 0.796 0.000000218 0.00296 2.63e-34 3.56e-30 0.894 -## 3 B CD4 T H2-Eb1 Cd4 4.02 0.996 2.51 0 0 2.63e-34 3.56e-30 0.93 -## 4 Mono CD8 T Ptprc Dpp4 0.438 0.164 0.301 0.0000352 0.477 6.63e- 4 1 e+ 0 0.867 -## 5 Mono CD4 T Ptprc Cd4 0.438 0.996 0.717 0.0000352 0.477 2.63e-34 3.56e-30 0.867 -## 6 NK CD4 T Ptprc Cd247 0.596 0.457 0.526 0.000000218 0.00296 5.61e- 4 1 e+ 0 0.894 -## # ℹ 39 more variables: pct_expressed_receiver , avg_ligand , avg_receptor , ligand_receptor_prod , lfc_pval_ligand , -## # p_val_ligand_adapted , scaled_lfc_ligand , scaled_p_val_ligand , scaled_lfc_pval_ligand , scaled_p_val_ligand_adapted , activity , -## # rank , activity_zscore , scaled_activity , lfc_pval_receptor , p_val_receptor_adapted , scaled_lfc_receptor , -## # scaled_p_val_receptor , scaled_lfc_pval_receptor , scaled_p_val_receptor_adapted , scaled_avg_exprs_ligand , -## # scaled_avg_exprs_receptor , lfc_ligand_group , p_val_ligand_group , lfc_pval_ligand_group , p_val_ligand_adapted_group , -## # scaled_lfc_ligand_group , scaled_p_val_ligand_group , scaled_lfc_pval_ligand_group , scaled_p_val_ligand_adapted_group , -## # lfc_receptor_group , p_val_receptor_group , lfc_pval_receptor_group , p_val_receptor_adapted_group , scaled_lfc_receptor_group , … +## # A tibble: 6 × 52 +## sender receiver ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor pct_expressed_sender pct_expressed_receiver avg_ligand +## +## 1 NK CD8 T Ptprc Dpp4 0.642 0.299 0.471 0.000000218 0.00296 6.63e- 4 1 e+ 0 0.894 0.148 16.6 +## 2 NK CD4 T Ptprc Cd4 0.642 1.71 1.18 0.000000218 0.00296 2.63e-34 3.56e-30 0.894 0.226 16.6 +## 3 Mono CD8 T Ptprc Dpp4 0.474 0.299 0.386 0.0000352 0.477 6.63e- 4 1 e+ 0 0.867 0.148 14.9 +## 4 B CD4 T H2-Eb1 Cd4 5.00 1.71 3.36 0 0 2.63e-34 3.56e-30 0.93 0.226 31.0 +## 5 Mono CD4 T Ptprc Cd4 0.474 1.71 1.09 0.0000352 0.477 2.63e-34 3.56e-30 0.867 0.226 14.9 +## 6 NK CD4 T Ptprc Cd247 0.642 0.599 0.620 0.000000218 0.00296 5.61e- 4 1 e+ 0 0.894 0.309 16.6 +## # ℹ 38 more variables: avg_receptor , ligand_receptor_prod , lfc_pval_ligand , p_val_adapted_ligand , scaled_lfc_ligand , scaled_p_val_ligand , +## # scaled_lfc_pval_ligand , scaled_p_val_adapted_ligand , activity , rank , activity_zscore , scaled_activity , lfc_pval_receptor , +## # p_val_adapted_receptor , scaled_lfc_receptor , scaled_p_val_receptor , scaled_lfc_pval_receptor , scaled_p_val_adapted_receptor , scaled_avg_exprs_ligand , +## # scaled_avg_exprs_receptor , lfc_ligand_group , p_val_ligand_group , lfc_pval_ligand_group , p_val_adapted_ligand_group , scaled_lfc_ligand_group , +## # scaled_p_val_ligand_group , scaled_lfc_pval_ligand_group , scaled_p_val_adapted_ligand_group , lfc_receptor_group , p_val_receptor_group , +## # lfc_pval_receptor_group , p_val_adapted_receptor_group , scaled_lfc_receptor_group , scaled_p_val_receptor_group , scaled_lfc_pval_receptor_group , +## # scaled_p_val_adapted_receptor_group , prioritization_score , prioritization_rank ``` ### Extra visualization of ligand-receptor pairs @@ -488,9 +496,11 @@ ligand-receptor-sender pair with `show_rankings`, as well as show all data points for context (`show_all_datapoints`). `true_color_range = TRUE` will adjust the limits of the color gradient to the min-max of the values, instead of the limit being from 0 to 1. -Note that the numbers displayed here are the rankings within the chosen -cell type and not across all receiver cell types (in case of multiple -receivers). +Note that the numbers displayed here are the rankings across all +receiver cell types (in case of multiple receivers), and by default the +`top_n` ligand-receptor pairs are shown despite the absolute ranking. To +show only pairs that have an absolute ranking within top_n across all +receivers, set `use_absolute_rank = TRUE`. ``` r receiver_oi <- "CD8 T" @@ -513,7 +523,10 @@ columns from the prioritization table (those with the `_ligand` or ``` r print(paste0("Column names that you can use are: ", paste0(prior_table %>% select(ends_with(c("_ligand", "_receptor", "_sender", "_receiver"))) %>% colnames() %>% str_remove("_ligand|_receptor|_sender|_receiver") %>% unique, collapse = ", "))) -## [1] "Column names that you can use are: lfc, p_val, p_adj, avg, lfc_pval, scaled_lfc, scaled_p_val, scaled_lfc_pval, scaled_avg_exprs, pct_expressed" +## [1] "Column names that you can use are: lfc, p_val, p_adj, avg, lfc_pval, p_val_adapted, scaled_lfc, scaled_p_val, scaled_lfc_pval, scaled_p_val_adapted, scaled_avg_exprs, pct_expressed" +``` + +``` r # Change size and color columns make_mushroom_plot(prior_table, top_n = 30, size = "pct_expressed", color = "scaled_avg_exprs") + @@ -525,64 +538,55 @@ make_mushroom_plot(prior_table, top_n = 30, size = "pct_expressed", color = "sca ``` r sessionInfo() -## R version 4.3.2 (2023-10-31) -## Platform: x86_64-redhat-linux-gnu (64-bit) +## R version 4.3.3 (2024-02-29) +## Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: CentOS Stream 8 ## ## Matrix products: default -## BLAS/LAPACK: /usr/lib64/libopenblaso-r0.3.15.so; LAPACK version 3.9.0 +## BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so; LAPACK version 3.9.0 ## ## locale: -## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 -## [6] LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C -## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C +## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 +## [8] LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C ## -## time zone: Asia/Bangkok +## time zone: Europe/Brussels ## tzcode source: system (glibc) ## ## attached base packages: ## [1] stats graphics grDevices utils datasets methods base ## ## other attached packages: -## [1] forcats_1.0.0 stringr_1.5.0 dplyr_1.1.4 purrr_1.0.2 readr_2.1.2 tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.4 -## [9] tidyverse_1.3.1 SeuratObject_5.0.1 Seurat_4.4.0 nichenetr_2.1.0 +## [1] Seurat_5.1.0 SeuratObject_5.0.2 sp_2.1-4 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 +## [11] tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0 nichenetr_2.2.0 ## ## loaded via a namespace (and not attached): -## [1] fs_1.6.3 matrixStats_1.2.0 spatstat.sparse_3.0-3 bitops_1.0-7 lubridate_1.9.3 httr_1.4.7 -## [7] RColorBrewer_1.1-3 doParallel_1.0.17 tools_4.3.2 sctransform_0.4.0 backports_1.4.1 utf8_1.2.4 -## [13] R6_2.5.1 lazyeval_0.2.2 uwot_0.1.16 GetoptLong_1.0.5 withr_2.5.2 sp_2.1-2 -## [19] gridExtra_2.3 fdrtool_1.2.17 progressr_0.14.0 cli_3.6.2 spatstat.explore_3.2-1 labeling_0.4.3 -## [25] spatstat.data_3.0-3 randomForest_4.7-1.1 proxy_0.4-27 ggridges_0.5.5 pbapply_1.7-2 foreign_0.8-85 -## [31] parallelly_1.36.0 limma_3.56.2 readxl_1.4.3 rstudioapi_0.15.0 gridGraphics_0.5-1 visNetwork_2.1.2 -## [37] generics_0.1.3 shape_1.4.6 ica_1.0-3 spatstat.random_3.2-2 car_3.1-2 Matrix_1.6-4 -## [43] fansi_1.0.6 S4Vectors_0.38.1 abind_1.4-5 lifecycle_1.0.4 yaml_2.3.8 carData_3.0-5 -## [49] recipes_1.0.7 Rtsne_0.17 grid_4.3.2 promises_1.2.1 crayon_1.5.2 miniUI_0.1.1.1 -## [55] lattice_0.21-9 haven_2.4.3 cowplot_1.1.2 pillar_1.9.0 knitr_1.45 ComplexHeatmap_2.16.0 -## [61] rjson_0.2.21 future.apply_1.11.0 codetools_0.2-19 leiden_0.3.9 glue_1.6.2 data.table_1.14.10 -## [67] vctrs_0.6.5 png_0.1-8 spam_2.10-0 cellranger_1.1.0 gtable_0.3.4 assertthat_0.2.1 -## [73] gower_1.0.1 xfun_0.41 mime_0.12 prodlim_2023.08.28 survival_3.5-7 timeDate_4032.109 -## [79] iterators_1.0.14 hardhat_1.3.0 lava_1.7.3 DiagrammeR_1.0.10 ellipsis_0.3.2 fitdistrplus_1.1-11 -## [85] ROCR_1.0-11 ipred_0.9-14 nlme_3.1-163 RcppAnnoy_0.0.21 irlba_2.3.5.1 KernSmooth_2.23-22 -## [91] rpart_4.1.21 colorspace_2.1-0 BiocGenerics_0.46.0 DBI_1.1.3 Hmisc_5.1-0 nnet_7.3-19 -## [97] tidyselect_1.2.0 compiler_4.3.2 rvest_1.0.2 htmlTable_2.4.1 xml2_1.3.6 plotly_4.10.0 -## [103] shadowtext_0.1.2 checkmate_2.3.1 scales_1.3.0 caTools_1.18.2 lmtest_0.9-40 digest_0.6.33 -## [109] goftest_1.2-3 spatstat.utils_3.0-4 rmarkdown_2.11 htmltools_0.5.7 pkgconfig_2.0.3 base64enc_0.1-3 -## [115] highr_0.10 dbplyr_2.1.1 fastmap_1.1.1 rlang_1.1.2 GlobalOptions_0.1.2 htmlwidgets_1.6.2 -## [121] shiny_1.7.1 farver_2.1.1 zoo_1.8-12 jsonlite_1.8.8 ModelMetrics_1.2.2.2 magrittr_2.0.3 -## [127] Formula_1.2-5 dotCall64_1.1-1 patchwork_1.1.3 munsell_0.5.0 Rcpp_1.0.11 ggnewscale_0.4.9 -## [133] reticulate_1.34.0 stringi_1.7.6 pROC_1.18.5 MASS_7.3-60 plyr_1.8.9 parallel_4.3.2 -## [139] listenv_0.9.0 ggrepel_0.9.4 deldir_2.0-2 splines_4.3.2 tensor_1.5 hms_1.1.3 -## [145] circlize_0.4.15 igraph_1.2.11 ggpubr_0.6.0 spatstat.geom_3.2-7 ggsignif_0.6.4 reshape2_1.4.4 -## [151] stats4_4.3.2 reprex_2.0.1 evaluate_0.23 modelr_0.1.8 tzdb_0.4.0 foreach_1.5.2 -## [157] tweenr_2.0.2 httpuv_1.6.13 RANN_2.6.1 polyclip_1.10-6 future_1.33.0 clue_0.3-64 -## [163] scattermore_1.2 ggforce_0.4.1 broom_0.7.12 xtable_1.8-4 e1071_1.7-14 rstatix_0.7.2 -## [169] later_1.3.2 viridisLite_0.4.2 class_7.3-22 IRanges_2.34.1 cluster_2.1.4 timechange_0.2.0 -## [175] globals_0.16.2 caret_6.0-94 +## [1] RcppAnnoy_0.0.22 splines_4.3.3 later_1.3.2 bitops_1.0-7 polyclip_1.10-6 hardhat_1.3.1 pROC_1.18.5 rpart_4.1.23 +## [9] fastDummies_1.7.3 lifecycle_1.0.4 rstatix_0.7.2 doParallel_1.0.17 globals_0.16.3 lattice_0.22-5 MASS_7.3-60.0.1 backports_1.4.1 +## [17] magrittr_2.0.3 limma_3.58.1 rmarkdown_2.27 Hmisc_5.1-2 plotly_4.10.4 yaml_2.3.8 httpuv_1.6.15 sctransform_0.4.1 +## [25] spam_2.10-0 spatstat.sparse_3.0-3 reticulate_1.37.0 cowplot_1.1.3 pbapply_1.7-2 RColorBrewer_1.1-3 abind_1.4-5 Rtsne_0.17 +## [33] presto_1.0.0 BiocGenerics_0.48.1 nnet_7.3-19 tweenr_2.0.3 ipred_0.9-14 circlize_0.4.16 lava_1.8.0 IRanges_2.36.0 +## [41] S4Vectors_0.40.2 ggrepel_0.9.5 irlba_2.3.5.1 listenv_0.9.1 spatstat.utils_3.0-4 goftest_1.2-3 RSpectra_0.16-1 spatstat.random_3.2-3 +## [49] fitdistrplus_1.1-11 parallelly_1.37.1 leiden_0.4.3.1 codetools_0.2-19 ggforce_0.4.2 tidyselect_1.2.1 shape_1.4.6.1 farver_2.1.2 +## [57] matrixStats_1.3.0 stats4_4.3.3 base64enc_0.1-3 spatstat.explore_3.2-7 jsonlite_1.8.8 caret_6.0-94 GetoptLong_1.0.5 e1071_1.7-14 +## [65] progressr_0.14.0 Formula_1.2-5 ggridges_0.5.6 survival_3.5-8 iterators_1.0.14 foreach_1.5.2 tools_4.3.3 ggnewscale_0.4.10 +## [73] ica_1.0-3 Rcpp_1.0.12 glue_1.7.0 prodlim_2023.08.28 gridExtra_2.3 xfun_0.44 withr_3.0.0 fastmap_1.2.0 +## [81] fansi_1.0.6 caTools_1.18.2 digest_0.6.35 gridGraphics_0.5-1 timechange_0.3.0 R6_2.5.1 mime_0.12 colorspace_2.1-0 +## [89] scattermore_1.2 tensor_1.5 spatstat.data_3.0-4 DiagrammeR_1.0.11 utf8_1.2.4 generics_0.1.3 data.table_1.15.4 recipes_1.0.10 +## [97] class_7.3-22 httr_1.4.7 htmlwidgets_1.6.4 uwot_0.2.2 ModelMetrics_1.2.2.2 pkgconfig_2.0.3 gtable_0.3.5 timeDate_4032.109 +## [105] ComplexHeatmap_2.18.0 lmtest_0.9-40 shadowtext_0.1.3 htmltools_0.5.8.1 carData_3.0-5 dotCall64_1.1-1 clue_0.3-65 scales_1.3.0 +## [113] png_0.1-8 gower_1.0.1 knitr_1.46 rstudioapi_0.16.0 tzdb_0.4.0 reshape2_1.4.4 rjson_0.2.21 checkmate_2.3.1 +## [121] visNetwork_2.1.2 nlme_3.1-164 proxy_0.4-27 zoo_1.8-12 GlobalOptions_0.1.2 KernSmooth_2.23-22 parallel_4.3.3 miniUI_0.1.1.1 +## [129] foreign_0.8-86 pillar_1.9.0 grid_4.3.3 vctrs_0.6.5 RANN_2.6.1 ggpubr_0.6.0 randomForest_4.7-1.1 promises_1.3.0 +## [137] car_3.1-2 xtable_1.8-4 cluster_2.1.6 htmlTable_2.4.2 evaluate_0.23 cli_3.6.2 compiler_4.3.3 rlang_1.1.3 +## [145] crayon_1.5.2 ggsignif_0.6.4 future.apply_1.11.2 labeling_0.4.3 fdrtool_1.2.17 plyr_1.8.9 stringi_1.8.4 viridisLite_0.4.2 +## [153] deldir_2.0-4 munsell_0.5.1 lazyeval_0.2.2 spatstat.geom_3.2-9 Matrix_1.6-5 RcppHNSW_0.6.0 hms_1.1.3 patchwork_1.2.0 +## [161] future_1.33.2 statmod_1.5.0 shiny_1.8.1.1 highr_0.10 ROCR_1.0-11 broom_1.0.5 igraph_2.0.3 ``` ### References -
+
diff --git a/vignettes/seurat_steps_prioritization_files/figure-gfm/lr-circos-1.png b/vignettes/seurat_steps_prioritization_files/figure-gfm/lr-circos-1.png index b2c79b5..400a49c 100644 Binary files a/vignettes/seurat_steps_prioritization_files/figure-gfm/lr-circos-1.png and b/vignettes/seurat_steps_prioritization_files/figure-gfm/lr-circos-1.png differ diff --git a/vignettes/seurat_steps_prioritization_files/figure-gfm/lr-circos-unused-1.png b/vignettes/seurat_steps_prioritization_files/figure-gfm/lr-circos-unused-1.png index bcbfffb..3972427 100644 Binary files a/vignettes/seurat_steps_prioritization_files/figure-gfm/lr-circos-unused-1.png and b/vignettes/seurat_steps_prioritization_files/figure-gfm/lr-circos-unused-1.png differ diff --git a/vignettes/seurat_steps_prioritization_files/figure-gfm/mushroom-plot-1-1.png b/vignettes/seurat_steps_prioritization_files/figure-gfm/mushroom-plot-1-1.png index 976c32c..5f3784c 100644 Binary files a/vignettes/seurat_steps_prioritization_files/figure-gfm/mushroom-plot-1-1.png and b/vignettes/seurat_steps_prioritization_files/figure-gfm/mushroom-plot-1-1.png differ diff --git a/vignettes/seurat_steps_prioritization_files/figure-gfm/mushroom-plot-2-1.png b/vignettes/seurat_steps_prioritization_files/figure-gfm/mushroom-plot-2-1.png index c662f46..1c6cb73 100644 Binary files a/vignettes/seurat_steps_prioritization_files/figure-gfm/mushroom-plot-2-1.png and b/vignettes/seurat_steps_prioritization_files/figure-gfm/mushroom-plot-2-1.png differ