In saeyslab/multinichenetr: MultiNicheNet: a flexible framework for differential cell-cell communication analysis from multi-sample multi-condition single-cell transcriptomics data

Prepare visualizations

library(tidyverse)
library(nichenetr)
library(multinichenetr)

multinichenet_output = readRDS("../output/MNN_BCA_preEvspreNE_final_output.rds")

sample_id = "sample_id"
group_id = "expansion_timepoint"
celltype_id = "subType"
covariates = NA 
batches = NA 

contrasts_oi = c("'PreE-PreNE','PreNE-PreE'")
contrast_tbl = tibble(contrast =
                        c("PreE-PreNE", "PreNE-PreE"),
                      group = c("PreE", "PreNE"))
logFC_threshold = 0.50
p_val_threshold = 0.05
fraction_cutoff = 0.05
top_n_target = 250

Make cell type names more readable

mutate_cond <- function(.data, condition, ..., envir = parent.frame()) {
  condition <- eval(substitute(condition), .data, envir)
  .data[condition, ] <- .data[condition, ] %>% mutate(...)
  .data
}

celltype_labels = c(
  "macrophages" = "Macrophages",
  "Fibroblast" = "Fibroblasts",
  "Mast_cell" = "Mast cells",
  "CD4T" = "CD4 T cells",
  "CD8T" = "CD8 T cells",
  "CD4REG" = "CD4 TREGs",
  "NK" = "NK cells",
  "gdT" = "gd T cells",
  "Cancer_cell" = "Cancer cells",
  "B_cell" = "B cells",
  "monocytes" = "Monocytes",
  "Endothelial_cell" = "Endothelial cells"
)

for(i in seq(length(celltype_labels))){
  old_name = celltype_labels[i] %>% names()
  new_name = celltype_labels[i]
  multinichenet_output$prioritization_tables$group_prioritization_tbl = multinichenet_output$prioritization_tables$group_prioritization_tbl %>% mutate_cond(sender == old_name, sender = new_name) %>% mutate_cond(receiver == old_name, receiver = new_name)
  multinichenet_output$prioritization_tables$group_prioritization_table_source = multinichenet_output$prioritization_tables$group_prioritization_table_source %>% mutate_cond(sender == old_name, sender = new_name) %>% mutate_cond(receiver == old_name, receiver = new_name)
  multinichenet_output$prioritization_tables$sample_prioritization_tbl = multinichenet_output$prioritization_tables$sample_prioritization_tbl %>% mutate_cond(sender == old_name, sender = new_name) %>% mutate_cond(receiver == old_name, receiver = new_name)
  multinichenet_output$prioritization_tables$ligand_activities_target_de_tbl = multinichenet_output$prioritization_tables$ligand_activities_target_de_tbl %>% mutate_cond(receiver == old_name, receiver = new_name)
  multinichenet_output$ligand_activities_targets_DEgenes$ligand_activities = multinichenet_output$ligand_activities_targets_DEgenes$ligand_activities %>% mutate_cond(receiver == old_name, receiver = new_name)
  multinichenet_output$ligand_activities_targets_DEgenes$de_genes_df = multinichenet_output$ligand_activities_targets_DEgenes$de_genes_df %>% mutate_cond(receiver == old_name, receiver = new_name)
  multinichenet_output$celltype_info$pb_df = multinichenet_output$celltype_info$pb_df %>% mutate(celltype = as.character(celltype)) %>% mutate_cond(celltype == old_name, celltype = new_name)
  multinichenet_output$lr_target_prior_cor = multinichenet_output$lr_target_prior_cor %>% mutate_cond(sender == old_name, sender = new_name) %>% mutate_cond(receiver == old_name, receiver = new_name)

} 

