In browaeysrobin/nichenetr: Modeling Intercellular Communication by Linking Ligands to Target Genes with NicheNet

knitr::opts_chunk$set(
  collapse = TRUE,
  # comment = "#>",
  warning = FALSE,
  message = FALSE
)

Following the Construction of NicheNet's ligand-target model vignette, we will now demonstrate how to use ligand-receptor reactions from LIANA to build the ligand-target model. LIANA is a framework that combines both resources and computational tools for ligand-receptor cell-cell communication inference (Dimitrov et al., 2022). As the NicheNet prior model is built by integrating ligand-receptor, signaling, and gene regulatory databases, each part can be replaced with external data sources. We will show how the first part, the ligand-receptor database, can be replaced with those from LIANA, and how to run the model afterward.

Important: Since LIANA also offers functions to calculate ligand-receptor interactions of interest, it is also possible to use them to select which ligands are of interest to do the ligand activity analysis. This is explained further in the LIANA vignette.

First, we will install LIANA:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

if (!requireNamespace("remotes", quietly = TRUE))
    install.packages("remotes")

remotes::install_github('saezlab/liana')

Load necessary packages.

library(liana)
library(nichenetr)
library(tidyverse)
library(Seurat)

To check which resources are present in LIANA, we can use the show_resources() function. These are then accessed via select_resource().

show_resources()

Next, we will calculate how much overlap there is between the ligands and receptors in the LIANA and NicheNet databases. If the overlap between LIANA receptors and NicheNet signaling network is too low, the integration will probably not work very well. Furthermore, The decomplexify() function of LIANA is crucial in our case, as we would like to separate receptors into their respective subunits.

# Load signaling networks of NicheNet
lr_network_human <- readRDS(url("https://zenodo.org/record/7074291/files/lr_network_human_21122021.rds"))
lr_network_mouse <- readRDS(url("https://zenodo.org/record/7074291/files/lr_network_mouse_21122021.rds"))

sig_network_human <- readRDS(url("https://zenodo.org/record/7074291/files/signaling_network_human_21122021.rds"))
sig_network_mouse <- readRDS(url("https://zenodo.org/record/7074291/files/signaling_network_mouse_21122021.rds"))

overlap_df <- lapply(show_resources()[-1], function(resource) {
  db <- select_resource(resource)[[1]] %>% decomplexify()

  lr_network <- paste0("lr_network_", ifelse(grepl("Mouse", resource), "mouse", "human"))
  sig_network <- paste0("sig_network_", ifelse(grepl("Mouse", resource), "mouse", "human"))

  data.frame(row.names = resource,
             n_ligands = length(unique(db$source_genesymbol)),
             n_receptors = length(unique(db$target_genesymbol)),
             n_ligands_overlap = length(intersect(db$source_genesymbol, get(lr_network)$from)),
             n_receptors_overlap_lr = length(intersect(db$target_genesymbol, get(lr_network)$to)),
             n_receptors_overlap_sig = length(intersect(db$target_genesymbol, get(sig_network)$from))
  ) %>% mutate(frac_ligands_overlap = n_ligands_overlap / n_ligands,
               frac_receptors_overlap_lr = n_receptors_overlap_lr / n_receptors,
               frac_receptors_overlap_sig = n_receptors_overlap_sig / n_receptors)

}) %>% do.call(rbind, .)

overlap_df

On average, ~90% of the ligands and receptors of LIANA databases are in the NicheNet LR network (frac_ligands_overlap, frac_receptors_overlap_lr), and almost all of the receptors in LIANA databases are present in the NicheNet signaling network (frac_receptors_overlap_sig). When using the "Consensus" database of LIANA, there are ~100 ligands that are not present in NicheNet; in contrast, there are 303 ligands in NicheNet that are not present in the LIANA consensus database.

To build the ligand-target model, we can use a very similar code to the Construction of NicheNet's ligand-target model vignette. Users can choose between replacing the NicheNet LR database entirely with LIANA's (replace_nichenet_lr = TRUE), or just adding the LIANA database as an additional data source, which may contain a lot of redundant information.

# Load liana consensus LR network, and rename columns to be compatible with nichenet function
# Users can choose here whether or not to replace the NicheNet LR network entirely, or just add the LIANA db to it
replace_nichenet_lr <- TRUE
liana_db <- select_resource("Consensus")[[1]] %>% decomplexify()
liana_db <- liana_db %>% rename(from = source_genesymbol, to = target_genesymbol) %>% select(from, to) %>% mutate(source = "liana") %>%
  {if (replace_nichenet_lr) (.) else (bind_rows(lr_network_human, .))}

# Change source weights dataframe (but in this case all source weights are 1)
source_weights <- source_weights_df %>%
  add_row(source = "liana", weight = 1, .before = 1)

# Load the gene regulatory network
gr_network_human <- readRDS(url("https://zenodo.org/record/7074291/files/gr_network_human_21122021.rds"))

