Update basic vignette of MNN v2

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 = "[email protected]",
   role = c("aut", "cre"))
 Maintainer: Robin Browaeys <[email protected]>
-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

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)

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)
+}
 

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(

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 -------------------------------------------------------------------------------------------------------------------------------------------