### read in NicheNet model
lr_network_all = readRDS("../../../NicheNet_V2/networks/data/ligand_receptor/lr_network_human_allInfo_30112033.rds") %>% mutate(ligand = convert_alias_to_symbols(ligand, organism = "human"), receptor = convert_alias_to_symbols(receptor, organism = "human"))

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("../../../NicheNet_V2/model_construction/models/ligand_target_matrix_nsga2r_final.rds")
colnames(ligand_target_matrix) = colnames(ligand_target_matrix) %>% convert_alias_to_symbols(organism = "human") %>% make.names()
rownames(ligand_target_matrix) = rownames(ligand_target_matrix) %>% convert_alias_to_symbols(organism = "human") %>% make.names()
lr_network = lr_network %>% filter(ligand %in% colnames(ligand_target_matrix))
ligand_target_matrix = ligand_target_matrix[, lr_network$ligand %>% unique()]

Interpreting the MultiNicheNet analysis output

Visualization of differential cell-cell interactions

Summarizing ChordDiagram circos plots

In a first instance, we will look at the broad overview of prioritized interactions via condition-specific Chordiagram circos plots. The aim of this visualizatin is to provide a summary of the top prioritized senderLigand-receiverReceptor interactions per condition (between all cell types or between cell type pairs of interest).

We will look here at the top 50 predictions across all contrasts, senders, and receivers of interest.

prioritized_tbl_oi_all = get_top_n_lr_pairs(
  multinichenet_output$prioritization_tables, 
  top_n = 50, 
  rank_per_group = FALSE
  )

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 = c(RColorBrewer::brewer.pal(n = length(senders_receivers), name = 'Spectral'), "black") %>% magrittr::set_names(senders_receivers)
colors_receiver = c(RColorBrewer::brewer.pal(n = length(senders_receivers), name = 'Spectral'), "black") %>% magrittr::set_names(senders_receivers)

circos_list = make_circos_group_comparison(prioritized_tbl_oi, colors_sender, colors_receiver)

Whereas these ChordDiagram circos plots show the most specific interactions per group, they don't give insights into the data behind these predictions. Because inspecting the data behind the prioritization is recommended to decide on which interactions to validate, we created several functionalities to do this.

Therefore we will now generate "interpretable bubble plots" that indicate the different prioritization criteria used in MultiNicheNet.

Interpretable bubble plots

In the next type of plots, we will visualize the following prioritization criteria used in MultiNicheNet: 1) differential expression of ligand and receptor: the per-sample scaled product of normalized ligand and receptor pseudobulk expression 2) the scaled ligand activities * 3) cell-type specificity of ligand and receptor.

As a further help for users to further prioritize, we also visualize: the condition-average of the fraction of cells expressing the ligand and receptor in the cell types of interest the level of curation of these LR pairs as defined by the Intercellular Communication part of the Omnipath database (https://omnipathdb.org/)

We will create this plot for PreE group specific interactions of the overall top50 interactions that we visualized in the Circos Chorddiagrams above:

group_oi = "PreE"

prioritized_tbl_oi_PreE = prioritized_tbl_oi_all %>% 
  filter(group == group_oi)

plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
  multinichenet_output$prioritization_tables, 
  prioritized_tbl_oi_PreE %>% inner_join(lr_network_all)
  )
plot_oi

Some notes about this plot: Samples that were left out of the DE analysis (because too few cells in that celltype-sample combination) 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 you can see, the HEBP1-FPR2 interaction does not have 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 (lr_network_all %>% filter(ligand == "HEBP1" & receptor == "FPR2")).

We encourage users to make these plots also for the other groups, like we will do now first for the PreNE group

group_oi = "PreNE"

prioritized_tbl_oi_PreNE_50 = prioritized_tbl_oi_all %>% 
  filter(group == group_oi) 

plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
  multinichenet_output$prioritization_tables, 
  prioritized_tbl_oi_PreNE_50 %>% inner_join(lr_network_all)
)
plot_oi

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:

We will illustrate this for the "Macrophages" cell type as receiver in the PreE group:

group_oi = "PreE"
prioritized_tbl_oi_PreE = get_top_n_lr_pairs(
  multinichenet_output$prioritization_tables, 
  50, 
  groups_oi = group_oi, 
  receivers_oi = "Macrophages"
  ) 

plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
  multinichenet_output$prioritization_tables, 
  prioritized_tbl_oi_PreE %>% inner_join(lr_network_all)
  )
plot_oi

And now as sender:

