diff --git a/DESCRIPTION b/DESCRIPTION index f5a9ae2..aa55f81 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,15 +1,12 @@ Package: multinichenetr Type: Package -Title: Differential Cell-Cell Communication analysis for scRNAseq data with complex multi-sample multi-group designs +Title: MultiNicheNet: a flexible framework for differential cell-cell communication analysis from multi-sample multi-condition single-cell transcriptomics data Version: 2.0.0 Author: person("Robin", "Browaeys", email = "robin.browaeys@ugent.be", role = c("aut", "cre")) Maintainer: Robin Browaeys -Description: This package allows you the investigate intercellular communication from a computational perspective. - It is an extension of the NicheNet framework (https://github.com/saeyslab/nichenetr) to better suit scRNAseq datasets with complex designs. - These datasets can contain multiple samples (e.g. patients) over different groups of interest (e.g. disease subtypes). - With MultiNicheNet, you can now better analyze the differences in cell-cell signaling between the different groups of interest. - MultiNicheNet will give a you a list of the ligand-receptor interactions that are most strongly differentially expressed between patients of the different groups, and also most active in the different groups as given by the NicheNet ligand activity. +Description: This package allows you the investigate differences in intercellular communication between multiple conditions of interest. It is a flexible framework that builds upon the NicheNet framework (https://github.com/saeyslab/nichenetr) to better suit scRNAseq datasets with complex designs. + These datasets can contain multiple samples (e.g. patients) over different groups of interest (e.g. disease subtypes). With MultiNicheNet, you can now better analyze the differences in cell-cell signaling between the different groups of interest. License: GPL-3 + file LICENSE Encoding: UTF-8 URL: https://github.com/browaeysrobin/multinichenetr @@ -59,7 +56,11 @@ Suggests: knitr, testthat, covr, - tidyverse + tidyverse, + BiocStyle, + rmarkdown +VignetteBuilder: + knitr Remotes: github::saeyslab/nichenetr RoxygenNote: 7.2.3 diff --git a/NAMESPACE b/NAMESPACE index b8617e2..99b5d9b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -55,6 +55,7 @@ export(prioritize_condition_specific_receiver) export(prioritize_condition_specific_sender) export(process_abund_info) export(process_abundance_expression_info) +export(process_geneset_data) export(process_info_to_ic) export(visualize_network) import(circlize) diff --git a/R/ligand_activities.R b/R/ligand_activities.R index bca8670..e309ab7 100644 --- a/R/ligand_activities.R +++ b/R/ligand_activities.R @@ -547,4 +547,73 @@ get_ligand_activities_targets_DEgenes_beta = function(receiver_de, receivers_oi, return(list(ligand_activities = ligand_activities, de_genes_df = de_genes_df)) } +#' @title process_geneset_data +#' +#' @description \code{process_geneset_data} Determine ratio's of geneset_oi vs background for a certain logFC/p-val thresholds setting. +#' @usage process_geneset_data(contrast_oi, receiver_de, logFC_threshold = 0.5, p_val_adj = FALSE, p_val_threshold = 0.05) +#' +#' @param contrast_oi Name of one of the contrasts in the celltype_DE / receiver_DE output tibble. +#' @inheritParams get_ligand_activities_targets_DEgenes +#' +#' @return Tibble indicating the nr of up- and down genes per contrast/cell type combination, and an indication whether this is in the recommended ranges +#' +#' @import dplyr +#' @import tidyr +#' +#' @examples +#' \dontrun{ +#' library(dplyr) +#' lr_network = readRDS(url("https://zenodo.org/record/3260758/files/lr_network.rds")) +#' lr_network = lr_network %>% dplyr::rename(ligand = from, receptor = to) %>% dplyr::distinct(ligand, receptor) +#' ligand_target_matrix = readRDS(url("https://zenodo.org/record/3260758/files/ligand_target_matrix.rds")) +#' sample_id = "tumor" +#' group_id = "pEMT" +#' celltype_id = "celltype" +#' batches = NA +#' contrasts_oi = c("'High-Low','Low-High'") +#' contrast_tbl = tibble(contrast = c("High-Low","Low-High"), group = c("High","Low")) +#' receivers_oi = SummarizedExperiment::colData(sce)[,celltype_id] %>% unique() +#' celltype_info = get_avg_frac_exprs_abund(sce = sce, sample_id = sample_id, celltype_id = celltype_id, group_id = group_id) +#' celltype_de = perform_muscat_de_analysis( +#' sce = sce, +#' sample_id = sample_id, +#' celltype_id = celltype_id, +#' group_id = group_id, +#' batches = batches, +#' contrasts = contrasts_oi) +#' receiver_de = celltype_de$de_output_tidy +#' geneset_assessment = contrast_tbl$contrast %>% +#' lapply(process_geneset_data, receiver_de) %>% bind_rows() +#' } +#' +#' @export +#' +process_geneset_data = function(contrast_oi, receiver_de, logFC_threshold = 0.5, p_val_adj = FALSE, p_val_threshold = 0.05){ + + requireNamespace("dplyr") + + celltype_de = receiver_de %>% filter(contrast == contrast_oi) + + background_df = celltype_de %>% group_by(cluster_id) %>% count() %>% ungroup() %>% rename(n_background = n) + if(p_val_adj == FALSE){ + geneset_oi_up_df = celltype_de %>% filter(logFC >= logFC_threshold & p_val <= p_val_threshold) %>% group_by(cluster_id) %>% count() %>% ungroup() %>% rename(n_geneset_up = n) + geneset_oi_down_df = celltype_de %>% filter(logFC <= -logFC_threshold & p_val <= p_val_threshold) %>% group_by(cluster_id) %>% count() %>% ungroup() %>% rename(n_geneset_down = n) + } else { + geneset_oi_up_df = celltype_de %>% filter(logFC >= logFC_threshold & p_adj <= p_val_threshold) %>% group_by(cluster_id) %>% count() %>% ungroup() %>% rename(n_geneset_up = n) + geneset_oi_down_df = celltype_de %>% filter(logFC <= -logFC_threshold & p_adj <= p_val_threshold) %>% group_by(cluster_id) %>% count() %>% ungroup() %>% rename(n_geneset_down = n) + } + + n_df = background_df %>% left_join(geneset_oi_up_df) %>% left_join(geneset_oi_down_df) %>% mutate(n_geneset_up = n_geneset_up %>% tidyr::replace_na(0), n_geneset_down = n_geneset_down %>% tidyr::replace_na(0)) + + geneset_df = n_df %>% mutate(prop_geneset_up = n_geneset_up/n_background, prop_geneset_down = n_geneset_down/n_background) + geneset_df = geneset_df %>% mutate( + in_range_up = prop_geneset_up >= 1/200 & prop_geneset_up <= 1/10, + in_range_down = prop_geneset_down >= 1/200 & prop_geneset_down <= 1/10, + contrast = contrast_oi + ) + + geneset_df = geneset_df %>% mutate(logFC_threshold = logFC_threshold, p_val_threshold = p_val_threshold, adjusted = p_val_adj) + + return(geneset_df) +} diff --git a/R/pipeline.R b/R/pipeline.R index 733572d..18c49d1 100644 --- a/R/pipeline.R +++ b/R/pipeline.R @@ -29,7 +29,7 @@ #' `contrast_tbl = tibble(contrast = c("A-(B+C+D)/3","B-(A+C+D)/3"), group = c("A","B"))` #' @param fraction_cutoff Cutoff indicating the minimum fraction of cells of a cell type in a specific sample that are necessary to consider a gene (e.g. ligand/receptor) as expressed in a sample. #' @param min_sample_prop Parameter to define the minimal required nr of samples in which a gene should be expressed in more than `fraction_cutoff` of cells in that sample (per cell type). This nr of samples is calculated as the `min_sample_prop` fraction of the nr of samples of the smallest group (after considering samples with n_cells >= `min_cells`. Default: `min_sample_prop = 0.50`. Examples: if there are 8 samples in the smallest group, there should be min_sample_prop*8 (= 4 in this example) samples with sufficient fraction of expressing cells. -#' @param scenario Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 2 settings are implemented: "regular" (default) and "lower_DE". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity. +#' @param scenario Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 3 settings are implemented: "regular" (default), "lower_DE", and "no_frac_LR_expr". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity. "no_frac_LR_expr" is the scenario that will exclude the criterion "fraction of samples expressing the LR pair'. This may be beneficial in case of few samples per group. #' @param ligand_activity_down Default: FALSE, downregulatory ligand activity is not considered for prioritization. TRUE: both up- and downregulatory activity are considered for prioritization. #' @param assay_oi_pb Indicates which information of the assay of interest should be used (counts, scaled data,...). Default: "counts". See `muscat::aggregateData`. #' @param fun_oi_pb Indicates way of doing the pseudobulking. Default: "sum". See `muscat::aggregateData`. @@ -157,7 +157,6 @@ multi_nichenet_analysis = function(sce, if(is.double(SummarizedExperiment::colData(sce)[,sample_id])){ stop("SummarizedExperiment::colData(sce)[,sample_id] should be a character vector or a factor") } - # if some of these are factors, and not all levels have syntactically valid names - prompt to change this if(is.factor(SummarizedExperiment::colData(sce)[,celltype_id])){ is_make_names = levels(SummarizedExperiment::colData(sce)[,celltype_id]) == make.names(levels(SummarizedExperiment::colData(sce)[,celltype_id])) @@ -421,7 +420,6 @@ multi_nichenet_analysis = function(sce, DE_info_emp = get_empirical_pvals(DE_info$celltype_de$de_output_tidy) } } - if(empirical_pval == FALSE){ if(findMarkers == TRUE){ celltype_de = DE_info$celltype_de_findmarkers @@ -585,7 +583,6 @@ multi_nichenet_analysis = function(sce, multinichenet_output = make_lite_output_condition_specific(multinichenet_output, top_n_LR = top_n_LR) - } else { print("There are no condition specific cell types in the data. MultiNicheNet analysis is performed in the regular way for all cell types.") multinichenet_output = list( diff --git a/R/prioritization.R b/R/prioritization.R index 0796881..0317a9b 100644 --- a/R/prioritization.R +++ b/R/prioritization.R @@ -105,7 +105,10 @@ generate_prioritization_tables = function(sender_receiver_info, sender_receiver_ } if(scenario == "lower_DE"){ prioritizing_weights = c("de_ligand" = 0.5,"de_receptor" = 0.5,"activity_scaled" = 2,"exprs_ligand" = 1,"exprs_receptor" = 1, "frac_exprs_ligand_receptor" = 1) - } + } + if(scenario == "no_frac_LR_expr"){ + prioritizing_weights = c("de_ligand" = 1,"de_receptor" = 1,"activity_scaled" = 1,"exprs_ligand" = 1,"exprs_receptor" = 1, "frac_exprs_ligand_receptor" = 0) + } # Group prioritization table ------------------------------------------------------------------------------------------------------------------------------------------- diff --git a/man/add_extra_criterion.Rd b/man/add_extra_criterion.Rd index c365f07..1f0336f 100644 --- a/man/add_extra_criterion.Rd +++ b/man/add_extra_criterion.Rd @@ -13,7 +13,7 @@ add_extra_criterion(prioritization_tables, new_criteria_tbl, regular_criteria_tb \item{regular_criteria_tbl}{tibble with 3 columns: criterion, weight, regularization_factor. See example code and vignette for usage.} -\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 2 settings are implemented: "regular" (default) and "lower_DE". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity.} +\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 3 settings are implemented: "regular" (default), "lower_DE", and "no_frac_LR_expr". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity. "no_frac_LR_expr" is the scenario that will exclude the criterion "fraction of samples expressing the LR pair'. This may be beneficial in case of few samples per group.} } \value{ prioritization_tables with updated aggregated prioritization score based on the new criteria (same output as `generate_prioritization_tables`) diff --git a/man/generate_prioritization_tables.Rd b/man/generate_prioritization_tables.Rd index ed35b14..1457e41 100644 --- a/man/generate_prioritization_tables.Rd +++ b/man/generate_prioritization_tables.Rd @@ -21,7 +21,7 @@ Example for `contrasts_oi = c("'A-(B+C+D)/3', 'B-(A+C+D)/3'")`: \item{grouping_tbl}{Data frame showing the groups of each sample (and batches per sample if applicable) (columns: sample and group; and if applicable all batches of interest)} -\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 2 settings are implemented: "regular" (default) and "lower_DE". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity.} +\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 3 settings are implemented: "regular" (default), "lower_DE", and "no_frac_LR_expr". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity. "no_frac_LR_expr" is the scenario that will exclude the criterion "fraction of samples expressing the LR pair'. This may be beneficial in case of few samples per group.} \item{fraction_cutoff}{Cutoff indicating the minimum fraction of cells of a cell type in a specific sample that are necessary to consider a gene (e.g. ligand/receptor) as expressed in a sample.} diff --git a/man/generate_prioritization_tables_condition_specific_celltypes_receiver.Rd b/man/generate_prioritization_tables_condition_specific_celltypes_receiver.Rd index 2efd362..554a618 100644 --- a/man/generate_prioritization_tables_condition_specific_celltypes_receiver.Rd +++ b/man/generate_prioritization_tables_condition_specific_celltypes_receiver.Rd @@ -21,7 +21,7 @@ Example for `contrasts_oi = c("'A-(B+C+D)/3', 'B-(A+C+D)/3'")`: \item{grouping_tbl}{Data frame showing the groups of each sample (and batches per sample if applicable) (columns: sample and group; and if applicable all batches of interest)} -\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 2 settings are implemented: "regular" (default) and "lower_DE". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity.} +\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 3 settings are implemented: "regular" (default), "lower_DE", and "no_frac_LR_expr". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity. "no_frac_LR_expr" is the scenario that will exclude the criterion "fraction of samples expressing the LR pair'. This may be beneficial in case of few samples per group.} \item{fraction_cutoff}{Cutoff indicating the minimum fraction of cells of a cell type in a specific sample that are necessary to consider a gene (e.g. ligand/receptor) as expressed in a sample.} diff --git a/man/generate_prioritization_tables_condition_specific_celltypes_sender.Rd b/man/generate_prioritization_tables_condition_specific_celltypes_sender.Rd index c65c3c6..6e65073 100644 --- a/man/generate_prioritization_tables_condition_specific_celltypes_sender.Rd +++ b/man/generate_prioritization_tables_condition_specific_celltypes_sender.Rd @@ -21,7 +21,7 @@ Example for `contrasts_oi = c("'A-(B+C+D)/3', 'B-(A+C+D)/3'")`: \item{grouping_tbl}{Data frame showing the groups of each sample (and batches per sample if applicable) (columns: sample and group; and if applicable all batches of interest)} -\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 2 settings are implemented: "regular" (default) and "lower_DE". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity.} +\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 3 settings are implemented: "regular" (default), "lower_DE", and "no_frac_LR_expr". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity. "no_frac_LR_expr" is the scenario that will exclude the criterion "fraction of samples expressing the LR pair'. This may be beneficial in case of few samples per group.} \item{fraction_cutoff}{Cutoff indicating the minimum fraction of cells of a cell type in a specific sample that are necessary to consider a gene (e.g. ligand/receptor) as expressed in a sample.} diff --git a/man/multi_nichenet_analysis.Rd b/man/multi_nichenet_analysis.Rd index 62c6ee2..cf0a613 100644 --- a/man/multi_nichenet_analysis.Rd +++ b/man/multi_nichenet_analysis.Rd @@ -47,7 +47,7 @@ Example for `contrasts_oi = c("'A-(B+C+D)/3', 'B-(A+C+D)/3'")`: \item{min_sample_prop}{Parameter to define the minimal required nr of samples in which a gene should be expressed in more than `fraction_cutoff` of cells in that sample (per cell type). This nr of samples is calculated as the `min_sample_prop` fraction of the nr of samples of the smallest group (after considering samples with n_cells >= `min_cells`. Default: `min_sample_prop = 0.50`. Examples: if there are 8 samples in the smallest group, there should be min_sample_prop*8 (= 4 in this example) samples with sufficient fraction of expressing cells.} -\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 2 settings are implemented: "regular" (default) and "lower_DE". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity.} +\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 3 settings are implemented: "regular" (default), "lower_DE", and "no_frac_LR_expr". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity. "no_frac_LR_expr" is the scenario that will exclude the criterion "fraction of samples expressing the LR pair'. This may be beneficial in case of few samples per group.} \item{ligand_activity_down}{Default: FALSE, downregulatory ligand activity is not considered for prioritization. TRUE: both up- and downregulatory activity are considered for prioritization.} diff --git a/man/multi_nichenet_analysis_sampleAgnostic.Rd b/man/multi_nichenet_analysis_sampleAgnostic.Rd index 5d324a7..37cecd3 100644 --- a/man/multi_nichenet_analysis_sampleAgnostic.Rd +++ b/man/multi_nichenet_analysis_sampleAgnostic.Rd @@ -47,7 +47,7 @@ Example for `contrasts_oi = c("'A-(B+C+D)/3', 'B-(A+C+D)/3'")`: \item{min_sample_prop}{Parameter to define the minimal required nr of samples in which a gene should be expressed in more than `fraction_cutoff` of cells in that sample (per cell type). This nr of samples is calculated as the `min_sample_prop` fraction of the nr of samples of the smallest group (after considering samples with n_cells >= `min_cells`. Default: `min_sample_prop = 0.50`. Examples: if there are 8 samples in the smallest group, there should be min_sample_prop*8 (= 4 in this example) samples with sufficient fraction of expressing cells.} -\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 2 settings are implemented: "regular" (default) and "lower_DE". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity.} +\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 3 settings are implemented: "regular" (default), "lower_DE", and "no_frac_LR_expr". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity. "no_frac_LR_expr" is the scenario that will exclude the criterion "fraction of samples expressing the LR pair'. This may be beneficial in case of few samples per group.} \item{ligand_activity_down}{Default: FALSE, downregulatory ligand activity is not considered for prioritization. TRUE: both up- and downregulatory activity are considered for prioritization.} diff --git a/man/prioritize_condition_specific_receiver.Rd b/man/prioritize_condition_specific_receiver.Rd index 7610eec..6cba11c 100644 --- a/man/prioritize_condition_specific_receiver.Rd +++ b/man/prioritize_condition_specific_receiver.Rd @@ -27,7 +27,7 @@ Example for `contrasts_oi = c("'A-(B+C+D)/3', 'B-(A+C+D)/3'")`: \item{ligand_activities_targets_DEgenes}{Output of `get_ligand_activities_targets_DEgenes`} -\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 2 settings are implemented: "regular" (default) and "lower_DE". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity.} +\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 3 settings are implemented: "regular" (default), "lower_DE", and "no_frac_LR_expr". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity. "no_frac_LR_expr" is the scenario that will exclude the criterion "fraction of samples expressing the LR pair'. This may be beneficial in case of few samples per group.} \item{ligand_activity_down}{For prioritization based on ligand activity: consider the max of up- and downregulation (`TRUE`) or consider only upregulated activity (`FALSE`, default from version 2 on).} } diff --git a/man/prioritize_condition_specific_sender.Rd b/man/prioritize_condition_specific_sender.Rd index d2aacc8..fe00b3f 100644 --- a/man/prioritize_condition_specific_sender.Rd +++ b/man/prioritize_condition_specific_sender.Rd @@ -27,7 +27,7 @@ Example for `contrasts_oi = c("'A-(B+C+D)/3', 'B-(A+C+D)/3'")`: \item{ligand_activities_targets_DEgenes}{Output of `get_ligand_activities_targets_DEgenes`} -\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 2 settings are implemented: "regular" (default) and "lower_DE". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity.} +\item{scenario}{Character vector indicating which prioritization weights should be used during the MultiNicheNet analysis. Currently 3 settings are implemented: "regular" (default), "lower_DE", and "no_frac_LR_expr". The setting "regular" is strongly recommended and gives each criterion equal weight. The setting "lower_DE" is recommended in cases your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). It halves the weight for DE criteria, and doubles the weight for ligand activity. "no_frac_LR_expr" is the scenario that will exclude the criterion "fraction of samples expressing the LR pair'. This may be beneficial in case of few samples per group.} \item{ligand_activity_down}{For prioritization based on ligand activity: consider the max of up- and downregulation (`TRUE`) or consider only upregulated activity (`FALSE`, default from version 2 on).} } diff --git a/man/process_geneset_data.Rd b/man/process_geneset_data.Rd new file mode 100644 index 0000000..0f57076 --- /dev/null +++ b/man/process_geneset_data.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ligand_activities.R +\name{process_geneset_data} +\alias{process_geneset_data} +\title{process_geneset_data} +\usage{ +process_geneset_data(contrast_oi, receiver_de, logFC_threshold = 0.5, p_val_adj = FALSE, p_val_threshold = 0.05) +} +\arguments{ +\item{contrast_oi}{Name of one of the contrasts in the celltype_DE / receiver_DE output tibble.} + +\item{receiver_de}{Differential expression analysis output for the receiver cell types. `de_output_tidy` slot of the output of `perform_muscat_de_analysis`.} + +\item{logFC_threshold}{For defining the gene set of interest for NicheNet ligand activity: what is the minimum logFC a gene should have to belong to this gene set? Default: 0.25/} + +\item{p_val_adj}{For defining the gene set of interest for NicheNet ligand activity: should we look at the p-value corrected for multiple testing? Default: FALSE.} + +\item{p_val_threshold}{For defining the gene set of interest for NicheNet ligand activity: what is the maximam p-value a gene should have to belong to this gene set? Default: 0.05.} +} +\value{ +Tibble indicating the nr of up- and down genes per contrast/cell type combination, and an indication whether this is in the recommended ranges +} +\description{ +\code{process_geneset_data} Determine ratio's of geneset_oi vs background for a certain logFC/p-val thresholds setting. +} +\examples{ +\dontrun{ +library(dplyr) +lr_network = readRDS(url("https://zenodo.org/record/3260758/files/lr_network.rds")) +lr_network = lr_network \%>\% dplyr::rename(ligand = from, receptor = to) \%>\% dplyr::distinct(ligand, receptor) +ligand_target_matrix = readRDS(url("https://zenodo.org/record/3260758/files/ligand_target_matrix.rds")) +sample_id = "tumor" +group_id = "pEMT" +celltype_id = "celltype" +batches = NA +contrasts_oi = c("'High-Low','Low-High'") +contrast_tbl = tibble(contrast = c("High-Low","Low-High"), group = c("High","Low")) +receivers_oi = SummarizedExperiment::colData(sce)[,celltype_id] \%>\% unique() +celltype_info = get_avg_frac_exprs_abund(sce = sce, sample_id = sample_id, celltype_id = celltype_id, group_id = group_id) +celltype_de = perform_muscat_de_analysis( + sce = sce, + sample_id = sample_id, + celltype_id = celltype_id, + group_id = group_id, + batches = batches, + contrasts = contrasts_oi) +receiver_de = celltype_de$de_output_tidy +geneset_assessment = contrast_tbl$contrast \%>\% +lapply(process_geneset_data, receiver_de) \%>\% bind_rows() +} + +} diff --git a/vignettes/basic_analysis_steps_MISC.Rmd b/vignettes/basic_analysis_steps_MISC.Rmd index 2f485dc..b1b2ac9 100644 --- a/vignettes/basic_analysis_steps_MISC.Rmd +++ b/vignettes/basic_analysis_steps_MISC.Rmd @@ -1,17 +1,20 @@ --- title: "MultiNicheNet analysis: MIS-C threewise comparison - step-by-step" author: "Robin Browaeys" -date: "2023-06-06" -output: rmarkdown::html_vignette +package: "`r BiocStyle::pkg_ver('multinichenetr')`" +output: + BiocStyle::html_document vignette: > %\VignetteIndexEntry{MultiNicheNet analysis: MIS-C threewise comparison - step-by-step} %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} + %\VignetteEncoding{UTF-8} --- - + ```{r setup, include = FALSE} knitr::opts_chunk$set( @@ -20,40 +23,21 @@ knitr::opts_chunk$set( warning = FALSE, message = FALSE ) +library(BiocStyle) ``` -In this vignette, you can learn how to perform an all-vs-all MultiNicheNet analysis. In this vignette, we start from one SingleCellExperiment object containing cells from both sender and receiver cell types and from different patients. - -A MultiNicheNet analysis can be performed if you have multi-sample, multi-group single-cell data. MultiNicheNet will look for cell-cell communication between the cell types in your data for each sample, and compare the cell-cell communication patterns between the groups of interest. Therefore, the absolute minimum of meta data you need to have, are following columns indicating for each cell: the **group**, **sample** and **cell type**. - -As example expression data of interacting cells, we will here use scRNAseq data of immune cells in MIS-C patients and healthy siblings from this paper of Hoste et al.: [TIM3+ TRBV11-2 T cells and IFNγ signature in patrolling monocytes and CD16+ NK cells delineate MIS-C](https://rupress.org/jem/article/219/2/e20211381/212918/TIM3-TRBV11-2-T-cells-and-IFN-signature-in) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.6362434.svg)](https://doi.org/10.5281/zenodo.6362434) -. MIS-C (multisystem inflammatory syndrome in children) is a novel rare immunodysregulation syndrome that can arise after SARS-CoV-2 infection in children. We will use NicheNet to explore immune cell crosstalk enriched in MIS-C compared to healthy siblings. - -In this vignette, we will prepare the data and analysis parameters, and then perform the MultiNicheNet analysis. - -The different steps of the MultiNicheNet analysis are the following: - -* 0. Preparation of the analysis: load packages, NicheNet LR network & ligand-target matrix, single-cell expression data, and define main settings of the MultiNicheNet analysis - -* 1. Extract cell type abundance and expression information from receiver and sender cell types, and link this expression information for ligands of the sender cell types to the corresponding receptors of the receiver cell types +In this vignette, you can learn how to perform a MultiNicheNet analysis to compare cell-cell communication between conditions of interest. A MultiNicheNet analysis can be performed if you have multi-sample, multi-condition/group single-cell data. We strongly recommend having at least 4 samples in each of the groups/conditions you want to compare. With less samples, the benefits of performing a pseudobulk-based DE analysis are less clear. For those datasets, you can check and run our alternative workflow that makes use of cell-level sample-agnostic differential expression tools. -* 2. Perform genome-wide differential expression analysis of receiver and sender cell types to define DE genes between the conditions of interest. Based on this analysis, we can define the logFC/p-value of ligands in senders and receptors in receivers, and define the set of affected target genes in the receiver. +As input you need a SingleCellExperiment object containing at least the raw count matrix and metadata providing the following information for each cell: the **group**, **sample** and **cell type**. -* 3. Predict NicheNet ligand activities and NicheNet ligand-target links based on these differential expression results +As example expression data of interacting cells, we will here use scRNAseq data of immune cells in MIS-C patients and healthy siblings from this paper of Hoste et al.: [TIM3+ TRBV11-2 T cells and IFNγ signature in patrolling monocytes and CD16+ NK cells delineate MIS-C](https://rupress.org/jem/article/219/2/e20211381/212918/TIM3-TRBV11-2-T-cells-and-IFN-signature-in) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.6362434.svg)](https://doi.org/10.5281/zenodo.6362434). +MIS-C (multisystem inflammatory syndrome in children) is a novel rare immunodysregulation syndrome that can arise after SARS-CoV-2 infection in children. We will use MultiNicheNet to explore immune cell crosstalk enriched in MIS-C compared to healthy siblings and adult COVID-19 patients. -* 4. Use the information collected above to prioritize all sender-ligand---receiver-receptor pairs. +In this vignette, we will first prepare the MultiNicheNet core analysis, then run the several steps in the MultiNicheNet core analysis, and finally interpret the output. -* 5. Calculate correlation in expression between ligand-receptor pairs and their predicted target genes +# Preparation of the MultiNicheNet core analysis -In this vignette, we will demonstrate all these steps one by one. - -After the MultiNicheNet analysis is done, we will explore the output of the analysis with different ways of visualization. - -# Step 0: Preparation of the analysis: load packages, NicheNet LR network & ligand-target matrix, single-cell expression data - -## Step 0.1: Load required packages and NicheNet ligand-receptor network and ligand-target matrix - -```{r} +```{r load-libs, message = FALSE, warning = FALSE} library(SingleCellExperiment) library(dplyr) library(ggplot2) @@ -61,6 +45,10 @@ library(nichenetr) library(multinichenetr) ``` +## Load NicheNet's ligand-receptor network and ligand-target matrix + +MultiNicheNet builds upon the NicheNet framework and uses the same prior knowledge networks (ligand-receptor network and ligand-target matrix, currently v2 version). + The Nichenet v2 networks and matrices for both mouse and human can be downloaded from Zenodo [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.7074291.svg)](https://doi.org/10.5281/zenodo.7074291). We will read these object in for human because our expression data is of human patients. @@ -68,91 +56,208 @@ Gene names are here made syntactically valid via `make.names()` to avoid the los ```{r} organism = "human" +``` + +```{r, results='hide'} +options(timeout = 120) + if(organism == "human"){ - lr_network_all = readRDS(url("https://zenodo.org/record/10229222/files/lr_network_human_allInfo_30112033.rds")) %>% mutate(ligand = convert_alias_to_symbols(ligand, organism = organism), receptor = convert_alias_to_symbols(receptor, organism = organism)) - lr_network_all = lr_network_all %>% mutate(ligand = make.names(ligand), receptor = make.names(receptor)) - lr_network = lr_network_all %>% distinct(ligand, receptor) - ligand_target_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final.rds")) - colnames(ligand_target_matrix) = colnames(ligand_target_matrix) %>% convert_alias_to_symbols(organism = organism) %>% make.names() - rownames(ligand_target_matrix) = rownames(ligand_target_matrix) %>% convert_alias_to_symbols(organism = organism) %>% make.names() + + lr_network_all = + readRDS(url( + "https://zenodo.org/record/10229222/files/lr_network_human_allInfo_30112033.rds" + )) %>% + mutate( + ligand = convert_alias_to_symbols(ligand, organism = organism), + receptor = convert_alias_to_symbols(receptor, organism = organism)) + + lr_network_all = lr_network_all %>% + mutate(ligand = make.names(ligand), receptor = make.names(receptor)) + + lr_network = lr_network_all %>% + distinct(ligand, receptor) + + ligand_target_matrix = readRDS(url( + "https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final.rds" + )) + + colnames(ligand_target_matrix) = colnames(ligand_target_matrix) %>% + convert_alias_to_symbols(organism = organism) %>% make.names() + rownames(ligand_target_matrix) = rownames(ligand_target_matrix) %>% + convert_alias_to_symbols(organism = organism) %>% make.names() + lr_network = lr_network %>% filter(ligand %in% colnames(ligand_target_matrix)) ligand_target_matrix = ligand_target_matrix[, lr_network$ligand %>% unique()] + } else if(organism == "mouse"){ - lr_network_all = readRDS(url("https://zenodo.org/record/10229222/files/lr_network_mouse_allInfo_30112033.rds")) %>% mutate(ligand = convert_alias_to_symbols(ligand, organism = organism), receptor = convert_alias_to_symbols(receptor, organism = organism)) - lr_network = lr_network_all %>% mutate(ligand = make.names(ligand), receptor = make.names(receptor)) %>% distinct(ligand, receptor) - ligand_target_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final_mouse.rds")) - colnames(ligand_target_matrix) = colnames(ligand_target_matrix) %>% convert_alias_to_symbols(organism = organism) %>% make.names() - rownames(ligand_target_matrix) = rownames(ligand_target_matrix) %>% convert_alias_to_symbols(organism = organism) %>% make.names() + + lr_network_all = readRDS(url( + "https://zenodo.org/record/10229222/files/lr_network_mouse_allInfo_30112033.rds" + )) %>% + mutate( + ligand = convert_alias_to_symbols(ligand, organism = organism), + receptor = convert_alias_to_symbols(receptor, organism = organism)) + + lr_network_all = lr_network_all %>% + mutate(ligand = make.names(ligand), receptor = make.names(receptor)) + lr_network = lr_network_all %>% + distinct(ligand, receptor) + + ligand_target_matrix = readRDS(url( + "https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final_mouse.rds" + )) + + colnames(ligand_target_matrix) = colnames(ligand_target_matrix) %>% + convert_alias_to_symbols(organism = organism) %>% make.names() + rownames(ligand_target_matrix) = rownames(ligand_target_matrix) %>% + convert_alias_to_symbols(organism = organism) %>% make.names() + lr_network = lr_network %>% filter(ligand %in% colnames(ligand_target_matrix)) ligand_target_matrix = ligand_target_matrix[, lr_network$ligand %>% unique()] + } ``` -## Step 0.2: Read in SingleCellExperiment Objects +## Read in SingleCellExperiment Object -In this vignette, sender and receiver cell types are in the same SingleCellExperiment object, which we will load here. In this vignette, we will load in a subset of the scRNAseq data of the MIS-C [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.8010790.svg)](https://doi.org/10.5281/zenodo.8010790). For the sake of demonstration, this subset only contains 3 cell types. These celltypes are some of the cell types that were found to be most interesting related to MIS-C according to Hoste et al. +In this vignette, we will load in a subset of the scRNAseq data of the MIS-C [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.8010790.svg)](https://doi.org/10.5281/zenodo.8010790). For the sake of demonstration, this subset only contains 3 cell types. These celltypes are some of the cell types that were found to be most interesting related to MIS-C according to Hoste et al. -If you start from a Seurat object, you can convert it easily to a SingleCellExperiment via `sce = Seurat::as.SingleCellExperiment(seurat_obj, assay = "RNA")`. +If you start from a Seurat object, you can convert it easily to a SingleCellExperiment object via `sce = Seurat::as.SingleCellExperiment(seurat_obj, assay = "RNA")`. Because the NicheNet 2.0. networks are in the most recent version of the official gene symbols, we will make sure that the gene symbols used in the expression data are also updated (= converted from their "aliases" to official gene symbols). Afterwards, we will make them again syntactically valid. -```{r} -sce = readRDS(url("https://zenodo.org/record/8010790/files/sce_subset_misc.rds")) +```{r, results='hide'} +sce = readRDS(url( + "https://zenodo.org/record/8010790/files/sce_subset_misc.rds" + )) sce = alias_to_symbol_SCE(sce, "human") %>% makenames_SCE() ``` -## Step 0.3: Prepare settings of the MultiNicheNet cell-cell communication analysis +## Prepare the settings of the MultiNicheNet cell-cell communication analysis + +In this step, we will formalize our research question into MultiNicheNet input arguments. ### Define in which metadata columns we can find the **group**, **sample** and **cell type** IDs -In this case study, we want to study differences in cell-cell communication patterns between MIS-C patients (M), their healthy siblings (S) and adult patients with severe covid (A). The meta data columns that indicate this disease status is `MIS.C.AgeTier`. +In this case study, we want to study differences in cell-cell communication patterns between MIS-C patients (M), their healthy siblings (S) and adult patients with severe covid (A). The meta data columns that indicate this disease status(=group/condition of interest) is `MIS.C.AgeTier`. Cell type annotations are indicated in the `Annotation_v2.0` column, and the sample is indicated by the `ShortID` column. -If your cells are annotated in multiple hierarchical levels, we recommend using a high level in the hierarchy. This for 2 reasons: 1) MultiNicheNet focuses on differential expression and not differential abundance, and 2) there should be sufficient cells per sample-celltype combination. - -If you would have batch effects or covariates you can correct for, you can define this here as well. - -Important: for categorical covariates and batches, there should be at least one sample for every group-batch combination. If one of your groups/conditions lacks a certain level of your batch, you won't be able to correct for the batch effect because the model is then not able to distinguish batch from group/condition effects. - -Important: the column names of group, sample, cell type, batches and covariates should be syntactically valid (`make.names`) - -Important: All group, sample, cell type, batch and covariate names should be syntactically valid as well (`make.names`) (eg through `SummarizedExperiment::colData(sce)$ShortID = SummarizedExperiment::colData(sce)$ShortID %>% make.names()`) +If your cells are annotated in multiple hierarchical levels, we recommend using a relatively high level in the hierarchy. This for 2 reasons: 1) MultiNicheNet focuses on differential expression and not differential abundance, and 2) there should be sufficient cells per sample-celltype combination (see later). ```{r} sample_id = "ShortID" group_id = "MIS.C.AgeTier" celltype_id = "Annotation_v2.0" +``` + +__Important__: It is required that each sample-id is uniquely assigned to only one condition/group of interest. See the vignettes about paired and multifactorial analysis to see how to define your analysis input when you have multiple samples (and conditions) per patient. + +If you would have batch effects or covariates you can correct for, you can define this here as well. However, this is not applicable to this dataset. Therefore we will use the following NA settings: + +```{r} covariates = NA batches = NA ``` -Sender and receiver cell types also need to be defined. Both are here all cell types in the dataset because we are interested in an All-vs-All analysis. +__Important__: for categorical covariates and batches, there should be at least one sample for every group-batch combination. If one of your groups/conditions lacks a certain level of your batch, you won't be able to correct for the batch effect because the model is then not able to distinguish batch from group/condition effects. + +__Important__: The column names of group, sample, cell type, batches and covariates should be syntactically valid (`make.names`) + +__Important__: All group, sample, cell type, batch and covariate names should be syntactically valid as well (`make.names`) (eg through `SummarizedExperiment::colData(sce)$ShortID = SummarizedExperiment::colData(sce)$ShortID %>% make.names()`) + +### Define the contrasts of interest. + +Here, we want to compare each patient group to the other groups, so the MIS-C (M) group vs healthy control siblings (S) and adult COVID19 patients (A) (= M vs S+A) and so on. We want to know which cell-cell communication patterns are specific for the M vs A+S group, the A vs M+S group and the S vs A+M group. + +To perform this comparison, we need to set the following contrasts: + +```{r} +contrasts_oi = c("'M-(S+A)/2','S-(M+A)/2','A-(S+M)/2'") +``` + +__Very Important__ Note the format to indicate the contrasts! This formatting should be adhered to very strictly, and white spaces are not allowed! Check `?get_DE_info` for explanation about how to define this well. The most important points are that: +*each contrast is surrounded by single quotation marks +*contrasts are separated by a comma without any white space +*all contrasts together are surrounded by double quotation marks. + +If you compare against two groups, you should divide by 2 (as demonstrated here), if you compare against three groups, you should divide by 3 and so on. + +For downstream visualizations and linking contrasts to their main condition, we also need to run the following: +This is necessary because we will also calculate cell-type+condition specificity of ligands and receptors. + +```{r} +contrast_tbl = tibble(contrast = + c("M-(S+A)/2","S-(M+A)/2", "A-(S+M)/2"), + group = c("M","S","A")) +``` + +If you want to compare only two groups (eg M vs S), you can use the following: +`contrasts_oi = c("'M-S','S-M'") ` +`contrast_tbl = tibble(contrast = c("M-S","S-M"), group = c("M","S"))` + +Other vignettes will demonstrate how to formalize different types of research questions. + +### Define the sender and receiver cell types of interest. + +If you want to focus the analysis on specific cell types (e.g. because you know which cell types reside in the same microenvironments based on spatial data), you can define this here. If you have sufficient computational resources and no specific idea of cell-type colocalzations, we recommend to consider all cell types as potential senders and receivers. Later on during analysis of the output it is still possible to zoom in on the cell types that interest you most, but your analysis is not biased to them. + +Here we will consider all cell types in the data: ```{r} senders_oi = SummarizedExperiment::colData(sce)[,celltype_id] %>% unique() receivers_oi = SummarizedExperiment::colData(sce)[,celltype_id] %>% unique() +sce = sce[, SummarizedExperiment::colData(sce)[,celltype_id] %in% + c(senders_oi, receivers_oi) + ] ``` -If the user wants it, it is possible to use only a subset of senders and receivers. Senders and receivers can be entirely different, but also overlapping, or the same. If you don't use all the cell types in your data, we recommend to continue with a subset of your data. +In case you would have samples in your data that do not belong to one of the groups/conditions of interest, we recommend removing them and only keeping conditions of interst: ```{r} -sce = sce[, SummarizedExperiment::colData(sce)[,celltype_id] %in% c(senders_oi, receivers_oi)] +conditions_keep = c("M", "S", "A") +sce = sce[, SummarizedExperiment::colData(sce)[,group_id] %in% + conditions_keep + ] ``` +# Running the MultiNicheNet core analysis + +Now we will run the core of a MultiNicheNet analysis. This analysis consists of the following steps: + +* 1. Cell-type filtering: determine which cell types are sufficiently present +* 2. Gene filtering: determine which genes are sufficiently expressed in each present cell type +* 3. Pseudobulk expression calculation: determine and normalize per-sample pseudobulk expression levels for each expressed gene in each present cell type +* 4. Differential expression (DE) analysis: determine which genes are differentially expressed +* 5. Ligand activity prediction: use the DE analysis output to predict the activity of ligands in receiver cell types and infer their potential target genes +* 6. Prioritization: rank cell-cell communication patterns through multi-criteria prioritization + +Following these steps, one can optionally +* 7. Calculate the across-samples expression correlation between ligand-receptor pairs and target genes +* 8. Prioritize communication patterns involving condition-specific cell types through an alternative prioritization scheme -Now we will go to the first real step of the MultiNicheNet analysis +After these steps, the output can be further explored as we will demonstrate in the "Downstream analysis of the MultiNicheNet output" section. -# Step 1: Extract cell type abundance and expression information from receiver and sender cell types, and link this expression information for ligands of the sender cell types to the corresponding receptors of the receiver cell types +In this vignette, we will demonstrate these steps one-by-one, which offers the most flexibility to the user to assess intermediary results. Other vignettes will demonstrate the use of the `multi_nichenet_analysis` wrapper function. -Since MultiNicheNet will infer group differences at the sample level for each cell type (currently via Muscat - pseudobulking + EdgeR), we need to have sufficient cells per sample of a cell type, and this for both groups. In the following analysis we will set this minimum number of cells per cell type per sample at 10 (recommended minimum). +## Cell-type filtering: determine which cell types are sufficiently present +In this step we will calculate and visualize cell type abundances. This will give an indication about which cell types will be retained in the analysis, and which cell types will be filtered out. + +Since MultiNicheNet will infer group differences at the sample level for each cell type (currently via Muscat - pseudobulking + EdgeR), we need to have sufficient cells per sample of a cell type, and this for all groups. In the following analysis we will set this minimum number of cells per cell type per sample at 10. Samples that have less than `min_cells` cells will be excluded from the analysis for that specific cell type. + ```{r} min_cells = 10 ``` -Now we will calculate abundance and expression information for each cell type / sample / group combination with the following functions. In the output of this function, you can also find some 'Cell type abundance diagnostic plots' that will the users which celltype-sample combinations will be left out later on for DE calculation because the nr of cells is lower than de defined minimum defined here above. If too many celltype-sample combinations don't pass this threshold, we recommend to define your cell types in a more general way (use one level higher of the cell type ontology hierarchy) (eg TH17 CD4T cells --> CD4T cells). +We recommend using `min_cells = 10`, except for datasets with several lowly abundant cell types of interest. For those datasets, we recommend using `min_cells = 5`. ```{r} -abundance_info = get_abundance_info(sce = sce, sample_id = sample_id, group_id = group_id, celltype_id = celltype_id, min_cells = min_cells, senders_oi = senders_oi, receivers_oi = receivers_oi, batches = batches) +abundance_info = get_abundance_info( + sce = sce, + sample_id = sample_id, group_id = group_id, celltype_id = celltype_id, + min_cells = min_cells, + senders_oi = senders_oi, receivers_oi = receivers_oi, + batches = batches + ) ``` First, we will check the cell type abundance diagnostic plots. @@ -161,70 +266,215 @@ First, we will check the cell type abundance diagnostic plots. The first plot visualizes the number of cells per celltype-sample combination, and indicates which combinations are removed during the DE analysis because there are less than `min_cells` in the celltype-sample combination. -```{r, fig.width=15, fig.height=12} +```{r} abundance_info$abund_plot_sample ``` -The red dotted line indicates the required minimum of cells as defined above in `min_cells`. We can see here that some sample-celltype combinations are left out. For the DE analysis in the next step, only cell types will be considered if there are at least two samples per group with a sufficient number of cells. -__Important__: Based on the cell type abundance diagnostics, we recommend users to change their analysis settings if required (eg changing cell type annotation level, batches, ...), before proceeding with the rest of the analysis. +The red dotted line indicates the required minimum of cells as defined above in `min_cells`. We can see here that some sample-celltype combinations are left out. For the DE analysis in the next step, only cell types will be considered if there are at least two samples per group with a sufficient number of cells. But as we can see here: all cell types will be considered for the analysis and there are no condition-specific cell types. -as we can see here: no condition-specific cell types. In case you would have them... +__Important__: Based on the cell type abundance diagnostics, we recommend users to change their analysis settings if required (eg changing cell type annotation level, batches, ...), before proceeding with the rest of the analysis. If too many celltype-sample combinations don't pass this threshold, we recommend to define your cell types in a more general way (use one level higher of the cell type ontology hierarchy) (eg TH17 CD4T cells --> CD4T cells) or use `min_cells = 5` if this would not be possible. -# Step 2: Perform genome-wide differential expression analysis of receiver and sender cell types to define DE genes between the conditions of interest. Based on this analysis, we can define the logFC/p-value of ligands in senders and receptors in receivers, and define the set of affected target genes in the receiver. +### Cell type filtering based on cell type abundance information + +Running the following block of code can help you determine which cell types are condition-specific and which cell types are absent. + +```{r} +sample_group_celltype_df = abundance_info$abundance_data %>% + filter(n > min_cells) %>% + ungroup() %>% + distinct(sample_id, group_id) %>% + cross_join( + abundance_info$abundance_data %>% + ungroup() %>% + distinct(celltype_id) + ) %>% + arrange(sample_id) -Now we will go over to the multi-group, multi-sample differential expression (DE) analysis (also called 'differential state' analysis by the developers of Muscat). +abundance_df = sample_group_celltype_df %>% left_join( + abundance_info$abundance_data %>% ungroup() + ) -### Define the contrasts and covariates of interest for the DE analysis. +abundance_df$n[is.na(abundance_df$n)] = 0 +abundance_df$keep[is.na(abundance_df$keep)] = FALSE +abundance_df_summarized = abundance_df %>% + mutate(keep = as.logical(keep)) %>% + group_by(group_id, celltype_id) %>% + summarise(samples_present = sum((keep))) -Here, we want to compare each patient group to the other groups, so the MIS-C (M) group vs healthy control siblins (S) and adult COVID19 patients (A) etcetera. -To do this comparison, we need to set the following contrasts: +celltypes_absent_one_condition = abundance_df_summarized %>% + filter(samples_present == 0) %>% pull(celltype_id) %>% unique() +# find truly condition-specific cell types by searching for cell types +# truely absent in at least one condition + +celltypes_present_one_condition = abundance_df_summarized %>% + filter(samples_present >= 2) %>% pull(celltype_id) %>% unique() +# require presence in at least 2 samples of one group so +# it is really present in at least one condition + +condition_specific_celltypes = intersect( + celltypes_absent_one_condition, + celltypes_present_one_condition) + +total_nr_conditions = SummarizedExperiment::colData(sce)[,group_id] %>% + unique() %>% length() + +absent_celltypes = abundance_df_summarized %>% + filter(samples_present < 2) %>% + group_by(celltype_id) %>% + count() %>% + filter(n == total_nr_conditions) %>% + pull(celltype_id) + +print("condition-specific celltypes:") +print(condition_specific_celltypes) + +print("absent celltypes:") +print(absent_celltypes) +``` +Absent cell types will be filtered out, condition-specific cell types can be filtered out if you as a user do not want to run the alternative workflow for condition-specific cell types in the optional step 8 of the core MultiNicheNet analysis. ```{r} -contrasts_oi = c("'M-(S+A)/2','S-(M+A)/2','A-(S+M)/2'") +analyse_condition_specific_celltypes = FALSE +``` + +```{r} +if(analyse_condition_specific_celltypes == TRUE){ + senders_oi = senders_oi %>% setdiff(absent_celltypes) + receivers_oi = receivers_oi %>% setdiff(absent_celltypes) +} else { + senders_oi = senders_oi %>% + setdiff(union(absent_celltypes, condition_specific_celltypes)) + receivers_oi = receivers_oi %>% + setdiff(union(absent_celltypes, condition_specific_celltypes)) +} + +sce = sce[, SummarizedExperiment::colData(sce)[,celltype_id] %in% + c(senders_oi, receivers_oi) + ] ``` -__Very Important__ Note the format to indicate the contrasts! This formatting should be adhered to very strictly, and white spaces are not allowed! Check `?get_DE_info` for explanation about how to define this well. The most important things are that: each contrast is surrounded by single quotation marks, contrasts are separated by a comma without any whitespace, and alle contrasts together are surrounded by double quotation marks. If you compare against two groups, you should divide by 2, if you compare against three groups, you should divide by 3 etcetera. +## Gene filtering: determine which genes are sufficiently expressed in each present cell type -For downstream visualizations and linking contrasts to their main group, you need to run the following: +Before running the DE analysis, we will determine which genes are not sufficiently expressed and should be filtered out. +We will perform gene filtering based on a similar procedure as used in `edgeR::filterByExpr`. However, we adapted this procedure to be more interpretable for single-cell datasets. + +For each cell type, we will consider genes expressed if they are expressed in at least a `min_sample_prop` fraction of samples in the condition with the lowest number of samples. By default, we set `min_sample_prop = 0.50`, which means that genes should be expressed in at least 2 samples if the group with lowest nr. of samples has 4 samples like this dataset. ```{r} -contrast_tbl = tibble(contrast = - c("M-(S+A)/2","S-(M+A)/2", "A-(S+M)/2"), - group = c("M","S","A")) +min_sample_prop = 0.50 ``` -If you want to compare only two groups (eg M vs S), you can do like this: -`contrasts_oi = c("'M-S','S-M'") ` -`contrast_tbl = tibble(contrast = c("M-S","S-M"), group = c("M","S"))` - -### Perform the DE analysis for each cell type. +But how do we define which genes are expressed in a sample? For this we will consider genes as expressed if they have non-zero expression values in a `fraction_cutoff` fraction of cells of that cell type in that sample. By default, we set `fraction_cutoff = 0.05`, which means that genes should show non-zero expression values in at least 5% of cells in a sample. -First define expressed genes ```{r} fraction_cutoff = 0.05 -min_sample_prop = 0.50 -frq_list = get_frac_exprs(sce = sce, sample_id = sample_id, celltype_id = celltype_id, group_id = group_id, batches = batches, min_cells = min_cells, fraction_cutoff = fraction_cutoff, min_sample_prop = min_sample_prop) ``` +We recommend using these default values unless there is specific interest in prioritizing (very) weakly expressed interactions. In that case, you could lower the value of `fraction_cutoff`. We explicitly recommend against using `fraction_cutoff > 0.10`. + +Now we will calculate the information required for gene filtering with the following command: + +```{r} +frq_list = get_frac_exprs( + sce = sce, + sample_id = sample_id, celltype_id = celltype_id, group_id = group_id, + batches = batches, + min_cells = min_cells, + fraction_cutoff = fraction_cutoff, min_sample_prop = min_sample_prop) +``` +Now only keep genes that are expressed by at least one cell type: + +```{r} +genes_oi = frq_list$expressed_df %>% + filter(expressed == TRUE) %>% pull(gene) %>% unique() +sce = sce[genes_oi, ] +``` + +## Pseudobulk expression calculation: determine and normalize per-sample pseudobulk expression levels for each expressed gene in each present cell type + +After filtering out absent cell types and genes, we will continue the analysis by calculating the different prioritization criteria that we will use to prioritize cell-cell communication patterns. + +First, we will determine and normalize per-sample pseudobulk expression levels for each expressed gene in each present cell type. The function `process_abundance_expression_info` will link this expression information for ligands of the sender cell types to the corresponding receptors of the receiver cell types. This will later on allow us to define the cell-type specicificy criteria for ligands and receptors. + +```{r} +abundance_expression_info = process_abundance_expression_info( + sce = sce, + sample_id = sample_id, group_id = group_id, celltype_id = celltype_id, + min_cells = min_cells, + senders_oi = senders_oi, receivers_oi = receivers_oi, + lr_network = lr_network, + batches = batches, + frq_list = frq_list, + abundance_info = abundance_info) +``` + +Normalized pseudobulk expression values per gene/celltype/sample can be inspected by: + ```{r} -DE_info = get_DE_info(sce = sce, sample_id = sample_id, group_id = group_id, celltype_id = celltype_id, batches = batches, covariates = covariates, contrasts_oi = contrasts_oi, min_cells = min_cells, expressed_df = frq_list$expressed_df) +abundance_expression_info$celltype_info$pb_df %>% head() +``` + +An average of these sample-level expression values per condition/group can be inspected by: + +```{r} +abundance_expression_info$celltype_info$pb_df_group %>% head() +``` + +Inspecting these values for ligand-receptor interactions can be done by: + +```{r} +abundance_expression_info$sender_receiver_info$pb_df %>% head() +abundance_expression_info$sender_receiver_info$pb_df_group %>% head() +``` + +## Differential expression (DE) analysis: determine which genes are differentially expressed + +In this step, we will perform genome-wide differential expression analysis of receiver and sender cell types to define DE genes between the conditions of interest (as formalized by the `contrasts_oi`). Based on this analysis, we later can define the levels of differential expression of ligands in senders and receptors in receivers, and define the set of affected target genes in the receiver cell types (which will be used for the ligand activity analysis). + +We will apply pseudobulking followed by EdgeR to perform multi-condition multi-sample differential expression (DE) analysis (also called 'differential state' analysis by the developers of Muscat). + +```{r} +DE_info = get_DE_info( + sce = sce, + sample_id = sample_id, group_id = group_id, celltype_id = celltype_id, + batches = batches, covariates = covariates, + contrasts_oi = contrasts_oi, + min_cells = min_cells, + expressed_df = frq_list$expressed_df) ``` ### Check DE results -Table with logFC and p-values for each gene-celltype-contrast: +Check DE output information in table with logFC and p-values for each gene-celltype-contrast: + +```{r} +DE_info$celltype_de$de_output_tidy %>% head() +``` +Evaluate the distributions of p-values: + +```{r} +DE_info$hist_pvals +``` + +These distributions look fine (uniform distribution, except peak at p-value <= 0.05), so we will continue using these regular p-values. In case these p-value distributions look irregular, you can estimate empirical p-values as we will demonstrate in another vignette. ```{r} -DE_info$celltype_de$de_output_tidy %>% arrange(p_adj) %>% head() +empirical_pval = FALSE ``` ```{r} -celltype_de = DE_info$celltype_de$de_output_tidy +if(empirical_pval == TRUE){ + DE_info_emp = get_empirical_pvals(DE_info$celltype_de$de_output_tidy) + celltype_de = DE_info_emp$de_output_tidy_emp %>% select(-p_val, -p_adj) %>% + rename(p_val = p_emp, p_adj = p_adj_emp) +} else { + celltype_de = DE_info$celltype_de$de_output_tidy +} ``` -In the next step, we will combine the DE information of senders and receivers by linking their ligands and receptors together: +### Combine DE information for ligand-senders and receptors-receivers -### Combine DE information for ligand-senders and receptors-receivers (similar to step1 - `abundance_expression_info$sender_receiver_info`) +To end this step, we will combine the DE information of senders and receivers by linking their ligands and receptors together based on the prior knowledge ligand-receptor network. ```{r} sender_receiver_de = combine_sender_receiver_de( @@ -240,114 +490,128 @@ sender_receiver_de = combine_sender_receiver_de( sender_receiver_de %>% head(20) ``` -# Calculate (pseudobulk) expression averages +## Ligand activity prediction: use the DE analysis output to predict the activity of ligands in receiver cell types and infer their potential target genes -```{r} -abundance_expression_info = process_abundance_expression_info(sce = sce, sample_id = sample_id, group_id = group_id, celltype_id = celltype_id, min_cells = min_cells, senders_oi = senders_oi, receivers_oi = receivers_oi, lr_network = lr_network, batches = batches, frq_list = frq_list, abundance_info = abundance_info) -``` +In this step, we will predict NicheNet ligand activities and NicheNet ligand-target links based on these differential expression results. We do this to prioritize interactions based on their predicted effect on a receiver cell type. We will assume that the most important group-specific interactions are those that lead to group-specific gene expression changes in a receiver cell type. + +Similarly to base NicheNet (https://github.com/saeyslab/nichenetr), we use the DE output to create a "geneset of interest": here we assume that DE genes within a cell type may be DE because of differential cell-cell communication processes. In the ligand activity prediction, we will assess the enrichment of target genes of ligands within this geneset of interest. In case high-probabiliy target genes of a ligand are enriched in this set compared to the background of expressed genes, we predict that this ligand may have a high activity. -# Step 3: Predict NicheNet ligand activities and NicheNet ligand-target links based on these differential expression results +Because the ligand activity analysis is an enrichment procedure, it is important that this geneset of interest should contain a sufficient but not too large number of genes. The ratio geneset_oi/background should ideally be between 1/200 and 1/10 (or close to these ratios). -## Define the parameters for the NicheNet ligand activity analysis +To determine the genesets of interest based on DE output, we need to define some logFC and/or p-value thresholds per cell type/contrast combination. In general, we recommend inspecting the nr. of DE genes for all cell types based on the default thresholds and adapting accordingly. By default, we will apply the p-value cutoff on the normal p-values, and not on the p-values corrected for multiple testing. This choice was made because most multi-sample single-cell transcriptomics datasets have just a few samples per group and we might have a lack of statistical power due to pseudobulking. But, if the smallest group >= 20 samples, we typically recommend using p_val_adj = TRUE. When the biological difference between the conditions is very large, we typically recommend increasing the logFC_threshold and/or using p_val_adj = TRUE. -Here, we need to define the thresholds that will be used to consider genes as differentially expressed or not (logFC, p-value, decision whether to use adjusted or normal p-value, minimum fraction of cells that should express a gene in at least one sample in a group, whether to use the normal p-values or empirical p-values). +### Assess geneset_oi-vs-background ratios for different DE output tresholds prior to the NicheNet ligand activity analysis -NicheNet ligand activity will then be calculated as the enrichment of predicted target genes of ligands in this set of DE genes compared to the genomic background. Here we choose for a minimum logFC of 0.50, maximum p-value of 0.05, and minimum fraction of expression of 0.05. +We will first inspect the geneset_oi-vs-background ratios for the default tresholds: ```{r} logFC_threshold = 0.50 p_val_threshold = 0.05 ``` -We will here choose for applying the p-value cutoff on the normal p-values, and not on the p-values corrected for multiple testing. This choice was made here because this dataset has only a few samples per group and we might have a lack of statistical power due to pseudobulking. In case of more samples per group, and a sufficient high number of DE genes per group-celltype (> 50), we would recommend using the adjusted p-values. - ```{r} -# p_val_adj = TRUE p_val_adj = FALSE ``` -For the NicheNet ligand-target inference, we also need to select which top n of the predicted target genes will be considered (here: top 250 targets per ligand). +```{r} +geneset_assessment = contrast_tbl$contrast %>% + lapply( + process_geneset_data, + celltype_de, logFC_threshold, p_val_adj, p_val_threshold + ) %>% + bind_rows() +geneset_assessment +``` +We can see here that for all cell type / contrast combinations, all geneset/background ratio's are within the recommended range (`in_range_up` and `in_range_down` columns). When these geneset/background ratio's would not be within the recommended ranges, we should interpret ligand activity results for these cell types with more caution, or use different thresholds (for these or all cell types). + +For the sake of demonstration, we will also calculate these ratio's in case we would use the adjusted p-value as threshold. ```{r} -top_n_target = 250 +geneset_assessment_adjustedPval = contrast_tbl$contrast %>% + lapply( + process_geneset_data, + celltype_de, logFC_threshold, p_val_adj = TRUE, p_val_threshold + ) %>% + bind_rows() +geneset_assessment_adjustedPval ``` +We can see here that for most cell type / contrast combinations, the geneset/background ratio's are not within the recommended range. Therefore, we will proceed with the default tresholds for the ligand activity analysis + +### Perform the ligand activity analysis and ligand-target inference -The NicheNet ligand activity analysis can be run in parallel for each receiver cell type, by changing the number of cores as defined here. This is only recommended if you have many receiver cell type. +After the ligand activity prediction, we will also infer the predicted target genes of these ligands in each contrast. For this ligand-target inference procedure, we also need to select which top n of the predicted target genes will be considered (here: top 250 targets per ligand). This parameter will not affect the ligand activity predictions. It will only affect ligand-target visualizations and construction of the intercellular regulatory network during the downstream analysis. We recommend users to test other settings in case they would be interested in exploring fewer, but more confident target genes, or vice versa. ```{r} -verbose = TRUE -cores_system = 8 -n.cores = min(cores_system, sender_receiver_de$receiver %>% unique() %>% length()) # use one core per receiver cell type +top_n_target = 250 ``` -## Run the NicheNet ligand activity analysis +The NicheNet ligand activity analysis can be run in parallel for each receiver cell type, by changing the number of cores as defined here. Using more cores will speed up the analysis at the cost of needing more memory. This is only recommended if you have many receiver cell types of interest. -(this might take some time) ```{r} -ligand_activities_targets_DEgenes = suppressMessages(suppressWarnings(get_ligand_activities_targets_DEgenes( - receiver_de = celltype_de, - receivers_oi = receivers_oi, - ligand_target_matrix = ligand_target_matrix, - logFC_threshold = logFC_threshold, - p_val_threshold = p_val_threshold, - p_val_adj = p_val_adj, - top_n_target = top_n_target, - verbose = verbose, - n.cores = n.cores -))) +verbose = TRUE +cores_system = 8 +n.cores = min(cores_system, celltype_de$cluster_id %>% unique() %>% length()) ``` -Check the DE genes used for the activity analysis +Running the ligand activity prediction will take some time (the more cell types and contrasts, the more time) ```{r} -ligand_activities_targets_DEgenes$de_genes_df %>% head(20) +ligand_activities_targets_DEgenes = suppressMessages(suppressWarnings( + get_ligand_activities_targets_DEgenes( + receiver_de = celltype_de, + receivers_oi = intersect(receivers_oi, celltype_de$cluster_id %>% unique()), + ligand_target_matrix = ligand_target_matrix, + logFC_threshold = logFC_threshold, + p_val_threshold = p_val_threshold, + p_val_adj = p_val_adj, + top_n_target = top_n_target, + verbose = verbose, + n.cores = n.cores + ) +)) ``` -Check the output of the activity analysis +You can check the output of the ligand activity and ligand-target inference here: ```{r} ligand_activities_targets_DEgenes$ligand_activities %>% head(20) ``` -# Step 4: Use the information collected above to prioritize all sender-ligand---receiver-receptor pairs. +## Prioritization: rank cell-cell communication patterns through multi-criteria prioritization -In the 3 previous steps, we calculated expression, differential expression and NicheNet activity information. Now we will combine these different types of information in one prioritization scheme. +In the previous steps, we calculated expression, differential expression and NicheNet ligand activity. In the final step, we will now combine all calculated information to rank all sender-ligand---receiver-receptor pairs according to group/condition specificity. We will use the following criteria to prioritize ligand-receptor interactions: -MultiNicheNet allows the user to define the weights of the following criteria to prioritize ligand-receptor interactions: +* Upregulation of the ligand in a sender cell type and/or upregulation of the receptor in a receiver cell type - in the condition of interest. +* Cell-type and condition specific expression of the ligand in the sender cell type and receptor in the receiver cell type (to mitigate the influence of upregulated but still relatively weakly expressed ligands/receptors). +* Sufficiently high expression levels of ligand and receptor in many samples of the same group. +* High NicheNet ligand activity, to further prioritize ligand-receptor pairs based on their predicted effect of the ligand-receptor interaction on the gene expression in the receiver cell type. -* Upregulation of the ligand in a sender cell type and/or upregulation of the receptor in a receiver cell type - in the condition of interest. : `de_ligand` and `de_receptor` -* Sufficiently high expression levels of ligand and receptor in many samples of the same group (to mitigate the influence of outlier samples). : `frac_exprs_ligand_receptor` -* Cell-type and condition specific expression of the ligand in the sender cell type and receptor in the receiver cell type (to mitigate the influence of upregulated but still relatively weakly expressed ligands/receptors) : `exprs_ligand` and `exprs_receptor` -* High NicheNet ligand activity, to further prioritize ligand-receptor pairs based on their predicted effect of the ligand-receptor interaction on the gene expression in the receiver cell type : `activity_scaled` -* High relative abundance of sender and/or receiver in the condition of interest: `abund_sender` and `abund_receiver` (experimental feature - not recommended to give non-zero weights for default analyses) +We will combine these prioritization criteria in a single aggregated prioritization score. In the default setting, we will weigh each of these criteria equally (`scenario = "regular"`). This setting is strongly recommended. However, we also provide some additional setting to accomodate different biological scenarios. The setting `scenario = "lower_DE"` halves the weight for DE criteria and doubles the weight for ligand activity. This is recommended in case your hypothesis is that the differential CCC patterns in your data are less likely to be driven by DE (eg in cases of differential migration into a niche). The setting `scenario = "no_frac_LR_expr"` ignores the criterion "Sufficiently high expression levels of ligand and receptor in many samples of the same group". This may be interesting for users that have data with a limited number of samples and don’t want to penalize interactions if they are not sufficiently expressed in some samples. -The different properties of the sender-ligand---receiver-receptor pairs can be weighted according to the user's preference and insight in the dataset at hand. +Finally, we still need to make one choice. For NicheNet ligand activity we can choose to prioritize ligands that only induce upregulation of target genes (`ligand_activity_down = FALSE`) or can lead potentially lead to both up- and downregulation (`ligand_activity_down = TRUE`). The benefit of `ligand_activity_down = FALSE` is ease of interpretability: prioritized ligand-receptor pairs will be upregulated in the condition of interest, just like their target genes. `ligand_activity_down = TRUE` can be harder to interpret because target genes of some interactions may be upregulated in the other conditions compared to the condition of interest. This is harder to interpret, but may help to pick up interactions that can also repress gene expression. -## Prepare grouping objects +Here we will choose for setting `ligand_activity_down = FALSE` and focus specifically on upregulating ligands. -Make necessary grouping data frame +```{r} +ligand_activity_down = FALSE +``` ```{r} -sender_receiver_tbl = sender_receiver_de %>% dplyr::distinct(sender, receiver) +sender_receiver_tbl = sender_receiver_de %>% distinct(sender, receiver) metadata_combined = SummarizedExperiment::colData(sce) %>% tibble::as_tibble() if(!is.na(batches)){ - grouping_tbl = metadata_combined[,c(sample_id, group_id, batches)] %>% tibble::as_tibble() %>% dplyr::distinct() + grouping_tbl = metadata_combined[,c(sample_id, group_id, batches)] %>% + tibble::as_tibble() %>% distinct() colnames(grouping_tbl) = c("sample","group",batches) } else { - grouping_tbl = metadata_combined[,c(sample_id, group_id)] %>% tibble::as_tibble() %>% dplyr::distinct() + grouping_tbl = metadata_combined[,c(sample_id, group_id)] %>% + tibble::as_tibble() %>% distinct() colnames(grouping_tbl) = c("sample","group") } -``` - -Crucial note: grouping_tbl: group should be the same as in the contrast_tbl, and as in the expression info tables! Rename accordingly if this would not be the case. If you followed the guidelines of this tutorial closely, there should be no problem. - -## Run the prioritization - -```{r} - prioritization_tables = suppressMessages(generate_prioritization_tables( +prioritization_tables = suppressMessages(generate_prioritization_tables( sender_receiver_info = abundance_expression_info$sender_receiver_info, sender_receiver_de = sender_receiver_de, ligand_activities_targets_DEgenes = ligand_activities_targets_DEgenes, @@ -358,7 +622,7 @@ Crucial note: grouping_tbl: group should be the same as in the contrast_tbl, and fraction_cutoff = fraction_cutoff, abundance_data_receiver = abundance_expression_info$abundance_data_receiver, abundance_data_sender = abundance_expression_info$abundance_data_sender, - ligand_activity_down = FALSE # use only upregulatory ligand activities to prioritize -- if you want to consider max of up and down: set this to TRUE + ligand_activity_down = ligand_activity_down )) ``` @@ -369,20 +633,35 @@ First: group-based summary table ```{r} prioritization_tables$group_prioritization_tbl %>% head(20) ``` +This table gives the final prioritization score of each interaction, and the values of the individual prioritization criteria. + +With this step, all required steps are finished. Now, we can optionally still run the following steps +* Calculate the across-samples expression correlation between ligand-receptor pairs and target genes +* Prioritize communication patterns involving condition-specific cell types through an alternative prioritization scheme +Here we will only focus on the expression correlation step: -# Step 5: Add information on prior knowledge and expression correlation between LR and target expression. +## Calculate the across-samples expression correlation between ligand-receptor pairs and target genes In multi-sample datasets, we have the opportunity to look whether expression of ligand-receptor across all samples is correlated with the expression of their by NicheNet predicted target genes. This is what we will do with the following line of code: ```{r} -lr_target_prior_cor = lr_target_prior_cor_inference(prioritization_tables$group_prioritization_tbl$receiver %>% unique(), abundance_expression_info, celltype_de, grouping_tbl, prioritization_tables, ligand_target_matrix, logFC_threshold = logFC_threshold, p_val_threshold = p_val_threshold, p_val_adj = p_val_adj) +lr_target_prior_cor = lr_target_prior_cor_inference( + receivers_oi = prioritization_tables$group_prioritization_tbl$receiver %>% unique(), + abundance_expression_info = abundance_expression_info, + celltype_de = celltype_de, + grouping_tbl = grouping_tbl, + prioritization_tables = prioritization_tables, + ligand_target_matrix = ligand_target_matrix, + logFC_threshold = logFC_threshold, + p_val_threshold = p_val_threshold, + p_val_adj = p_val_adj + ) ``` -# Save all the output of MultiNicheNet +## Save all the output of MultiNicheNet -To avoid needing to redo the analysis later. -All the output written down here is sufficient to make all in-built downstream visualizations. +To avoid needing to redo the analysis later, we will here to save an output object that contains all information to perform all downstream analyses. ```{r} path = "./" @@ -406,24 +685,30 @@ if(save == TRUE){ } ``` +# Interpreting the MultiNicheNet analysis output -# Visualization of the results of the cell-cell communication analysis +## Visualization of differential cell-cell interactions -In a first instance, we will look at the broad overview of prioritized interactions via condition-specific Circos plots. +### Summarizing ChordDiagram circos plots -## Circos plot of top-prioritized links +In a first instance, we will look at the broad overview of prioritized interactions via condition-specific Chordiagram circos plots. We will look here at the top 50 predictions across all contrasts, senders, and receivers of interest. ```{r} -prioritized_tbl_oi_all = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 50, rank_per_group = FALSE) +prioritized_tbl_oi_all = get_top_n_lr_pairs( + multinichenet_output$prioritization_tables, + top_n = 50, + rank_per_group = FALSE + ) ``` - -```{r, fig.width=15, fig.height=12} -prioritized_tbl_oi = multinichenet_output$prioritization_tables$group_prioritization_tbl %>% +```{r} +prioritized_tbl_oi = + multinichenet_output$prioritization_tables$group_prioritization_tbl %>% filter(id %in% prioritized_tbl_oi_all$id) %>% - distinct(id, sender, receiver, ligand, receptor, group) %>% left_join(prioritized_tbl_oi_all) + distinct(id, sender, receiver, ligand, receptor, group) %>% + left_join(prioritized_tbl_oi_all) prioritized_tbl_oi$prioritization_score[is.na(prioritized_tbl_oi$prioritization_score)] = 0 senders_receivers = union(prioritized_tbl_oi$sender %>% unique(), prioritized_tbl_oi$receiver %>% unique()) %>% sort() @@ -434,9 +719,11 @@ colors_receiver = RColorBrewer::brewer.pal(n = length(senders_receivers), name = circos_list = make_circos_group_comparison(prioritized_tbl_oi, colors_sender, colors_receiver) ``` -## Visualization of scaled ligand-receptor pseudobulk products and ligand activity +### Interpretable bubble plots + +Whereas these ChordDiagrams show the most specific interactions per group, they don't give insights into the data behind these predictions. Therefore we will now look at visualizations that indicate the different prioritization criteria used in MultiNicheNet. -Now we will visualize per sample the scaled product of ligand and receptor expression. Samples that were left out of the DE analysis are indicated with a smaller dot (this helps to indicate the samples that did not contribute to the calculation of the logFC, and thus not contributed to the final prioritization) +In the next type of plots, we will 1) visualize the per-sample scaled product of normalized ligand and receptor pseudobulk expression, 2) visualize the scaled ligand activities, 3) cell-type specificity. We will now check the top 50 interactions specific for the MIS-C group @@ -444,47 +731,71 @@ We will now check the top 50 interactions specific for the MIS-C group group_oi = "M" ``` -```{r, fig.height=13, fig.width=15} -prioritized_tbl_oi_M_50 = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 50, groups_oi = group_oi) +```{r} +prioritized_tbl_oi_M_50 = get_top_n_lr_pairs( + multinichenet_output$prioritization_tables, + top_n = 50, + groups_oi = group_oi) ``` -```{r, fig.height=13, fig.width=18} -plot_oi = make_sample_lr_prod_activity_plots(multinichenet_output$prioritization_tables, prioritized_tbl_oi_M_50) +```{r, fig.height=13, fig.width=16} +plot_oi = make_sample_lr_prod_activity_plots( + multinichenet_output$prioritization_tables, + prioritized_tbl_oi_M_50) plot_oi ``` +Samples that were left out of the DE analysis are indicated with a smaller dot (this helps to indicate the samples that did not contribute to the calculation of the logFC, and thus not contributed to the final prioritization) As a further help for further prioritization, we can assess the level of curation of these LR pairs as defined by the Intercellular Communication part of the Omnipath database ```{r} -prioritized_tbl_oi_M_50_omnipath = prioritized_tbl_oi_M_50 %>% inner_join(lr_network_all) +prioritized_tbl_oi_M_50_omnipath = prioritized_tbl_oi_M_50 %>% + inner_join(lr_network_all) ``` +Now we add this to the bubble plot visualization: ```{r, fig.height=13, fig.width=16} -plot_oi = make_sample_lr_prod_activity_plots_Omnipath(multinichenet_output$prioritization_tables, prioritized_tbl_oi_M_50_omnipath) +plot_oi = make_sample_lr_prod_activity_plots_Omnipath( + multinichenet_output$prioritization_tables, + prioritized_tbl_oi_M_50_omnipath) plot_oi ``` -As you can see, the HEBP1-FPR2 interaction has no Omnipath DB scores. This is because this LR pair was not documented by the Omnipath LR database. Instead it was documented by the original NicheNet LR network (source: Guide2Pharmacology) as can be seen in the table. - - +As you can see, the HEBP1-FPR2 interaction has no Omnipath DB scores. This is because this LR pair was not documented by the Omnipath LR database. Instead it was documented by the original NicheNet LR network (source: Guide2Pharmacology) as can be seen in the table. Further note: Typically, there are way more than 50 differentially expressed and active ligand-receptor pairs per group across all sender-receiver combinations. Therefore it might be useful to zoom in on specific cell types as senders/receivers: Eg M_Monocyte_CD16 as receiver: -```{r, fig.height=13, fig.width=15} -prioritized_tbl_oi_M_50 = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 50, groups_oi = group_oi, receivers_oi = "M_Monocyte_CD16") +```{r} +prioritized_tbl_oi_M_50 = get_top_n_lr_pairs( + multinichenet_output$prioritization_tables, + 50, + groups_oi = group_oi, + receivers_oi = "M_Monocyte_CD16") +``` -plot_oi = make_sample_lr_prod_activity_plots(multinichenet_output$prioritization_tables, prioritized_tbl_oi_M_50) +```{r, fig.height=13, fig.width=16} +plot_oi = make_sample_lr_prod_activity_plots_Omnipath( + multinichenet_output$prioritization_tables, + prioritized_tbl_oi_M_50 %>% inner_join(lr_network_all)) plot_oi ``` Eg M_Monocyte_CD16 as sender: -```{r, fig.height=13, fig.width=15} -prioritized_tbl_oi_M_50 = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 50, groups_oi = group_oi, senders_oi = "M_Monocyte_CD16") +```{r} +prioritized_tbl_oi_M_50 = get_top_n_lr_pairs( + multinichenet_output$prioritization_tables, + 50, + groups_oi = group_oi, + senders_oi = "M_Monocyte_CD16") +``` -plot_oi = make_sample_lr_prod_activity_plots(multinichenet_output$prioritization_tables, prioritized_tbl_oi_M_50) +```{r, fig.height=13, fig.width=16} +plot_oi = make_sample_lr_prod_activity_plots_Omnipath( + multinichenet_output$prioritization_tables, + prioritized_tbl_oi_M_50 %>% inner_join(lr_network_all)) plot_oi ``` @@ -494,70 +805,125 @@ You can make these plots also for the other groups, like we will illustrate now group_oi = "S" ``` -```{r, fig.height=13, fig.width=15} -prioritized_tbl_oi_S_50 = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 50, groups_oi = group_oi) - -plot_oi = make_sample_lr_prod_activity_plots(multinichenet_output$prioritization_tables, prioritized_tbl_oi_S_50) +```{r, fig.height=13, fig.width=18} +prioritized_tbl_oi_S_50 = get_top_n_lr_pairs( + multinichenet_output$prioritization_tables, + 50, + groups_oi = group_oi) + +plot_oi = make_sample_lr_prod_activity_plots_Omnipath( + multinichenet_output$prioritization_tables, + prioritized_tbl_oi_S_50 %>% inner_join(lr_network_all)) plot_oi ``` __Note__: Use `make_sample_lr_prod_activity_batch_plots` if you have batches and want to visualize them on this plot! -## Visualization of ligand-activity, ligand-target links, and target gene expression +## Visualization of differential ligand-target links + +### Without filtering of target genes based on LR-target expression correlation In another type of plot, we can visualize the ligand activities for a group-receiver combination, and show the predicted ligand-target links, and also the expression of the predicted target genes across samples. -For this, we now need to define a receiver cell type of interest. As example, we will take `L_T_TIM3._CD38._HLADR.` cells as receiver, and look at the top 20 senderLigand-receiverReceptor pairs with these cells as receiver. +For this, we now need to define a receiver cell type of interest. As example, we will take `L_T_TIM3._CD38._HLADR.` cells as receiver, and look at the top 10 senderLigand-receiverReceptor pairs with these cells as receiver. ```{r} group_oi = "M" -receiver_oi = "L_T_TIM3._CD38._HLADR." -prioritized_tbl_oi_M_10 = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 10, groups_oi = group_oi, receivers_oi = receiver_oi) -``` - -```{r, fig.width=20, fig.height=10} -combined_plot = make_ligand_activity_target_plot(group_oi, receiver_oi, prioritized_tbl_oi_M_10, multinichenet_output$prioritization_tables, multinichenet_output$ligand_activities_targets_DEgenes, contrast_tbl, multinichenet_output$grouping_tbl, multinichenet_output$celltype_info, ligand_target_matrix, plot_legend = FALSE) +receiver_oi = "M_Monocyte_CD16" +prioritized_tbl_oi_M_10 = get_top_n_lr_pairs( + multinichenet_output$prioritization_tables, + 10, + groups_oi = group_oi, + receivers_oi = receiver_oi) +``` + +```{r, fig.width=20, fig.height=7} +combined_plot = make_ligand_activity_target_plot( + group_oi, + receiver_oi, + prioritized_tbl_oi_M_10, + multinichenet_output$prioritization_tables, + multinichenet_output$ligand_activities_targets_DEgenes, contrast_tbl, + multinichenet_output$grouping_tbl, + multinichenet_output$celltype_info, + ligand_target_matrix, + plot_legend = FALSE) combined_plot ``` -__Note__ Use `make_DEgene_dotplot_pseudobulk_batch` if you want to indicate the batch of each sample to the plot - What if there is a specific ligand you are interested in? ```{r} group_oi = "M" -receiver_oi = "L_T_TIM3._CD38._HLADR." +receiver_oi = "M_Monocyte_CD16" ligands_oi = c("IFNG","IL15") -prioritized_tbl_ligands_oi = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 10000, groups_oi = group_oi, receivers_oi = receiver_oi) %>% filter(ligand %in% ligands_oi) # ligands should still be in the output tables of course -``` - -```{r, fig.width=20, fig.height=10} -combined_plot = make_ligand_activity_target_plot(group_oi, receiver_oi, prioritized_tbl_ligands_oi, multinichenet_output$prioritization_tables, multinichenet_output$ligand_activities_targets_DEgenes, contrast_tbl, multinichenet_output$grouping_tbl, multinichenet_output$celltype_info, ligand_target_matrix, plot_legend = FALSE) +prioritized_tbl_ligands_oi = get_top_n_lr_pairs( + multinichenet_output$prioritization_tables, + 10000, + groups_oi = group_oi, + receivers_oi = receiver_oi + ) %>% filter(ligand %in% ligands_oi) # ligands should still be in the output tables of course +``` + +```{r, fig.width=20, fig.height=7} +combined_plot = make_ligand_activity_target_plot( + group_oi, + receiver_oi, + prioritized_tbl_ligands_oi, + multinichenet_output$prioritization_tables, + multinichenet_output$ligand_activities_targets_DEgenes, + contrast_tbl, + multinichenet_output$grouping_tbl, + multinichenet_output$celltype_info, + ligand_target_matrix, + plot_legend = FALSE) combined_plot ``` -## Visualization of expression-correlated target genes of ligand-receptor pairs +### With filtering of target genes based on LR-target expression correlation -Before, we had calculated the correlation in expression between ligand-receptor pairs and DE genes. Now we will filter out correlated ligand-receptor --> target links that both show high expression correlation (spearman or activity > 0.50 in this example) and have some prior knowledge to support their link. +In the previous plots, target genes were shown that are predicted as target gene of ligands based on prior knowledge. However, we can use the multi-sample nature of this data to filter target genes based on expression correlation between the upstream ligand-receptor pair and the downstream target gene. We will filter out correlated ligand-receptor --> target links that both show high expression correlation (spearman or pearson correlation > 0.50 in this example) and have some prior knowledge to support their link. Note that you can only make these visualization if you ran step 7 of the core MultiNicheNet analysis. ```{r} group_oi = "M" receiver_oi = "M_Monocyte_CD16" -lr_target_prior_cor_filtered = multinichenet_output$lr_target_prior_cor %>% inner_join(multinichenet_output$ligand_activities_targets_DEgenes$ligand_activities %>% distinct(ligand, target, direction_regulation, contrast)) %>% inner_join(contrast_tbl) %>% filter(group == group_oi, receiver == receiver_oi) -lr_target_prior_cor_filtered_up = lr_target_prior_cor_filtered %>% filter(direction_regulation == "up") %>% filter( (rank_of_target < top_n_target) & (pearson > 0.50 | spearman > 0.50)) -lr_target_prior_cor_filtered_down = lr_target_prior_cor_filtered %>% filter(direction_regulation == "down") %>% filter( (rank_of_target < top_n_target) & (pearson < -0.50 | spearman < -0.50)) # downregulation -- negative correlation -lr_target_prior_cor_filtered = bind_rows(lr_target_prior_cor_filtered_up, lr_target_prior_cor_filtered_down) +lr_target_prior_cor_filtered = multinichenet_output$lr_target_prior_cor %>% + inner_join( + multinichenet_output$ligand_activities_targets_DEgenes$ligand_activities %>% + distinct(ligand, target, direction_regulation, contrast) + ) %>% + inner_join(contrast_tbl) %>% filter(group == group_oi, receiver == receiver_oi) + +lr_target_prior_cor_filtered_up = lr_target_prior_cor_filtered %>% + filter(direction_regulation == "up") %>% + filter( (rank_of_target < top_n_target) & (pearson > 0.50 | spearman > 0.50)) +lr_target_prior_cor_filtered_down = lr_target_prior_cor_filtered %>% + filter(direction_regulation == "down") %>% + filter( (rank_of_target < top_n_target) & (pearson < -0.50 | spearman < -0.50)) # downregulation -- negative correlation +lr_target_prior_cor_filtered = bind_rows( + lr_target_prior_cor_filtered_up, + lr_target_prior_cor_filtered_down) ``` Now we will visualize the top correlated target genes for the LR pairs that are also in the top 50 LR pairs discriminating the groups from each other: ```{r} -prioritized_tbl_oi = get_top_n_lr_pairs(prioritization_tables, 50, groups_oi = group_oi, receivers_oi = receiver_oi) +prioritized_tbl_oi = get_top_n_lr_pairs( + prioritization_tables, + 50, + groups_oi = group_oi, + receivers_oi = receiver_oi) ``` - -```{r, fig.width=26, fig.height=16} -lr_target_correlation_plot = make_lr_target_correlation_plot(multinichenet_output$prioritization_tables, prioritized_tbl_oi, lr_target_prior_cor_filtered , multinichenet_output$grouping_tbl, multinichenet_output$celltype_info, receiver_oi,plot_legend = FALSE) +```{r, fig.width=28, fig.height=16} +lr_target_correlation_plot = make_lr_target_correlation_plot( + multinichenet_output$prioritization_tables, + prioritized_tbl_oi, + lr_target_prior_cor_filtered , + multinichenet_output$grouping_tbl, + multinichenet_output$celltype_info, + receiver_oi, + plot_legend = FALSE) lr_target_correlation_plot$combined_plot ``` @@ -568,27 +934,52 @@ ligand_oi = "IFNG" receptor_oi = "IFNGR2" sender_oi = "L_T_TIM3._CD38._HLADR." receiver_oi = "M_Monocyte_CD16" -lr_target_scatter_plot = make_lr_target_scatter_plot(multinichenet_output$prioritization_tables, ligand_oi, receptor_oi, sender_oi, receiver_oi, multinichenet_output$celltype_info, multinichenet_output$grouping_tbl, lr_target_prior_cor_filtered) +lr_target_scatter_plot = make_lr_target_scatter_plot( + multinichenet_output$prioritization_tables, + ligand_oi, receptor_oi, sender_oi, receiver_oi, + multinichenet_output$celltype_info, + multinichenet_output$grouping_tbl, + lr_target_prior_cor_filtered) lr_target_scatter_plot ``` -## Intercellular regulatory network systems view +## Intercellular regulatory network inference and visualization In the plots before, we demonstrated that some DE genes have both expression correlation and prior knowledge support to be downstream of ligand-receptor pairs. Interestingly, some target genes can be ligands or receptors themselves. This illustrates that cells can send signals to other cells, who as a response to these signals produce signals themselves to feedback to the original sender cells, or who will effect other cell types. -As last plot, we can generate a 'systems' view of these intercellular feedback and cascade processes than can be occuring between the different cell populations involved. In this plot, we will draw links between ligands of sender cell types their ligand/receptor-annotated target genes in receiver cell types. So links are ligand-target links (= gene regulatory links) and not ligand-receptor protein-protein interactions! - -```{r} -prioritized_tbl_oi = get_top_n_lr_pairs(prioritization_tables, 250, rank_per_group = FALSE) - -lr_target_prior_cor_filtered = prioritization_tables$group_prioritization_tbl$group %>% unique() %>% lapply(function(group_oi){ - lr_target_prior_cor_filtered = multinichenet_output$lr_target_prior_cor %>% inner_join(multinichenet_output$ligand_activities_targets_DEgenes$ligand_activities %>% distinct(ligand, target, direction_regulation, contrast)) %>% inner_join(contrast_tbl) %>% filter(group == group_oi) - lr_target_prior_cor_filtered_up = lr_target_prior_cor_filtered %>% filter(direction_regulation == "up") %>% filter( (rank_of_target < top_n_target) & (pearson > 0.50 | spearman > 0.50)) - lr_target_prior_cor_filtered_down = lr_target_prior_cor_filtered %>% filter(direction_regulation == "down") %>% filter( (rank_of_target < top_n_target) & (pearson < -0.50 | spearman < -0.50)) - lr_target_prior_cor_filtered = bind_rows(lr_target_prior_cor_filtered_up, lr_target_prior_cor_filtered_down) +As last plot, we can generate a 'systems' view of these intercellular feedback and cascade processes than can be occuring between the different cell populations involved. In this plot, we will draw links between ligands of sender cell types their ligand/receptor-annotated target genes in receiver cell types. So links are ligand-target links (= gene regulatory links) and not ligand-receptor protein-protein interactions! We will infer this intercellular regulatory network here for the top250 interactions. + +```{r} +prioritized_tbl_oi = get_top_n_lr_pairs( + prioritization_tables, + 250, + rank_per_group = FALSE) + +lr_target_prior_cor_filtered = + prioritization_tables$group_prioritization_tbl$group %>% unique() %>% + lapply(function(group_oi){ + lr_target_prior_cor_filtered = multinichenet_output$lr_target_prior_cor %>% + inner_join( + multinichenet_output$ligand_activities_targets_DEgenes$ligand_activities %>% + distinct(ligand, target, direction_regulation, contrast) + ) %>% + inner_join(contrast_tbl) %>% filter(group == group_oi) + + lr_target_prior_cor_filtered_up = lr_target_prior_cor_filtered %>% + filter(direction_regulation == "up") %>% + filter( (rank_of_target < top_n_target) & (pearson > 0.50 | spearman > 0.50)) + + lr_target_prior_cor_filtered_down = lr_target_prior_cor_filtered %>% + filter(direction_regulation == "down") %>% + filter( (rank_of_target < top_n_target) & (pearson < -0.50 | spearman < -0.50)) + lr_target_prior_cor_filtered = bind_rows( + lr_target_prior_cor_filtered_up, + lr_target_prior_cor_filtered_down + ) }) %>% bind_rows() -lr_target_df = lr_target_prior_cor_filtered %>% distinct(group, sender, receiver, ligand, receptor, id, target, direction_regulation) +lr_target_df = lr_target_prior_cor_filtered %>% + distinct(group, sender, receiver, ligand, receptor, id, target, direction_regulation) ``` ```{r} @@ -603,55 +994,14 @@ network_graph = visualize_network(network, colors_sender) network_graph$plot ``` -An interesting follow-up question would be to find nodes that are at the top of the hierarchy in this network, and could thus explain the cascade of intercellular communication events we observe. - -Using PageRank on the opposite-directionality graph can help identify important nodes that act as sources in the hierarchical structure. By running PageRank on the reversed graph (where edges point in the opposite direction), you can potentially pinpoint nodes that have significant influence or act as sources of information in your network. - -We will exemplify for the M-group - -```{r,} -library(igraph) - -# Create the normal graph (normal directionality) -normal_edges <- get.data.frame(network_graph$network_igraph) %>% filter(group == "M") %>% .[, c("from", "to")] -normal_graph <- graph_from_data_frame(normal_edges, directed = TRUE) -# Calculate PageRank on the normal graph -pagerank_scores <- page_rank(normal_graph)$vector %>% sort(decreasing = TRUE) -pr_score_tbl = tibble(node = names(pagerank_scores), score_normal = pagerank_scores) %>% arrange(-score_normal) - -# Create the reversed graph (opposite directionality) -reversed_edges <- get.data.frame(network_graph$network_igraph) %>% filter(group == "M") %>% .[, c("to", "from")] # Swap 'from' and 'to' columns -reversed_graph <- graph_from_data_frame(reversed_edges, directed = TRUE) -# Calculate PageRank on the reversed graph -reversed_pagerank_scores <- page_rank(reversed_graph)$vector %>% sort(decreasing = TRUE) - -pr_score_tbl = pr_score_tbl %>% inner_join( - tibble(node = names(reversed_pagerank_scores), score_reverse = reversed_pagerank_scores) -) %>% arrange(-score_reverse) -pr_score_tbl -``` - -Let's zoom now in on the M-group network to check this - -```{r} -network = infer_intercellular_regulatory_network(lr_target_df %>% filter(group == "M"), prioritized_tbl_oi %>% filter(group == "M")) -network$links %>% head() -network$nodes %>% head() -``` - -```{r, fig.width=30, fig.height=20} -colors_sender["L_T_TIM3._CD38._HLADR."] = "pink" # the original yellow background with white font is not very readable -network_graph = visualize_network(network, colors_sender) -network_graph$plot -``` - Interestingly, we can also use this network to further prioritize differential CCC interactions. Here we will assume that the most important LR interactions are the ones that are involved in this intercellular regulatory network. We can get these interactions as follows: ```{r} network$prioritized_lr_interactions ``` ```{r, fig.width=30, fig.height=12} -prioritized_tbl_oi_network = prioritized_tbl_oi %>% inner_join(network$prioritized_lr_interactions) +prioritized_tbl_oi_network = prioritized_tbl_oi %>% inner_join( + network$prioritized_lr_interactions) prioritized_tbl_oi_network ``` @@ -660,18 +1010,12 @@ Visualize now the expression and activity of these interactions for the M group group_oi = "M" ``` -```{r, fig.height=16, fig.width=18} +```{r, fig.height=19, fig.width=18} prioritized_tbl_oi_M = prioritized_tbl_oi_network %>% filter(group == group_oi) -plot_oi = make_sample_lr_prod_activity_plots(multinichenet_output$prioritization_tables, prioritized_tbl_oi_M) +plot_oi = make_sample_lr_prod_activity_plots_Omnipath( + multinichenet_output$prioritization_tables, + prioritized_tbl_oi_M %>% inner_join(lr_network_all) + ) plot_oi ``` - -## Investingating the data sources behind the ligand-receptor interactions - -Now link the prioritized interactions to all the information of the ligand-receptor network. This can help to check which are the specific data sources behind the inferred interactions, which can aid in prioritizing interactions for further validation. - -Check OmnipathR for more information with respect to how to interpret the different metrics specific for the interactions coming from the Omnipath database. -```{r} -prioritized_tbl_oi_M %>% inner_join(lr_network_all) %>% arrange(-n_resources) # sort interactions based on nr of resources that document them -``` diff --git a/vignettes/basic_analysis_steps_MISC_old.Rmd b/vignettes/basic_analysis_steps_MISC_old.Rmd new file mode 100644 index 0000000..2f485dc --- /dev/null +++ b/vignettes/basic_analysis_steps_MISC_old.Rmd @@ -0,0 +1,677 @@ +--- +title: "MultiNicheNet analysis: MIS-C threewise comparison - step-by-step" +author: "Robin Browaeys" +date: "2023-06-06" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{MultiNicheNet analysis: MIS-C threewise comparison - step-by-step} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + + + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + # comment = "#>", + warning = FALSE, + message = FALSE +) +``` + +In this vignette, you can learn how to perform an all-vs-all MultiNicheNet analysis. In this vignette, we start from one SingleCellExperiment object containing cells from both sender and receiver cell types and from different patients. + +A MultiNicheNet analysis can be performed if you have multi-sample, multi-group single-cell data. MultiNicheNet will look for cell-cell communication between the cell types in your data for each sample, and compare the cell-cell communication patterns between the groups of interest. Therefore, the absolute minimum of meta data you need to have, are following columns indicating for each cell: the **group**, **sample** and **cell type**. + +As example expression data of interacting cells, we will here use scRNAseq data of immune cells in MIS-C patients and healthy siblings from this paper of Hoste et al.: [TIM3+ TRBV11-2 T cells and IFNγ signature in patrolling monocytes and CD16+ NK cells delineate MIS-C](https://rupress.org/jem/article/219/2/e20211381/212918/TIM3-TRBV11-2-T-cells-and-IFN-signature-in) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.6362434.svg)](https://doi.org/10.5281/zenodo.6362434) +. MIS-C (multisystem inflammatory syndrome in children) is a novel rare immunodysregulation syndrome that can arise after SARS-CoV-2 infection in children. We will use NicheNet to explore immune cell crosstalk enriched in MIS-C compared to healthy siblings. + +In this vignette, we will prepare the data and analysis parameters, and then perform the MultiNicheNet analysis. + +The different steps of the MultiNicheNet analysis are the following: + +* 0. Preparation of the analysis: load packages, NicheNet LR network & ligand-target matrix, single-cell expression data, and define main settings of the MultiNicheNet analysis + +* 1. Extract cell type abundance and expression information from receiver and sender cell types, and link this expression information for ligands of the sender cell types to the corresponding receptors of the receiver cell types + +* 2. Perform genome-wide differential expression analysis of receiver and sender cell types to define DE genes between the conditions of interest. Based on this analysis, we can define the logFC/p-value of ligands in senders and receptors in receivers, and define the set of affected target genes in the receiver. + +* 3. Predict NicheNet ligand activities and NicheNet ligand-target links based on these differential expression results + +* 4. Use the information collected above to prioritize all sender-ligand---receiver-receptor pairs. + +* 5. Calculate correlation in expression between ligand-receptor pairs and their predicted target genes + +In this vignette, we will demonstrate all these steps one by one. + +After the MultiNicheNet analysis is done, we will explore the output of the analysis with different ways of visualization. + +# Step 0: Preparation of the analysis: load packages, NicheNet LR network & ligand-target matrix, single-cell expression data + +## Step 0.1: Load required packages and NicheNet ligand-receptor network and ligand-target matrix + +```{r} +library(SingleCellExperiment) +library(dplyr) +library(ggplot2) +library(nichenetr) +library(multinichenetr) +``` + +The Nichenet v2 networks and matrices for both mouse and human can be downloaded from Zenodo [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.7074291.svg)](https://doi.org/10.5281/zenodo.7074291). + +We will read these object in for human because our expression data is of human patients. +Gene names are here made syntactically valid via `make.names()` to avoid the loss of genes (eg H2-M3) in downstream visualizations. + +```{r} +organism = "human" +if(organism == "human"){ + lr_network_all = readRDS(url("https://zenodo.org/record/10229222/files/lr_network_human_allInfo_30112033.rds")) %>% mutate(ligand = convert_alias_to_symbols(ligand, organism = organism), receptor = convert_alias_to_symbols(receptor, organism = organism)) + lr_network_all = lr_network_all %>% mutate(ligand = make.names(ligand), receptor = make.names(receptor)) + lr_network = lr_network_all %>% distinct(ligand, receptor) + ligand_target_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final.rds")) + colnames(ligand_target_matrix) = colnames(ligand_target_matrix) %>% convert_alias_to_symbols(organism = organism) %>% make.names() + rownames(ligand_target_matrix) = rownames(ligand_target_matrix) %>% convert_alias_to_symbols(organism = organism) %>% make.names() + lr_network = lr_network %>% filter(ligand %in% colnames(ligand_target_matrix)) + ligand_target_matrix = ligand_target_matrix[, lr_network$ligand %>% unique()] +} else if(organism == "mouse"){ + lr_network_all = readRDS(url("https://zenodo.org/record/10229222/files/lr_network_mouse_allInfo_30112033.rds")) %>% mutate(ligand = convert_alias_to_symbols(ligand, organism = organism), receptor = convert_alias_to_symbols(receptor, organism = organism)) + lr_network = lr_network_all %>% mutate(ligand = make.names(ligand), receptor = make.names(receptor)) %>% distinct(ligand, receptor) + ligand_target_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final_mouse.rds")) + colnames(ligand_target_matrix) = colnames(ligand_target_matrix) %>% convert_alias_to_symbols(organism = organism) %>% make.names() + rownames(ligand_target_matrix) = rownames(ligand_target_matrix) %>% convert_alias_to_symbols(organism = organism) %>% make.names() + lr_network = lr_network %>% filter(ligand %in% colnames(ligand_target_matrix)) + ligand_target_matrix = ligand_target_matrix[, lr_network$ligand %>% unique()] +} +``` + +## Step 0.2: Read in SingleCellExperiment Objects + +In this vignette, sender and receiver cell types are in the same SingleCellExperiment object, which we will load here. In this vignette, we will load in a subset of the scRNAseq data of the MIS-C [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.8010790.svg)](https://doi.org/10.5281/zenodo.8010790). For the sake of demonstration, this subset only contains 3 cell types. These celltypes are some of the cell types that were found to be most interesting related to MIS-C according to Hoste et al. + +If you start from a Seurat object, you can convert it easily to a SingleCellExperiment via `sce = Seurat::as.SingleCellExperiment(seurat_obj, assay = "RNA")`. + +Because the NicheNet 2.0. networks are in the most recent version of the official gene symbols, we will make sure that the gene symbols used in the expression data are also updated (= converted from their "aliases" to official gene symbols). Afterwards, we will make them again syntactically valid. + +```{r} +sce = readRDS(url("https://zenodo.org/record/8010790/files/sce_subset_misc.rds")) +sce = alias_to_symbol_SCE(sce, "human") %>% makenames_SCE() +``` + +## Step 0.3: Prepare settings of the MultiNicheNet cell-cell communication analysis + +### Define in which metadata columns we can find the **group**, **sample** and **cell type** IDs + +In this case study, we want to study differences in cell-cell communication patterns between MIS-C patients (M), their healthy siblings (S) and adult patients with severe covid (A). The meta data columns that indicate this disease status is `MIS.C.AgeTier`. + +Cell type annotations are indicated in the `Annotation_v2.0` column, and the sample is indicated by the `ShortID` column. +If your cells are annotated in multiple hierarchical levels, we recommend using a high level in the hierarchy. This for 2 reasons: 1) MultiNicheNet focuses on differential expression and not differential abundance, and 2) there should be sufficient cells per sample-celltype combination. + +If you would have batch effects or covariates you can correct for, you can define this here as well. + +Important: for categorical covariates and batches, there should be at least one sample for every group-batch combination. If one of your groups/conditions lacks a certain level of your batch, you won't be able to correct for the batch effect because the model is then not able to distinguish batch from group/condition effects. + +Important: the column names of group, sample, cell type, batches and covariates should be syntactically valid (`make.names`) + +Important: All group, sample, cell type, batch and covariate names should be syntactically valid as well (`make.names`) (eg through `SummarizedExperiment::colData(sce)$ShortID = SummarizedExperiment::colData(sce)$ShortID %>% make.names()`) + +```{r} +sample_id = "ShortID" +group_id = "MIS.C.AgeTier" +celltype_id = "Annotation_v2.0" +covariates = NA +batches = NA +``` + +Sender and receiver cell types also need to be defined. Both are here all cell types in the dataset because we are interested in an All-vs-All analysis. + +```{r} +senders_oi = SummarizedExperiment::colData(sce)[,celltype_id] %>% unique() +receivers_oi = SummarizedExperiment::colData(sce)[,celltype_id] %>% unique() +``` + +If the user wants it, it is possible to use only a subset of senders and receivers. Senders and receivers can be entirely different, but also overlapping, or the same. If you don't use all the cell types in your data, we recommend to continue with a subset of your data. + +```{r} +sce = sce[, SummarizedExperiment::colData(sce)[,celltype_id] %in% c(senders_oi, receivers_oi)] +``` + +Now we will go to the first real step of the MultiNicheNet analysis + +# Step 1: Extract cell type abundance and expression information from receiver and sender cell types, and link this expression information for ligands of the sender cell types to the corresponding receptors of the receiver cell types + +Since MultiNicheNet will infer group differences at the sample level for each cell type (currently via Muscat - pseudobulking + EdgeR), we need to have sufficient cells per sample of a cell type, and this for both groups. In the following analysis we will set this minimum number of cells per cell type per sample at 10 (recommended minimum). + +```{r} +min_cells = 10 +``` + +Now we will calculate abundance and expression information for each cell type / sample / group combination with the following functions. In the output of this function, you can also find some 'Cell type abundance diagnostic plots' that will the users which celltype-sample combinations will be left out later on for DE calculation because the nr of cells is lower than de defined minimum defined here above. If too many celltype-sample combinations don't pass this threshold, we recommend to define your cell types in a more general way (use one level higher of the cell type ontology hierarchy) (eg TH17 CD4T cells --> CD4T cells). + +```{r} +abundance_info = get_abundance_info(sce = sce, sample_id = sample_id, group_id = group_id, celltype_id = celltype_id, min_cells = min_cells, senders_oi = senders_oi, receivers_oi = receivers_oi, batches = batches) +``` + +First, we will check the cell type abundance diagnostic plots. + +### Interpretation of cell type abundance information + +The first plot visualizes the number of cells per celltype-sample combination, and indicates which combinations are removed during the DE analysis because there are less than `min_cells` in the celltype-sample combination. + +```{r, fig.width=15, fig.height=12} +abundance_info$abund_plot_sample +``` +The red dotted line indicates the required minimum of cells as defined above in `min_cells`. We can see here that some sample-celltype combinations are left out. For the DE analysis in the next step, only cell types will be considered if there are at least two samples per group with a sufficient number of cells. + +__Important__: Based on the cell type abundance diagnostics, we recommend users to change their analysis settings if required (eg changing cell type annotation level, batches, ...), before proceeding with the rest of the analysis. + +as we can see here: no condition-specific cell types. In case you would have them... + +# Step 2: Perform genome-wide differential expression analysis of receiver and sender cell types to define DE genes between the conditions of interest. Based on this analysis, we can define the logFC/p-value of ligands in senders and receptors in receivers, and define the set of affected target genes in the receiver. + +Now we will go over to the multi-group, multi-sample differential expression (DE) analysis (also called 'differential state' analysis by the developers of Muscat). + +### Define the contrasts and covariates of interest for the DE analysis. + +Here, we want to compare each patient group to the other groups, so the MIS-C (M) group vs healthy control siblins (S) and adult COVID19 patients (A) etcetera. +To do this comparison, we need to set the following contrasts: + +```{r} +contrasts_oi = c("'M-(S+A)/2','S-(M+A)/2','A-(S+M)/2'") +``` + +__Very Important__ Note the format to indicate the contrasts! This formatting should be adhered to very strictly, and white spaces are not allowed! Check `?get_DE_info` for explanation about how to define this well. The most important things are that: each contrast is surrounded by single quotation marks, contrasts are separated by a comma without any whitespace, and alle contrasts together are surrounded by double quotation marks. If you compare against two groups, you should divide by 2, if you compare against three groups, you should divide by 3 etcetera. + +For downstream visualizations and linking contrasts to their main group, you need to run the following: + +```{r} +contrast_tbl = tibble(contrast = + c("M-(S+A)/2","S-(M+A)/2", "A-(S+M)/2"), + group = c("M","S","A")) +``` + +If you want to compare only two groups (eg M vs S), you can do like this: +`contrasts_oi = c("'M-S','S-M'") ` +`contrast_tbl = tibble(contrast = c("M-S","S-M"), group = c("M","S"))` + +### Perform the DE analysis for each cell type. + +First define expressed genes +```{r} +fraction_cutoff = 0.05 +min_sample_prop = 0.50 +frq_list = get_frac_exprs(sce = sce, sample_id = sample_id, celltype_id = celltype_id, group_id = group_id, batches = batches, min_cells = min_cells, fraction_cutoff = fraction_cutoff, min_sample_prop = min_sample_prop) +``` + +```{r} +DE_info = get_DE_info(sce = sce, sample_id = sample_id, group_id = group_id, celltype_id = celltype_id, batches = batches, covariates = covariates, contrasts_oi = contrasts_oi, min_cells = min_cells, expressed_df = frq_list$expressed_df) +``` + +### Check DE results + +Table with logFC and p-values for each gene-celltype-contrast: + +```{r} +DE_info$celltype_de$de_output_tidy %>% arrange(p_adj) %>% head() +``` + +```{r} +celltype_de = DE_info$celltype_de$de_output_tidy +``` + +In the next step, we will combine the DE information of senders and receivers by linking their ligands and receptors together: + +### Combine DE information for ligand-senders and receptors-receivers (similar to step1 - `abundance_expression_info$sender_receiver_info`) + +```{r} +sender_receiver_de = combine_sender_receiver_de( + sender_de = celltype_de, + receiver_de = celltype_de, + senders_oi = senders_oi, + receivers_oi = receivers_oi, + lr_network = lr_network +) +``` + +```{r} +sender_receiver_de %>% head(20) +``` + +# Calculate (pseudobulk) expression averages + +```{r} +abundance_expression_info = process_abundance_expression_info(sce = sce, sample_id = sample_id, group_id = group_id, celltype_id = celltype_id, min_cells = min_cells, senders_oi = senders_oi, receivers_oi = receivers_oi, lr_network = lr_network, batches = batches, frq_list = frq_list, abundance_info = abundance_info) +``` + +# Step 3: Predict NicheNet ligand activities and NicheNet ligand-target links based on these differential expression results + +## Define the parameters for the NicheNet ligand activity analysis + +Here, we need to define the thresholds that will be used to consider genes as differentially expressed or not (logFC, p-value, decision whether to use adjusted or normal p-value, minimum fraction of cells that should express a gene in at least one sample in a group, whether to use the normal p-values or empirical p-values). + +NicheNet ligand activity will then be calculated as the enrichment of predicted target genes of ligands in this set of DE genes compared to the genomic background. Here we choose for a minimum logFC of 0.50, maximum p-value of 0.05, and minimum fraction of expression of 0.05. + +```{r} +logFC_threshold = 0.50 +p_val_threshold = 0.05 +``` + +We will here choose for applying the p-value cutoff on the normal p-values, and not on the p-values corrected for multiple testing. This choice was made here because this dataset has only a few samples per group and we might have a lack of statistical power due to pseudobulking. In case of more samples per group, and a sufficient high number of DE genes per group-celltype (> 50), we would recommend using the adjusted p-values. + +```{r} +# p_val_adj = TRUE +p_val_adj = FALSE +``` + +For the NicheNet ligand-target inference, we also need to select which top n of the predicted target genes will be considered (here: top 250 targets per ligand). + +```{r} +top_n_target = 250 +``` + +The NicheNet ligand activity analysis can be run in parallel for each receiver cell type, by changing the number of cores as defined here. This is only recommended if you have many receiver cell type. + +```{r} +verbose = TRUE +cores_system = 8 +n.cores = min(cores_system, sender_receiver_de$receiver %>% unique() %>% length()) # use one core per receiver cell type +``` + +## Run the NicheNet ligand activity analysis + +(this might take some time) +```{r} +ligand_activities_targets_DEgenes = suppressMessages(suppressWarnings(get_ligand_activities_targets_DEgenes( + receiver_de = celltype_de, + receivers_oi = receivers_oi, + ligand_target_matrix = ligand_target_matrix, + logFC_threshold = logFC_threshold, + p_val_threshold = p_val_threshold, + p_val_adj = p_val_adj, + top_n_target = top_n_target, + verbose = verbose, + n.cores = n.cores +))) +``` + +Check the DE genes used for the activity analysis + +```{r} +ligand_activities_targets_DEgenes$de_genes_df %>% head(20) +``` + +Check the output of the activity analysis + +```{r} +ligand_activities_targets_DEgenes$ligand_activities %>% head(20) +``` + +# Step 4: Use the information collected above to prioritize all sender-ligand---receiver-receptor pairs. + +In the 3 previous steps, we calculated expression, differential expression and NicheNet activity information. Now we will combine these different types of information in one prioritization scheme. + +MultiNicheNet allows the user to define the weights of the following criteria to prioritize ligand-receptor interactions: + +* Upregulation of the ligand in a sender cell type and/or upregulation of the receptor in a receiver cell type - in the condition of interest. : `de_ligand` and `de_receptor` +* Sufficiently high expression levels of ligand and receptor in many samples of the same group (to mitigate the influence of outlier samples). : `frac_exprs_ligand_receptor` +* Cell-type and condition specific expression of the ligand in the sender cell type and receptor in the receiver cell type (to mitigate the influence of upregulated but still relatively weakly expressed ligands/receptors) : `exprs_ligand` and `exprs_receptor` +* High NicheNet ligand activity, to further prioritize ligand-receptor pairs based on their predicted effect of the ligand-receptor interaction on the gene expression in the receiver cell type : `activity_scaled` +* High relative abundance of sender and/or receiver in the condition of interest: `abund_sender` and `abund_receiver` (experimental feature - not recommended to give non-zero weights for default analyses) + +The different properties of the sender-ligand---receiver-receptor pairs can be weighted according to the user's preference and insight in the dataset at hand. + +## Prepare grouping objects + +Make necessary grouping data frame + +```{r} +sender_receiver_tbl = sender_receiver_de %>% dplyr::distinct(sender, receiver) + +metadata_combined = SummarizedExperiment::colData(sce) %>% tibble::as_tibble() + +if(!is.na(batches)){ + grouping_tbl = metadata_combined[,c(sample_id, group_id, batches)] %>% tibble::as_tibble() %>% dplyr::distinct() + colnames(grouping_tbl) = c("sample","group",batches) +} else { + grouping_tbl = metadata_combined[,c(sample_id, group_id)] %>% tibble::as_tibble() %>% dplyr::distinct() + colnames(grouping_tbl) = c("sample","group") +} + +``` + +Crucial note: grouping_tbl: group should be the same as in the contrast_tbl, and as in the expression info tables! Rename accordingly if this would not be the case. If you followed the guidelines of this tutorial closely, there should be no problem. + +## Run the prioritization + +```{r} + prioritization_tables = suppressMessages(generate_prioritization_tables( + sender_receiver_info = abundance_expression_info$sender_receiver_info, + sender_receiver_de = sender_receiver_de, + ligand_activities_targets_DEgenes = ligand_activities_targets_DEgenes, + contrast_tbl = contrast_tbl, + sender_receiver_tbl = sender_receiver_tbl, + grouping_tbl = grouping_tbl, + scenario = "regular", # all prioritization criteria will be weighted equally + fraction_cutoff = fraction_cutoff, + abundance_data_receiver = abundance_expression_info$abundance_data_receiver, + abundance_data_sender = abundance_expression_info$abundance_data_sender, + ligand_activity_down = FALSE # use only upregulatory ligand activities to prioritize -- if you want to consider max of up and down: set this to TRUE + )) +``` + +Check the output tables + +First: group-based summary table + +```{r} +prioritization_tables$group_prioritization_tbl %>% head(20) +``` + + +# Step 5: Add information on prior knowledge and expression correlation between LR and target expression. + +In multi-sample datasets, we have the opportunity to look whether expression of ligand-receptor across all samples is correlated with the expression of their by NicheNet predicted target genes. This is what we will do with the following line of code: + +```{r} +lr_target_prior_cor = lr_target_prior_cor_inference(prioritization_tables$group_prioritization_tbl$receiver %>% unique(), abundance_expression_info, celltype_de, grouping_tbl, prioritization_tables, ligand_target_matrix, logFC_threshold = logFC_threshold, p_val_threshold = p_val_threshold, p_val_adj = p_val_adj) +``` + +# Save all the output of MultiNicheNet + +To avoid needing to redo the analysis later. +All the output written down here is sufficient to make all in-built downstream visualizations. + +```{r} +path = "./" + +multinichenet_output = list( + celltype_info = abundance_expression_info$celltype_info, + celltype_de = celltype_de, + sender_receiver_info = abundance_expression_info$sender_receiver_info, + sender_receiver_de = sender_receiver_de, + ligand_activities_targets_DEgenes = ligand_activities_targets_DEgenes, + prioritization_tables = prioritization_tables, + grouping_tbl = grouping_tbl, + lr_target_prior_cor = lr_target_prior_cor + ) +multinichenet_output = make_lite_output(multinichenet_output) + +save = FALSE +if(save == TRUE){ + saveRDS(multinichenet_output, paste0(path, "multinichenet_output.rds")) + +} +``` + + +# Visualization of the results of the cell-cell communication analysis + +In a first instance, we will look at the broad overview of prioritized interactions via condition-specific Circos plots. + +## Circos plot of top-prioritized links + +We will look here at the top 50 predictions across all contrasts, senders, and receivers of interest. + +```{r} +prioritized_tbl_oi_all = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 50, rank_per_group = FALSE) +``` + + +```{r, fig.width=15, fig.height=12} +prioritized_tbl_oi = multinichenet_output$prioritization_tables$group_prioritization_tbl %>% + filter(id %in% prioritized_tbl_oi_all$id) %>% + distinct(id, sender, receiver, ligand, receptor, group) %>% left_join(prioritized_tbl_oi_all) +prioritized_tbl_oi$prioritization_score[is.na(prioritized_tbl_oi$prioritization_score)] = 0 + +senders_receivers = union(prioritized_tbl_oi$sender %>% unique(), prioritized_tbl_oi$receiver %>% unique()) %>% sort() + +colors_sender = RColorBrewer::brewer.pal(n = length(senders_receivers), name = 'Spectral') %>% magrittr::set_names(senders_receivers) +colors_receiver = RColorBrewer::brewer.pal(n = length(senders_receivers), name = 'Spectral') %>% magrittr::set_names(senders_receivers) + +circos_list = make_circos_group_comparison(prioritized_tbl_oi, colors_sender, colors_receiver) +``` + +## Visualization of scaled ligand-receptor pseudobulk products and ligand activity + +Now we will visualize per sample the scaled product of ligand and receptor expression. Samples that were left out of the DE analysis are indicated with a smaller dot (this helps to indicate the samples that did not contribute to the calculation of the logFC, and thus not contributed to the final prioritization) + +We will now check the top 50 interactions specific for the MIS-C group + +```{r} +group_oi = "M" +``` + +```{r, fig.height=13, fig.width=15} +prioritized_tbl_oi_M_50 = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 50, groups_oi = group_oi) +``` + +```{r, fig.height=13, fig.width=18} +plot_oi = make_sample_lr_prod_activity_plots(multinichenet_output$prioritization_tables, prioritized_tbl_oi_M_50) +plot_oi +``` + +As a further help for further prioritization, we can assess the level of curation of these LR pairs as defined by the Intercellular Communication part of the Omnipath database + +```{r} +prioritized_tbl_oi_M_50_omnipath = prioritized_tbl_oi_M_50 %>% inner_join(lr_network_all) +``` + +```{r, fig.height=13, fig.width=16} +plot_oi = make_sample_lr_prod_activity_plots_Omnipath(multinichenet_output$prioritization_tables, prioritized_tbl_oi_M_50_omnipath) +plot_oi +``` + +As you can see, the HEBP1-FPR2 interaction has no Omnipath DB scores. This is because this LR pair was not documented by the Omnipath LR database. Instead it was documented by the original NicheNet LR network (source: Guide2Pharmacology) as can be seen in the table. + + + +Further note: Typically, there are way more than 50 differentially expressed and active ligand-receptor pairs per group across all sender-receiver combinations. Therefore it might be useful to zoom in on specific cell types as senders/receivers: + +Eg M_Monocyte_CD16 as receiver: + +```{r, fig.height=13, fig.width=15} +prioritized_tbl_oi_M_50 = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 50, groups_oi = group_oi, receivers_oi = "M_Monocyte_CD16") + +plot_oi = make_sample_lr_prod_activity_plots(multinichenet_output$prioritization_tables, prioritized_tbl_oi_M_50) +plot_oi +``` + +Eg M_Monocyte_CD16 as sender: + +```{r, fig.height=13, fig.width=15} +prioritized_tbl_oi_M_50 = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 50, groups_oi = group_oi, senders_oi = "M_Monocyte_CD16") + +plot_oi = make_sample_lr_prod_activity_plots(multinichenet_output$prioritization_tables, prioritized_tbl_oi_M_50) +plot_oi +``` + +You can make these plots also for the other groups, like we will illustrate now for the S group + +```{r} +group_oi = "S" +``` + +```{r, fig.height=13, fig.width=15} +prioritized_tbl_oi_S_50 = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 50, groups_oi = group_oi) + +plot_oi = make_sample_lr_prod_activity_plots(multinichenet_output$prioritization_tables, prioritized_tbl_oi_S_50) +plot_oi +``` + +__Note__: Use `make_sample_lr_prod_activity_batch_plots` if you have batches and want to visualize them on this plot! + +## Visualization of ligand-activity, ligand-target links, and target gene expression + +In another type of plot, we can visualize the ligand activities for a group-receiver combination, and show the predicted ligand-target links, and also the expression of the predicted target genes across samples. + +For this, we now need to define a receiver cell type of interest. As example, we will take `L_T_TIM3._CD38._HLADR.` cells as receiver, and look at the top 20 senderLigand-receiverReceptor pairs with these cells as receiver. + +```{r} +group_oi = "M" +receiver_oi = "L_T_TIM3._CD38._HLADR." +prioritized_tbl_oi_M_10 = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 10, groups_oi = group_oi, receivers_oi = receiver_oi) +``` + +```{r, fig.width=20, fig.height=10} +combined_plot = make_ligand_activity_target_plot(group_oi, receiver_oi, prioritized_tbl_oi_M_10, multinichenet_output$prioritization_tables, multinichenet_output$ligand_activities_targets_DEgenes, contrast_tbl, multinichenet_output$grouping_tbl, multinichenet_output$celltype_info, ligand_target_matrix, plot_legend = FALSE) +combined_plot +``` + +__Note__ Use `make_DEgene_dotplot_pseudobulk_batch` if you want to indicate the batch of each sample to the plot + +What if there is a specific ligand you are interested in? + +```{r} +group_oi = "M" +receiver_oi = "L_T_TIM3._CD38._HLADR." +ligands_oi = c("IFNG","IL15") +prioritized_tbl_ligands_oi = get_top_n_lr_pairs(multinichenet_output$prioritization_tables, 10000, groups_oi = group_oi, receivers_oi = receiver_oi) %>% filter(ligand %in% ligands_oi) # ligands should still be in the output tables of course +``` + +```{r, fig.width=20, fig.height=10} +combined_plot = make_ligand_activity_target_plot(group_oi, receiver_oi, prioritized_tbl_ligands_oi, multinichenet_output$prioritization_tables, multinichenet_output$ligand_activities_targets_DEgenes, contrast_tbl, multinichenet_output$grouping_tbl, multinichenet_output$celltype_info, ligand_target_matrix, plot_legend = FALSE) +combined_plot +``` + +## Visualization of expression-correlated target genes of ligand-receptor pairs + +Before, we had calculated the correlation in expression between ligand-receptor pairs and DE genes. Now we will filter out correlated ligand-receptor --> target links that both show high expression correlation (spearman or activity > 0.50 in this example) and have some prior knowledge to support their link. + +```{r} +group_oi = "M" +receiver_oi = "M_Monocyte_CD16" +lr_target_prior_cor_filtered = multinichenet_output$lr_target_prior_cor %>% inner_join(multinichenet_output$ligand_activities_targets_DEgenes$ligand_activities %>% distinct(ligand, target, direction_regulation, contrast)) %>% inner_join(contrast_tbl) %>% filter(group == group_oi, receiver == receiver_oi) +lr_target_prior_cor_filtered_up = lr_target_prior_cor_filtered %>% filter(direction_regulation == "up") %>% filter( (rank_of_target < top_n_target) & (pearson > 0.50 | spearman > 0.50)) +lr_target_prior_cor_filtered_down = lr_target_prior_cor_filtered %>% filter(direction_regulation == "down") %>% filter( (rank_of_target < top_n_target) & (pearson < -0.50 | spearman < -0.50)) # downregulation -- negative correlation +lr_target_prior_cor_filtered = bind_rows(lr_target_prior_cor_filtered_up, lr_target_prior_cor_filtered_down) +``` + +Now we will visualize the top correlated target genes for the LR pairs that are also in the top 50 LR pairs discriminating the groups from each other: + +```{r} +prioritized_tbl_oi = get_top_n_lr_pairs(prioritization_tables, 50, groups_oi = group_oi, receivers_oi = receiver_oi) +``` + + +```{r, fig.width=26, fig.height=16} +lr_target_correlation_plot = make_lr_target_correlation_plot(multinichenet_output$prioritization_tables, prioritized_tbl_oi, lr_target_prior_cor_filtered , multinichenet_output$grouping_tbl, multinichenet_output$celltype_info, receiver_oi,plot_legend = FALSE) +lr_target_correlation_plot$combined_plot +``` + +You can also visualize the expression correlation in the following way for a selected LR pair and their targets: + +```{r, fig.width=21, fig.height=6} +ligand_oi = "IFNG" +receptor_oi = "IFNGR2" +sender_oi = "L_T_TIM3._CD38._HLADR." +receiver_oi = "M_Monocyte_CD16" +lr_target_scatter_plot = make_lr_target_scatter_plot(multinichenet_output$prioritization_tables, ligand_oi, receptor_oi, sender_oi, receiver_oi, multinichenet_output$celltype_info, multinichenet_output$grouping_tbl, lr_target_prior_cor_filtered) +lr_target_scatter_plot +``` + +## Intercellular regulatory network systems view + +In the plots before, we demonstrated that some DE genes have both expression correlation and prior knowledge support to be downstream of ligand-receptor pairs. Interestingly, some target genes can be ligands or receptors themselves. This illustrates that cells can send signals to other cells, who as a response to these signals produce signals themselves to feedback to the original sender cells, or who will effect other cell types. + +As last plot, we can generate a 'systems' view of these intercellular feedback and cascade processes than can be occuring between the different cell populations involved. In this plot, we will draw links between ligands of sender cell types their ligand/receptor-annotated target genes in receiver cell types. So links are ligand-target links (= gene regulatory links) and not ligand-receptor protein-protein interactions! + +```{r} +prioritized_tbl_oi = get_top_n_lr_pairs(prioritization_tables, 250, rank_per_group = FALSE) + +lr_target_prior_cor_filtered = prioritization_tables$group_prioritization_tbl$group %>% unique() %>% lapply(function(group_oi){ + lr_target_prior_cor_filtered = multinichenet_output$lr_target_prior_cor %>% inner_join(multinichenet_output$ligand_activities_targets_DEgenes$ligand_activities %>% distinct(ligand, target, direction_regulation, contrast)) %>% inner_join(contrast_tbl) %>% filter(group == group_oi) + lr_target_prior_cor_filtered_up = lr_target_prior_cor_filtered %>% filter(direction_regulation == "up") %>% filter( (rank_of_target < top_n_target) & (pearson > 0.50 | spearman > 0.50)) + lr_target_prior_cor_filtered_down = lr_target_prior_cor_filtered %>% filter(direction_regulation == "down") %>% filter( (rank_of_target < top_n_target) & (pearson < -0.50 | spearman < -0.50)) + lr_target_prior_cor_filtered = bind_rows(lr_target_prior_cor_filtered_up, lr_target_prior_cor_filtered_down) +}) %>% bind_rows() + +lr_target_df = lr_target_prior_cor_filtered %>% distinct(group, sender, receiver, ligand, receptor, id, target, direction_regulation) +``` + +```{r} +network = infer_intercellular_regulatory_network(lr_target_df, prioritized_tbl_oi) +network$links %>% head() +network$nodes %>% head() +``` + +```{r, fig.width=30, fig.height=12} +colors_sender["L_T_TIM3._CD38._HLADR."] = "pink" # the original yellow background with white font is not very readable +network_graph = visualize_network(network, colors_sender) +network_graph$plot +``` + +An interesting follow-up question would be to find nodes that are at the top of the hierarchy in this network, and could thus explain the cascade of intercellular communication events we observe. + +Using PageRank on the opposite-directionality graph can help identify important nodes that act as sources in the hierarchical structure. By running PageRank on the reversed graph (where edges point in the opposite direction), you can potentially pinpoint nodes that have significant influence or act as sources of information in your network. + +We will exemplify for the M-group + +```{r,} +library(igraph) + +# Create the normal graph (normal directionality) +normal_edges <- get.data.frame(network_graph$network_igraph) %>% filter(group == "M") %>% .[, c("from", "to")] +normal_graph <- graph_from_data_frame(normal_edges, directed = TRUE) +# Calculate PageRank on the normal graph +pagerank_scores <- page_rank(normal_graph)$vector %>% sort(decreasing = TRUE) +pr_score_tbl = tibble(node = names(pagerank_scores), score_normal = pagerank_scores) %>% arrange(-score_normal) + +# Create the reversed graph (opposite directionality) +reversed_edges <- get.data.frame(network_graph$network_igraph) %>% filter(group == "M") %>% .[, c("to", "from")] # Swap 'from' and 'to' columns +reversed_graph <- graph_from_data_frame(reversed_edges, directed = TRUE) +# Calculate PageRank on the reversed graph +reversed_pagerank_scores <- page_rank(reversed_graph)$vector %>% sort(decreasing = TRUE) + +pr_score_tbl = pr_score_tbl %>% inner_join( + tibble(node = names(reversed_pagerank_scores), score_reverse = reversed_pagerank_scores) +) %>% arrange(-score_reverse) +pr_score_tbl +``` + +Let's zoom now in on the M-group network to check this + +```{r} +network = infer_intercellular_regulatory_network(lr_target_df %>% filter(group == "M"), prioritized_tbl_oi %>% filter(group == "M")) +network$links %>% head() +network$nodes %>% head() +``` + +```{r, fig.width=30, fig.height=20} +colors_sender["L_T_TIM3._CD38._HLADR."] = "pink" # the original yellow background with white font is not very readable +network_graph = visualize_network(network, colors_sender) +network_graph$plot +``` + +Interestingly, we can also use this network to further prioritize differential CCC interactions. Here we will assume that the most important LR interactions are the ones that are involved in this intercellular regulatory network. We can get these interactions as follows: +```{r} +network$prioritized_lr_interactions +``` + +```{r, fig.width=30, fig.height=12} +prioritized_tbl_oi_network = prioritized_tbl_oi %>% inner_join(network$prioritized_lr_interactions) +prioritized_tbl_oi_network +``` + +Visualize now the expression and activity of these interactions for the M group +```{r} +group_oi = "M" +``` + +```{r, fig.height=16, fig.width=18} +prioritized_tbl_oi_M = prioritized_tbl_oi_network %>% filter(group == group_oi) + +plot_oi = make_sample_lr_prod_activity_plots(multinichenet_output$prioritization_tables, prioritized_tbl_oi_M) +plot_oi +``` + +## Investingating the data sources behind the ligand-receptor interactions + +Now link the prioritized interactions to all the information of the ligand-receptor network. This can help to check which are the specific data sources behind the inferred interactions, which can aid in prioritizing interactions for further validation. + +Check OmnipathR for more information with respect to how to interpret the different metrics specific for the interactions coming from the Omnipath database. +```{r} +prioritized_tbl_oi_M %>% inner_join(lr_network_all) %>% arrange(-n_resources) # sort interactions based on nr of resources that document them +```