# Aggregate the individual data sources in a weighted manner to obtain a weighted integrated signaling network
weighted_networks <- construct_weighted_networks(lr_network = liana_db,
                                                 sig_network = sig_network_human,
                                                 gr_network = gr_network_human,
                                                 source_weights_df = source_weights)

# downweigh the importance of signaling and gene regulatory hubs - use the optimized parameters of this
weighted_networks <- apply_hub_corrections(weighted_networks = weighted_networks,
                                          lr_sig_hub = hyperparameter_list %>% filter(parameter == "lr_sig_hub") %>% pull(avg_weight),
                                          gr_hub = hyperparameter_list %>% filter(parameter == "gr_hub") %>% pull(avg_weight))

# in this example we will calculate target gene regulatory potential scores for TNF and the ligand combination TNF+IL6
ligands <- list("TNF",c("TNF","IL6"))
ligand_target_matrix_liana <- construct_ligand_target_matrix(weighted_networks = weighted_networks, ligands = ligands, algorithm = "PPR",
                                                      damping_factor = hyperparameter_list %>% filter(parameter == "damping_factor") %>% pull(avg_weight),
                                                      ltf_cutoff = hyperparameter_list %>% filter(parameter == "ltf_cutoff") %>% pull(avg_weight))

Running NicheNet on the LIANA ligand-target model

In this section compare the results between using the LIANA LR network and NicheNet one in a typical NicheNet analysis. As this is mouse scRNA-seq data, we will build the model using mouse networks ("MouseConsensus" resource in LIANA). Furthermore, instead of using the same source weights for all data sources as in the previous section, we will use the optimized data source weights. Here, we will use the same data source weight for the LIANA model as the one that was computed for the NicheNet LR network. This may not be the optimum weight, but it requires a lot of runtime to optimize this parameter (see Parameter optimization vignette for more information).

# Typical nichenet analysis
# Load seurat object
seuratObj <- readRDS(url("https://zenodo.org/record/3531889/files/seuratObj.rds"))
seuratObj <- UpdateSeuratObject(seuratObj)
seuratObj <- alias_to_symbol_seurat(seuratObj, "mouse") # convert gene names

# Load LIANA mouse consensus ligand-receptor network
lr_network_liana <- select_resource("MouseConsensus")[[1]] %>% decomplexify()
lr_network_liana <- lr_network_liana %>% rename(from = source_genesymbol, to = target_genesymbol) %>% select(from, to) %>% mutate(source = "liana")

# Define receiver cell type
receiver = "CD8 T"
expressed_genes_receiver = get_expressed_genes(receiver, seuratObj, pct = 0.10)

# Define sender cell types
sender_celltypes <- c("CD4 T","Treg", "Mono", "NK", "B", "DC")

# Get expressed genes in the sender
list_expressed_genes_sender <- sender_celltypes %>% unique() %>% lapply(get_expressed_genes, seuratObj, 0.10) # lapply to get the expressed genes of every sender cell type separately here
expressed_genes_sender <- list_expressed_genes_sender %>% unlist() %>% unique()

# Get gene set of interest through DE analysis
seurat_obj_receiver <- subset(seuratObj, idents = receiver)
condition_oi <- "LCMV"; condition_reference <- "SS"
DE_table_receiver <- FindMarkers(object = seurat_obj_receiver,
                                 ident.1 = condition_oi, ident.2 = condition_reference,
                                 group.by = "aggregate", 
                                 min.pct = 0.10) %>% rownames_to_column("gene")
geneset_oi <- DE_table_receiver %>% filter(p_val_adj <= 0.05 & abs(avg_log2FC) >= 0.25) %>% pull(gene)

# Get potential ligands
ligands <- lr_network_liana %>% pull(from) %>% unique()
receptors <- lr_network_liana %>% pull(to) %>% unique()
expressed_ligands <- intersect(ligands,expressed_genes_sender)
expressed_receptors <- intersect(receptors,expressed_genes_receiver)
potential_ligands <- lr_network_liana %>% filter(from %in% expressed_ligands & to %in% expressed_receptors) %>% pull(from) %>% unique()
potential_ligands

# Constructing the ligand-target matrix
gr_network_mouse <- readRDS(url("https://zenodo.org/records/7074291/file/gr_network_mouse_21122021.rds"))

# Define optimum weight as the one calculated for the NicheNet LR network
optim_weight <- optimized_source_weights_df %>% filter(source == "nichenet_verschueren") %>% pull(avg_weight) # 0.2788104
source_weights <- optimized_source_weights_df %>%
  add_row(source = "liana", avg_weight = optim_weight, .before = 1) %>% rename(weight=avg_weight)

# Aggregate the individual data sources in a weighted manner to obtain a weighted integrated signaling network
weighted_networks_liana <- construct_weighted_networks(lr_network = lr_network_liana, sig_network = sig_network_mouse, gr_network = gr_network_mouse, source_weights_df = source_weights)

# Downweigh the importance of signaling and gene regulatory hubs - use the optimized parameters of this
weighted_networks_liana <- apply_hub_corrections(weighted_networks = weighted_networks_liana,
                                           lr_sig_hub = hyperparameter_list %>% filter(parameter == "lr_sig_hub") %>% pull(avg_weight),
                                           gr_hub = hyperparameter_list %>% filter(parameter == "gr_hub") %>% pull(avg_weight))