prioritized_tbl_oi_PreE = get_top_n_lr_pairs(
  multinichenet_output$prioritization_tables, 
  50, 
  groups_oi = group_oi, 
  senders_oi = "Macrophages")

plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
  multinichenet_output$prioritization_tables, 
  prioritized_tbl_oi_PreE %>% inner_join(lr_network_all))
plot_oi

These two types of plots created above (Circos ChordDiagram and Interpretable Bubble Plot) for the most strongly prioritized interactions are the types of plot you should always create and inspect as an end-user.

The plots that we will discuss in the rest of the vignette are more optional, and can help to dive more deeply in the data. They are however not as necessary as the plots above.

So, let's now continue with more detailed plots and downstream functionalities:

Intercellular regulatory network inference and visualization

In the plots above, we showed some of the prioritized interactions, and focused on their expression and activity. These interactions were visualized as independent interactions. However, they are likely not functioning independently in a complex multicellular biological system: cells can send signals to other cells, who as a response to these signals produce extracellular signals themselves to give feedback to the original sender cells, or to propogate the signal to other cell types ("cascade"). In other words: ligands from cell type A may induce the expression of ligands and receptors in cell type B. These ligands and receptors can then be involved in other interactions towards cell type A and interactions towards cell type C. Etc.

Because one of the elements of MultiNicheNet is the ligand activity and ligand-target inference part of NicheNet, we can actually infer the predicted ligand/receptor-encoding target genes of prioritized ligand-receptor interactions. And as a result, we can get this type of functional insight in the biological system of interest, which we will demonstrate now.

First, we will showcase how to do this by considering target genes supported by NicheNet's prior knowledge solely

Without filtering of target genes based on LR-target expression correlation (for demonstration purposes only)

First: get the target genes of prioritized ligand-receptor pairs (here focused on the overall top50 prioritized LR pairs that were visualized in the Circos ChordDiagrams above)

lr_target_prior = prioritized_tbl_oi_all %>% inner_join(
        multinichenet_output$ligand_activities_targets_DEgenes$ligand_activities %>%
          distinct(ligand, target, direction_regulation, contrast) %>% inner_join(contrast_tbl) %>% ungroup() 
        ) 
lr_target_df = lr_target_prior %>% distinct(group, sender, receiver, ligand, receptor, id, target, direction_regulation) 

Second, subset on ligands/receptors as target genes

lr_target_df %>% filter(target %in% union(lr_network$ligand, lr_network$receptor))

Whereas these code blocks are just to demonstrate that this type of information is available in MultiNicheNet, the next block of code will infer the systems-wide intercellular regulatory network automatically:

network = infer_intercellular_regulatory_network(lr_target_df, prioritized_tbl_oi_all)
network$links %>% head()
network$nodes %>% head()

And this network can be visualized here in R by running:

colors_sender["Endothelial cells"] = "pink" # the  original yellow background with white font is not very readable
colors_sender["Fibroblasts"] = "limegreen" # the  original yellow background with white font is not very readable
network_graph = visualize_network(network, colors_sender)
network_graph$plot

As you can see here: we can see see here that several prioritized ligands seem to be regulated by other prioritized ligands! But, it may be challenging sometimes to discern individual links when several interactions are shown. Therefore, inspection of the underlying data tables (network$links and network$nodes) may be necessary to discern individual interactions. It is also suggested to export these data tables into more sophisticated network visualization tools (e.g., CytoScape) for better inspection of this network.

To inspect interactions involving specific ligands, such as IFNG as example, we can run the following code:

network$nodes %>% filter(gene == "IFNG")

IFNG as regulating ligand:

network$links %>% filter(sender_ligand == "CD4 T cells_IFNG" & direction_regulation == "up" & group == "PreE")

IFNG as regulated target:

network$links %>% filter(receiver_target == "CD4 T cells_IFNG" & direction_regulation == "up" & group == "PreE")

