Fix issue disconcordance wrapper/step-by-step

Author: browaeysrobin

.DS_Store

@@ -1,7 +1,7 @@
 Package: multinichenetr
 Type: Package
 Title: MultiNicheNet: a flexible framework for differential cell-cell communication analysis from multi-sample multi-condition single-cell transcriptomics data
-Version: 2.0.0
+Version: 2.0.1
 Author: person("Robin", "Browaeys", email = "robin.browaeys@ugent.be",
   role = c("aut", "cre"))
 Maintainer: Robin Browaeys <robin.browaeys@ugent.be>

DESCRIPTION

@@ -355,11 +355,19 @@ multi_nichenet_analysis = function(sce,
   # do not do this -- this could give errors if only interested in one contrast but multiple groups
 
   if(verbose == TRUE){
-    print("Make diagnostic abundance plots + define expressed genes")
+    print("Cell type & sample filtering")
   }
   ## check abundance info
 
-  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)
 
   ## check for condition-specific cell types
   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)
@@ -377,6 +385,9 @@ multi_nichenet_analysis = function(sce,
   print("condition-specific celltypes:")
   print(condition_specific_celltypes)
 
+  condition_specific_celltypes_senders = condition_specific_celltypes %>% intersect(senders_oi)
+  condition_specific_celltypes_receivers = condition_specific_celltypes %>% intersect(receivers_oi)
+  
   print("absent celltypes:")
   print(absent_celltypes)
 
@@ -388,9 +399,28 @@ multi_nichenet_analysis = function(sce,
   sce = sce[, SummarizedExperiment::colData(sce)[,celltype_id] %in% retained_celltypes]
 
   ## define expressed genes
+  if(verbose == TRUE){
+    print("Gene filtering")
+  }
   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)
 
-  ### Perform the DE analysis ----------------------------------------------------------------
+  ## filter out non-expressed genes
+  genes_oi = frq_list$expressed_df %>% 
+    filter(expressed == TRUE) %>% pull(gene) %>% unique() 
+  sce = sce[genes_oi, ]
+  
+  ## calculate PB expresion
+  if(verbose == TRUE){
+    print("Calculate normalized average and pseudobulk expression")
+  }
+  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 = union(senders_oi, condition_specific_celltypes_senders), 
+    receivers_oi = union(receivers_oi, condition_specific_celltypes_receivers), 
+    lr_network = lr_network, batches = batches, frq_list = frq_list, abundance_info = abundance_info
+    )
+  
+  ## Perform the DE analysis ----------------------------------------------------------------
 
   if(verbose == TRUE){
     print("Calculate differential expression for all cell types")
@@ -432,9 +462,9 @@ multi_nichenet_analysis = function(sce,
 
   # print(celltype_de %>% dplyr::group_by(cluster_id, contrast) %>% dplyr::filter(p_adj <= p_val_threshold & abs(logFC) >= logFC_threshold) %>% dplyr::count() %>% dplyr::arrange(-n))
 
-  senders_oi = celltype_de$cluster_id %>% unique()
-  receivers_oi = celltype_de$cluster_id %>% unique()
-  genes_oi = celltype_de$gene %>% unique()
+  senders_oi = intersect(celltype_de$cluster_id %>% unique(), senders_oi)
+  receivers_oi = intersect(celltype_de$cluster_id %>% unique(), receivers_oi)
+  genes_oi = intersect(genes_oi, celltype_de$gene %>% unique())
 
   retained_celltypes = c(senders_oi, receivers_oi) %>% unique()
   retained_celltypes = c(retained_celltypes, condition_specific_celltypes) 
@@ -453,12 +483,6 @@ multi_nichenet_analysis = function(sce,
   ))
   sender_receiver_tbl = sender_receiver_de %>% dplyr::distinct(sender, receiver)
 
-  ### Receiver abundance plots + Calculate expression information
-  if(verbose == TRUE){
-    print("Calculate normalized average and pseudobulk expression")
-  }
-  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 = union(senders_oi, condition_specific_celltypes), receivers_oi = union(receivers_oi, condition_specific_celltypes), lr_network = lr_network, batches = batches, frq_list = frq_list, abundance_info = abundance_info)
-  
   metadata_combined = SummarizedExperiment::colData(sce) %>% tibble::as_tibble()
 
   if(!is.na(batches)){
@@ -470,10 +494,13 @@ multi_nichenet_analysis = function(sce,
   }
 
   rm(sce)
-  ### Use the DE analysis for defining the DE genes in the receiver cell type and perform NicheNet ligand activity and ligand-target inference ----------------------------------------------------------------
+  ## Use the DE analysis for defining the DE genes in the receiver cell type and perform NicheNet ligand activity and ligand-target inference ----------------------------------------------------------------
   if(verbose == TRUE){
     print("Calculate NicheNet ligand activities and ligand-target links")
   }
+  
+  n.cores = min(n.cores, celltype_de$cluster_id %>% unique() %>% length()) 
+  
   ligand_activities_targets_DEgenes = suppressMessages(suppressWarnings(get_ligand_activities_targets_DEgenes(
     receiver_de = celltype_de,
     receivers_oi = receivers_oi,
@@ -486,7 +513,8 @@ multi_nichenet_analysis = function(sce,
     n.cores = n.cores
   )))
 
-  ### Combine the three types of information calculated above to prioritize ligand-receptor interactions ----------------------------------------------------------------
+  
+  ## Combine the three types of information calculated above to prioritize ligand-receptor interactions ----------------------------------------------------------------
   if(verbose == TRUE){
     print("Combine all the information in prioritization tables")
   }
@@ -505,7 +533,7 @@ multi_nichenet_analysis = function(sce,
     ligand_activity_down = ligand_activity_down # use only upregulatory ligand activities to prioritize
   ))
 
-  # Prepare Unsupervised analysis of samples! ------------------------------------------------------------------------------------------------------------
+  ## Prepare Unsupervised analysis of samples! ------------------------------------------------------------------------------------------------------------
 
   if(return_lr_prod_matrix == TRUE){
 

R/pipeline.R

@@ -33,10 +33,9 @@ In this vignette, you can learn how to perform a MultiNicheNet analysis when you
 As the best solution, we recommend - if possible and biologically sensible to annotate your cells at one level higher in the cell type annotation hierarchy (eg sub cell type can be seen as a different cell state). Hereby, you will not have this problem of condition-specific cell types anymore. 
 
 A a second solution, we propose the procedure as showcased in this vignette. This consists of 3 consecutive workflows: 
-
 *1) General Workflow: This is the regular MultiNicheNet approach. In the general workflow, condition-specific cell types are excluded from the analysis. This approach focuses on achieving the best prioritization for other cell types that are present across different conditions, allowing for a standard DE analysis. The exclusion of condition-specific cell types simplifies the analysis but at the expense of potentially missing out on valuable insights that these unique cell types could provide. Therefore you can go further with the next two workflows: 
 *2) Workflow A (Sender Cell Types): This workflow adapts the analysis to include condition-specific cell types as sender cell types. It circumvents the need for DE analysis of the ligand by only focusing the ligand-prioritization on ligand cell type specificity and ligand activity. 
-*3) Workflow B (Receiver Cell Types): Workflow B allows for the inclusion of condition-specific cell types as receiver cell types by forgoing the DE analysis of the receiver, meaning that we cannot consider receptor differential expresion and ligand activity. Instead, receptor prioritization is based solely on the receptor cell type specificity. 
+*3) Workflow B (Receiver Cell Types): Workflow B allows for the inclusion of condition-specific cell types as receiver cell types by forgoing the DE analysis of the receiver, meaning that we cannot consider receptor differential expresion and ligand activity. Instead, receptor prioritization is based solely on the receptor cell type specificity.
 
 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 and patients with adult COVID-19. 
@@ -60,7 +59,6 @@ The Nichenet v2 networks and matrices for both mouse and human can be downloaded
 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"
 ```

vignettes/condition_specific_celltype_MISC.Rmd

@@ -420,6 +420,9 @@ We will now check the top 50 interactions specific for the MIS-C group
 
 ```{r}
 group_oi = "M"
+
+print(prioritized_tbl_oi_all)
+
 ```
 
 ```{r}