# Construct ligand-target matrix using only the potential ligands to save time
ligand_target_matrix_liana <- construct_ligand_target_matrix(weighted_networks = weighted_networks_liana, ligands = as.list(potential_ligands), algorithm = "PPR",
                                                       damping_factor = hyperparameter_list %>% filter(parameter == "damping_factor") %>% pull(avg_weight),
                                                       ltf_cutoff = hyperparameter_list %>% filter(parameter == "ltf_cutoff") %>% pull(avg_weight))

# Filter background genes and the gene set of interest to only the ones in the ligand-target matrix
background_expressed_genes <- expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix_liana)]
geneset_oi <- geneset_oi %>% .[. %in% rownames(ligand_target_matrix_liana)]

# Perform ligand activity analysis
ligand_activities <- predict_ligand_activities(geneset = geneset_oi, background_expressed_genes = background_expressed_genes,
                                              ligand_target_matrix = ligand_target_matrix_liana, potential_ligands = potential_ligands)
ligand_activities <- ligand_activities %>% arrange(-aupr_corrected) %>% mutate(rank = rank(desc(aupr_corrected)))
ligand_activities

Compare results with NicheNet ligand-target model

Run NicheNet using the wrapper function.

# Use the aggregate function
ligand_target_matrix <- readRDS(url("https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final_mouse.rds"))
weighted_networks <- readRDS(url("https://zenodo.org/record/7074291/files/weighted_networks_nsga2r_final_mouse.rds"))

nichenet_output <- nichenet_seuratobj_aggregate(
  seurat_obj = seuratObj,
  receiver = "CD8 T",
  condition_colname = "aggregate", condition_oi = "LCMV", condition_reference = "SS",
  sender = c("CD4 T","Treg", "Mono", "NK", "B", "DC"),
  ligand_target_matrix = ligand_target_matrix,
  lr_network = lr_network_mouse %>% distinct(from, to),
  weighted_networks = weighted_networks)

nichenet_output$ligand_activities

Compare the results between LIANA and NicheNet. Here we see that half of the top 20 ligands between the two are the same.

intersect(ligand_activities[1:20,]$test_ligand, nichenet_output$ligand_activities[1:20,]$test_ligand)

Below is a comparison of the rankings. Some of the NicheNet top-ranked ligands are not present in the LIANA LR network, such as Ptprc and some H2- ligands. Nonetheless, the LIANA network seems to have made some new links that results in new ligands appearing, such as Lck, Ccl5, and Crlf2.

rankings_df <- bind_rows(ligand_activities %>% select(test_ligand, rank) %>% mutate(db = "LIANA"),
                         nichenet_output$ligand_activities %>% select(test_ligand, rank) %>% mutate(db = "NicheNet"))
rankings_df <- rankings_df %>% group_by(db) %>% mutate(new_rank = 1:n()) %>%
  group_by(db, rank) %>% mutate(ties = n() > 1, ties = factor(ties, levels = c(TRUE, FALSE))) %>%
  ungroup() %>% mutate(x_label = case_when(db == "NicheNet" ~ 0.8, db == "LIANA" ~ 2.2),
                       db = factor(db, levels = c("NicheNet", "LIANA")))

# Top 20
top_n <- 20
p1 <- ggplot(rankings_df %>% filter(rank <= top_n),
       aes(x=db, y=new_rank, label=test_ligand, group = test_ligand, color = ties)) +
  geom_point() +
  geom_line() +
  geom_text(aes(x = x_label), show.legend = FALSE) +
  scale_color_manual(values = c("tomato", "black")) +
  scale_y_reverse() + theme_classic() +
  labs(y = "Ligand rankings", color = "Tied rankings", title = "Top 20 ligands") +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size=10), legend.position = "none")

# All ligands
p2 <- ggplot(rankings_df, aes(x=db, y=new_rank, label=test_ligand, group = test_ligand, color = ties)) +
  geom_point() +
  geom_line() +
  scale_color_manual(values = c("tomato", "black")) +
  scale_y_reverse() + theme_classic() +
  labs(y = "Ligand rankings", color = "Tied rankings", title = "All ligands") +
  theme(axis.title.x = element_blank(), axis.text.x = element_text(size=10))


cowplot::plot_grid(p1, p2)

References

Dimitrov, D., Türei, D., Garrido-Rodriguez M., Burmedi P.L., Nagai, J.S., Boys, C., Flores, R.O.R., Kim, H., Szalai, B., Costa, I.G., Valdeolivas, A., Dugourd, A. and Saez-Rodriguez, J. Comparison of methods and resources for cell-cell communication inference from single-cell RNA-Seq data. Nat Commun 13, 3224 (2022). https://doi.org/10.1038/s41467-022-30755-0

browaeysrobin/nichenetr documentation built on Sept. 9, 2024, 11:25 a.m.