Ligand- and receptor-encoding target genes that were shown here are predicted as target genes of ligands based on prior knowledge. However, it is uncertain whether they are also potentially active in the system under study: e.g., it is possible that some genes are regulated by their upstream ligand only in cell types that are not studied in this context. To increase the chance that inferred ligand-target links are potentially active, 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. This is under the assumption that target genes that show across-sample expression correlation with their upstream ligand-receptor pairs may be more likely to be true active target genes than target genes that don’t show this pattern. This correlation was calculated in the (optional) step 7 of the MultiNicheNet analysis.

In the next subsection of the inference of intercellular regulator networks, we will showcase how to consider target genes that are both supported by NicheNet's prior knowledge and expression correlation.

With filtering of target genes based on LR-target expression correlation (recommended for analysis in practice)

Now, we will filter out correlated ligand-receptor --> target links that both show high expression correlation (pearson correlation > 0.33 in this example) and have some prior knowledge to support their link.

lr_target_prior_cor_filtered = 
  multinichenet_output$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.33))

    lr_target_prior_cor_filtered_down = lr_target_prior_cor_filtered %>% 
      filter(direction_regulation == "down") %>% 
      filter( (rank_of_target < top_n_target) & (pearson < -0.33))
    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) 

network = infer_intercellular_regulatory_network(lr_target_df, prioritized_tbl_oi_all)
network$links %>% head()
network$nodes %>% head()

network_graph = visualize_network(network, colors_sender)
network_graph$plot

As can be expected, we see fewer links here than in the previously generated intercellular regulatory network. The links that are not present anymore in this network are those ligand-target links that are not supported by high across-sample expression correlation. In conclusion, the links visualized here are the most trustworthy ones, since they are both supported by prior knowledge and expression correlation.

Interestingly, ligands/receptors visualized in this network can be considered as additionally prioritized because they are not only a prioritized ligand/receptor but also a target gene of another prioritized ligand-receptor interaction! So, we can also use this network to further prioritize differential CCC interactions. We can get these interactions as follows:

network$prioritized_lr_interactions

prioritized_tbl_oi_network = prioritized_tbl_oi_all %>% inner_join(
  network$prioritized_lr_interactions)
prioritized_tbl_oi_network

Visualize now the expression and activity of these interactions for the PreE group

group_oi = "PreE"

prioritized_tbl_oi_M = prioritized_tbl_oi_network %>% filter(group == group_oi)

plot_oi = make_sample_lr_prod_activity_plots_Omnipath(
  multinichenet_output$prioritization_tables, 
  prioritized_tbl_oi_M %>% inner_join(lr_network_all)
  )
plot_oi

To summarize: this interpretable bubble plot is an important and helpful plot because: 1) these LR interactions are all in the overall top50 of condition-specific interactions 2) they are a likely interaction inducing one or more other prioritized LR interaction and/or they are regulated by one or more other prioritized LR interactions. Because of this, interactions in this plot may be interesting candidates for follow-up experimental validation.

Note: These networks were generated by only looking at the top50 interactions overall. In practice, we encourage users to explore more hits than the top50, certainly if many cell type pairs are considered in the analysis.

All the previous were informative for interactions where both the sender and receiver cell types are captured in the data and where ligand and receptor are sufficiently expressed at the RNA level. However, these two conditions are not always fulfilled and some interesting cell-cell communication signals may be missed as a consequence. Can we still have an idea about these potentially missed interactions? Yes, we can.

Visualize sender-agnostic ligand activities for each receiver-group combination

In the next type of plot, we plot all the ligand activities (both scaled and absolute activities) of each receiver-condition combination. This can give us some insights in active signaling pathways across conditions. Note that we can thus show top ligands based on ligand activity - irrespective and agnostic of expression in sender. Benefits of this analysis are the possibility to infer the activity of ligands that are expressed by cell types that are not in your single-cell dataset or that are hard to pick up at the RNA level.

The following block of code will show how to visualize the activities for the top5 ligands for each receiver cell type - condition combination:

ligands_oi = multinichenet_output$prioritization_tables$ligand_activities_target_de_tbl %>% 
  inner_join(contrast_tbl) %>% 
  group_by(group, receiver) %>% filter(direction_regulation == "up") %>% 
  distinct(ligand, receiver, group, activity) %>% 
  top_n(5, activity) %>% 
  pull(ligand) %>% unique()

plot_oi = make_ligand_activity_plots(
  multinichenet_output$prioritization_tables, 
  ligands_oi, 
  contrast_tbl,
  widths = NULL)
plot_oi

Note you can replace the automatically determined ligands_oi by any set of ligands that are of interest to you.

With this plot/downstream analysis, we end the overview of visualizations that can help you in finding interesting hypotheses about important differential ligand-receptor interactions in your data. In case you ended up with a shortlist of interactions for further checks and potential experimental validation, we recommend going over the visualizations that are introduced in the next section. They are some additional "sound checks" for your shortlist of interactions. However, we don't recommend generating these plots before having thoroughly analyzed and inspected all the previous visualizations. Only go further now if you understood all the previous steps to avoid getting more overwhelmed.

Deep Dive into the data

Visualization of differential ligand-target links

Even though the interpretable bubble plots already provide a lot of information, they do not visualize the specific target genes downstream of the prioritized interactions. Hereby, we still miss some interesting functional information and we cannot assess whether high activity values may be due to a reasonable number of specific target genes or not. Therefore we will now go over some visualizations to inspect target genes downstream of prioritized ligand-receptor interactions.

Without filtering of target genes based on LR-target expression correlation: Ligand activity - target gene combination plots

In this 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 M_Monocyte_CD16 cells as receiver, and look at the top 10 senderLigand-receiverReceptor pairs with these cells as receiver.

group_oi = "PreE"
receiver_oi = "Macrophages"
prioritized_tbl_oi_M_10 = get_top_n_lr_pairs(
  multinichenet_output$prioritization_tables, 
  10, 
  groups_oi = group_oi, 
  receivers_oi = receiver_oi)

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

One observation we can make here is that several genes upregulated in the M-group are indeed high-confident target genes of IFNG (dark purple - high regulatory potential scores). Most of these genes are also potential target genes of TNF, but some specific genes are present as well.

Whereas this plot just showed the top ligands for a certain receiver-contrast, you can also zoom in on specific ligands of interest. As example, we will look at IFNG and IL15:

group_oi = "PreE"
receiver_oi = "Macrophages"
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

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

In summary, these "Ligand activity - target gene combination plots" show how well ligand-target links are supported by general prior knowledge, but not whether they are likely to be active in the system under study. That's what we will look at now.

With filtering of target genes based on LR-target expression correlation: Ligand activity - target gene - expression correlation combination plot

In the previous plot, 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 (pearson correlation > 0.33 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.

group_oi = "PreE"
receiver_oi = "Macrophages"
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.33)) # replace pearson by spearman if you want to filter on the spearman correlation
lr_target_prior_cor_filtered_down = lr_target_prior_cor_filtered %>% 
  filter(direction_regulation == "down") %>% 
  filter( (rank_of_target < top_n_target) & (pearson < -0.33)) # downregulation -- negative correlation - # replace pearson by spearman if you want to filter on the spearman 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:

prioritized_tbl_oi = get_top_n_lr_pairs(
  multinichenet_output$prioritization_tables, 
  50, 
  groups_oi = group_oi, 
  receivers_oi = receiver_oi)

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

This visualization can help users assess whether ligand-target links that are supported by general prior knowledge, are also potentially active in the system under study: target genes that show across-sample expression correlation with their upstream ligand-receptor pairs may be more likely true target genes than target genes that don’t show this pattern.

Even though this plot indicates the strength of the correlation between ligand-receptor expression and target gene expression, it’s hard to assess the pattern of correlation. To help users evaluate whether high correlation values are not due to artifacts, we provide the following LR-target expression scatter plot visualization for a selected LR pair and their targets:

ligand_oi = "IFNG"
receptor_oi = "IFNGR2"
sender_oi = "CD4 T cells"
receiver_oi = "Macrophages"
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

Visualization of ligand-to-target signaling paths

The next type of "sound check" visualization will visualize potential signaling paths between ligands and target genes of interest. In addition to this visualization, we also get a network table documenting the underlying data source(s) behind each of the links shown in this graph. This analysis can help users to assess the trustworthiness of ligand-target predictions. This is strongly recommended before going into experimental validation of ligand-target links.

This inference of 'prior knowledge' ligand-receptor-to-target signaling paths is done similarly to the workflow described in the nichenetr package https://github.com/saeyslab/nichenetr/blob/master/vignettes/ligand_target_signaling_path.md

First read in the required networks:

organism = "human"
if(organism == "human"){
  sig_network = readRDS(url("https://zenodo.org/record/7074291/files/signaling_network_human_21122021.rds")) %>% 
    mutate(from = make.names(from), to = make.names(to))

  gr_network = readRDS(url("https://zenodo.org/record/7074291/files/gr_network_human_21122021.rds")) %>% 
    mutate(from = make.names(from), to = make.names(to))

  ligand_tf_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_tf_matrix_nsga2r_final.rds"))
  colnames(ligand_tf_matrix) = colnames(ligand_tf_matrix) %>% make.names()
  rownames(ligand_tf_matrix) = rownames(ligand_tf_matrix) %>% make.names()

  weighted_networks = readRDS(url("https://zenodo.org/record/7074291/files/weighted_networks_nsga2r_final.rds"))
  weighted_networks$lr_sig = weighted_networks$lr_sig %>% mutate(from = make.names(from), to = make.names(to))
  weighted_networks$gr = weighted_networks$gr %>% mutate(from = make.names(from), to = make.names(to))

} else if(organism == "mouse"){
  sig_network = readRDS(url("https://zenodo.org/record/7074291/files/signaling_network_mouse_21122021.rds")) %>% 
    mutate(from = make.names(from), to = make.names(to))

  gr_network = readRDS(url("https://zenodo.org/record/7074291/files/gr_network_mouse_21122021.rds")) %>% 
    mutate(from = make.names(from), to = make.names(to))

  ligand_tf_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_tf_matrix_nsga2r_final_mouse.rds"))
  colnames(ligand_tf_matrix) = colnames(ligand_tf_matrix) %>% make.names()
  rownames(ligand_tf_matrix) = rownames(ligand_tf_matrix) %>% make.names()

  weighted_networks = readRDS(url("https://zenodo.org/record/7074291/files/weighted_networks_nsga2r_final_mouse.rds"))
  weighted_networks$lr_sig = weighted_networks$lr_sig %>% mutate(from = make.names(from), to = make.names(to))
  weighted_networks$gr = weighted_networks$gr %>% mutate(from = make.names(from), to = make.names(to))
}

Define which ligand and target genes you want to focus on: An interesting possiblity would be to focus on expression-correlated target genes downstream of IFNG, that are also encoding for prioritized ligands. To get these target genes, we rerun the following code IFNG as regulating ligand:

network$links %>% filter(sender_ligand == "CD4 T cells_IFNG" & direction_regulation == "up" & group == "PreE")

ligand_oi = "IFNG"
receptor_oi = "IFNGR2"
targets_all = c("CD274","CD47","CXCL10","CXCL11","CXCL9","TNFSF13B")

active_signaling_network = nichenetr::get_ligand_signaling_path_with_receptor(
  ligand_tf_matrix = ligand_tf_matrix, 
  ligands_all = ligand_oi, 
  receptors_all = receptor_oi, 
  targets_all = targets_all, 
  weighted_networks = weighted_networks, 
  top_n_regulators = 2
  )

data_source_network = nichenetr::infer_supporting_datasources(
  signaling_graph_list = active_signaling_network,
  lr_network = lr_network %>% dplyr::rename(from = ligand, to = receptor), 
  sig_network = sig_network, 
  gr_network = gr_network
  )

active_signaling_network_min_max = active_signaling_network
active_signaling_network_min_max$sig = active_signaling_network_min_max$sig %>% mutate(weight = ((weight-min(weight))/(max(weight)-min(weight))) + 0.75)
active_signaling_network_min_max$gr = active_signaling_network_min_max$gr %>% mutate(weight = ((weight-min(weight))/(max(weight)-min(weight))) + 0.75)
colors = c("ligand" = "purple", "receptor" = "orange", "target" = "royalblue", "mediator" = "grey60")
#ggraph_signaling_path = suppressWarnings(make_ggraph_signaling_path(active_signaling_network_min_max, colors, ligand_oi, receptor_oi, targets_all))
ggraph_signaling_path = make_ggraph_signaling_path(
  active_signaling_network_min_max, 
  colors, 
  ligand_oi, 
  receptor_oi, 
  targets_all)
ggraph_signaling_path$plot

As mentioned, we can also inspect the network table documenting the underlying data source(s) behind each of the links shown in this graph. This analysis can help users to assess the trustworthiness of ligand-target predictions.

data_source_network %>% head()

Visualize top DE genes for a cell type of interest

Finally, we provide some visualizations to just inspect the DE results that were generated during the MultiNicheNet analysis.

group_oi = "PreE"
receiver_oi = "Macrophages" 
DE_genes = multinichenet_output$ligand_activities_targets_DEgenes$de_genes_df %>% 
  inner_join(contrast_tbl) %>% 
  filter(group == group_oi) %>% 
  arrange(p_val) %>% 
  filter(
    receiver == receiver_oi & 
      logFC > 2 & 
      p_val <= 0.05 &
      contrast == contrast_tbl %>% filter(group == group_oi) %>% pull(contrast)) %>% 
  pull(gene) %>% unique()

p_target = make_DEgene_dotplot_pseudobulk(
  genes_oi = DE_genes, 
  celltype_info = multinichenet_output$celltype_info, 
  prioritization_tables = multinichenet_output$prioritization_tables, 
  celltype_oi = receiver_oi, 
  multinichenet_output$grouping_tbl)
p_target$pseudobulk_plot + ggtitle("DE genes (pseudobulk expression)")
p_target$singlecell_plot + ggtitle("DE genes (single-cell expression)")

Among these DE genes, you may be most interested in ligands or receptors

Ligands:

group_oi = "PreE"
receiver_oi = "Macrophages"
DE_genes = multinichenet_output$ligand_activities_targets_DEgenes$de_genes_df %>% 
  inner_join(contrast_tbl) %>% 
  filter(group == group_oi) %>% 
  arrange(p_val) %>% 
  filter(
    receiver == receiver_oi & 
      logFC > 1 & 
      p_val <= 0.05 &
      contrast == contrast_tbl %>% filter(group == group_oi) %>% pull(contrast)) %>% 
  pull(gene) %>% unique()
DE_genes = DE_genes %>% intersect(lr_network$ligand)
p_target = make_DEgene_dotplot_pseudobulk(
  genes_oi = DE_genes, 
  celltype_info = multinichenet_output$celltype_info, 
  prioritization_tables = multinichenet_output$prioritization_tables, 
  celltype_oi = receiver_oi, 
  multinichenet_output$grouping_tbl)
p_target$pseudobulk_plot + ggtitle("DE ligands (pseudobulk expression)")
p_target$singlecell_plot + ggtitle("DE ligands (single-cell expression)")

Receptors:

group_oi = "PreE"
receiver_oi = "Macrophages"
DE_genes = multinichenet_output$ligand_activities_targets_DEgenes$de_genes_df %>% 
  inner_join(contrast_tbl) %>% 
  filter(group == group_oi) %>% 
  arrange(p_val) %>% 
  filter(
    receiver == receiver_oi & 
      logFC > 1 & 
      p_val <= 0.05 &
      contrast == contrast_tbl %>% filter(group == group_oi) %>% pull(contrast)) %>% 
  pull(gene) %>% unique()
DE_genes = DE_genes %>% intersect(lr_network$receptor)
p_target = make_DEgene_dotplot_pseudobulk(
  genes_oi = DE_genes, 
  celltype_info = multinichenet_output$celltype_info, 
  prioritization_tables = multinichenet_output$prioritization_tables, 
  celltype_oi = receiver_oi, 
  multinichenet_output$grouping_tbl)
p_target$pseudobulk_plot + ggtitle("DE receptors (pseudobulk expression)")
p_target$singlecell_plot + ggtitle("DE receptors (single-cell expression)")

saeyslab/multinichenetr documentation built on Jan. 15, 2025, 7:55 p.m.