Robin Browaeys & Chananchida Sang-aram 2023-10-02
In this vignette, we will extend the basic NicheNet analysis analysis from Perform NicheNet analysis starting from a Seurat object: step-by-step analysis by incorporating gene expression as part of the prioritization This is a generalization of the Differential NicheNet and MultiNicheNet approach. While the original NicheNet only ranks ligands based on the ligand activity analysis, it is now also possible to prioritize ligands based on cell type and condition specificity of the ligand and receptor.
We will again make use of mouse NICHE-seq data to explore intercellular communication in the T cell area in the inguinal lymph node before and 72 hours after lymphocytic choriomeningitis virus (LCMV) infection (Medaglia et al. 2017). The used ligand-target matrix and the Seurat object of the processed NICHE-seq single-cell data can be downloaded from Zenodo.
Make sure you understand the different steps in a NicheNet analysis that are described in the basic vignette before proceeding with this vignette.
Load required packages, read in the Seurat object with processed expression data of interacting cells and NicheNet’s ligand-target prior model, ligand-receptor network and weighted integrated networks.
library(nichenetr) # Please update to v2.0.6
library(Seurat)
library(SeuratObject)
library(tidyverse)# Read Seurat object
seuratObj <- readRDS(url("https://zenodo.org/record/3531889/files/seuratObj.rds"))
seuratObj <- UpdateSeuratObject(seuratObj)
seuratObj <- alias_to_symbol_seurat(seuratObj, "mouse")
# Load in networks
lr_network <- readRDS(url("https://zenodo.org/record/7074291/files/lr_network_mouse_21122021.rds"))
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"))
lr_network <- lr_network %>% distinct(from, to)We will use the sender-focused approach here.
# 1. Define set of potential ligands
receiver <- "CD8 T"
expressed_genes_receiver <- get_expressed_genes(receiver, seuratObj, pct = 0.05)
sender_celltypes <- c("CD4 T", "Treg", "Mono", "NK", "B", "DC")
list_expressed_genes_sender <- sender_celltypes %>% unique() %>% lapply(get_expressed_genes, seuratObj, 0.05)
expressed_genes_sender <- list_expressed_genes_sender %>% unlist() %>% unique()
all_ligands <- unique(lr_network$from)
all_receptors <- unique(lr_network$to)
expressed_ligands <- intersect(all_ligands, expressed_genes_sender)
expressed_receptors <- intersect(all_receptors, expressed_genes_receiver)
potential_ligands <- lr_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors) %>%
pull(from) %>% unique()
# 2. Define the gene set of interest
condition_oi <- "LCMV"
condition_reference <- "SS"
seurat_obj_receiver <- subset(seuratObj, idents = receiver)
DE_table_receiver <- FindMarkers(object = seurat_obj_receiver,
ident.1 = condition_oi, ident.2 = condition_reference,
group.by = "aggregate",
min.pct = 0.05) %>% rownames_to_column("gene")
geneset_oi <- DE_table_receiver %>% filter(p_val_adj <= 0.05 & abs(avg_log2FC) >= 0.25) %>% pull(gene)
geneset_oi <- geneset_oi %>% .[. %in% rownames(ligand_target_matrix)]
# 3. Define background genes
background_expressed_genes <- expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]
# 4. Perform NicheNet ligand activity analysis
ligand_activities <- predict_ligand_activities(geneset = geneset_oi,
background_expressed_genes = background_expressed_genes,
ligand_target_matrix = ligand_target_matrix,
potential_ligands = potential_ligands)
ligand_activities <- ligand_activities %>% arrange(-aupr_corrected) %>%
mutate(rank = rank(desc(aupr_corrected)))We will prioritize ligand-receptor pairs based on the following criteria (with their corresponding weight names):
- Upregulation of the ligand in a sender cell type compared to other
cell types:
de_ligand - Upregulation of the receptor in a receiver cell type:
de_receptor - Average expression of the ligand in the sender cell type:
exprs_ligand - Average expression of the receptor in the receiver cell type:
exprs_receptor - Condition-specificity of the ligand across all cell types:
ligand_condition_specificity - Condition-specificity of the receptor across all cell types:
receptor_condition_specificity
This means that we will have to calculate:
- Differential expression of the ligand/receptor in a sender/receiver cell type
- The average expression of each ligand/receptor in each sender/receiver cell type
- Differential expression of the ligand/receptor between the two conditions
We provide a wrapper function generate_info_tables that will calculate
all these values for you. This function returns a list with three
dataframes:
sender_receiver_de: differential expression of the ligand and receptor in the sender-receiver cell type pair. These were first calculated separately (i.e., DE of ligand in sender cell type, DE of receptor in receiver cell type based on FindAllMarkers) and then combined based on possible interactions from the lr_network.sender_receiver_info: the average expression of the ligand and receptor in sender-receiver cell type pairslr_condition_de: differential expression of the ligand and receptor between the two conditions across all cell types.
Note that cell type specificity (i.e., the first four conditions) is calculated only in the condition of interest.
The “scenario” argument can be either “case_control” or “one_condition”. In “case_control” scenario, condition specificity is calculated.
lr_network_filtered <- lr_network %>% filter(from %in% expressed_ligands & to %in% expressed_receptors)
info_tables <- generate_info_tables(seuratObj,
celltype_colname = "celltype",
senders_oi = sender_celltypes,
receivers_oi = receiver,
lr_network = lr_network_filtered,
condition_colname = "aggregate",
condition_oi = condition_oi,
condition_reference = condition_reference,
scenario = "case_control")
names(info_tables)
## [1] "sender_receiver_de" "sender_receiver_info" "lr_condition_de"info_tables$sender_receiver_de %>% head()
## sender receiver ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor pct_expressed_sender pct_expressed_receiver
## 1 DC CD8 T H2-M2 Cd8a 11.002412 2.3838066 6.693109 1.017174e-272 1.377355e-268 5.250531e-206 7.109745e-202 0.429 0.659
## 2 DC CD8 T H2-M2 Klrd1 11.002412 0.9199196 5.961166 1.017174e-272 1.377355e-268 6.104465e-17 8.266056e-13 0.429 0.185
## 3 DC CD8 T Ccl22 Dpp4 9.920608 0.2991720 5.109890 1.590801e-296 2.154103e-292 6.628900e-04 1.000000e+00 0.500 0.148
## 4 DC CD8 T Vsig10 Il6st 10.070530 0.1411494 5.105840 2.637179e-194 3.571005e-190 1.470347e-02 1.000000e+00 0.286 0.090
## 5 DC CD8 T Ccl22 Ccr7 9.920608 0.1468652 5.033737 1.590801e-296 2.154103e-292 5.070025e-05 6.865321e-01 0.500 0.320
## 6 DC CD8 T Cxcl16 Cxcr6 8.101436 1.8384579 4.969947 1.138617e-243 1.541801e-239 5.987787e-21 8.108063e-17 0.929 0.089info_tables$sender_receiver_info %>% head()
## # A tibble: 6 × 7
## sender receiver ligand receptor avg_ligand avg_receptor ligand_receptor_prod
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 DC Mono B2m Tap1 216. 8.59 1856.
## 2 DC NK B2m Klrd1 216. 7.43 1607.
## 3 DC B B2m Tap1 216. 7.35 1588.
## 4 DC Treg B2m Tap1 216. 7.18 1552.
## 5 Mono Mono B2m Tap1 158. 8.59 1353.
## 6 DC DC B2m Tap1 216. 5.91 1277.info_tables$lr_condition_de %>% head()
## ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor
## 1 Cxcl11 Dpp4 7.197344 0.7345098 3.965927 0.0001621364 1 1.170731e-06 1.585287e-02
## 2 Sirpb1c Cd47 6.236414 0.7474147 3.491914 0.0006820290 1 8.720485e-23 1.180841e-18
## 3 Cxcl11 Cxcr3 7.197344 -1.1317386 3.032803 0.0001621364 1 1.918372e-06 2.597667e-02
## 4 Ccl22 Dpp4 5.075469 0.7345098 2.904989 0.0863610523 1 1.170731e-06 1.585287e-02
## 5 F13a1 Itga4 5.436884 0.1228459 2.779865 0.0299628836 1 6.837926e-02 1.000000e+00
## 6 Vcan Sell 5.234169 0.3254999 2.779835 0.0423593686 1 7.148719e-07 9.680080e-03Next, we generate the prioritization table. This table contains the
rankings of ligand-receptor pairs based on the different criteria. We
provide two scenarios: case_control and one_condition. In the
“case_control” scenario, all weights are set to 1. If “one_condition”,
the weights are set to 0 for condition specificity and 1 for the
remaining criteria. Users can also provide their own weights using the
prioritizing_weights argument.
prior_table <- generate_prioritization_tables(info_tables$sender_receiver_info,
info_tables$sender_receiver_de,
ligand_activities,
info_tables$lr_condition_de,
scenario = "case_control")
prior_table %>% head
## # A tibble: 6 × 52
## sender receiver ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor pct_expressed_sender pct_expressed_receiver avg_ligand
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NK CD8 T Ptprc Dpp4 0.642 0.299 0.471 2.18e- 7 2.96e- 3 0.000663 1 0.894 0.148 16.6
## 2 Mono CD8 T Ptprc Dpp4 0.474 0.299 0.386 3.52e- 5 4.77e- 1 0.000663 1 0.867 0.148 14.9
## 3 Treg CD8 T Ptprc Dpp4 0.307 0.299 0.303 1.44e- 2 1 e+ 0 0.000663 1 0.685 0.148 13.2
## 4 Mono CD8 T Cxcl10 Dpp4 4.86 0.299 2.58 2.53e-79 3.43e-75 0.000663 1 0.867 0.148 54.8
## 5 B CD8 T Ptprc Dpp4 0.201 0.299 0.250 2.08e- 2 1 e+ 0 0.000663 1 0.669 0.148 12.3
## 6 Mono CD8 T Ebi3 Il6st 4.00 0.141 2.07 9.77e-49 1.32e-44 0.0147 1 0.147 0.09 0.546
## # ℹ 38 more variables: avg_receptor <dbl>, ligand_receptor_prod <dbl>, lfc_pval_ligand <dbl>, p_val_adapted_ligand <dbl>, scaled_lfc_ligand <dbl>, scaled_p_val_ligand <dbl>,
## # scaled_lfc_pval_ligand <dbl>, scaled_p_val_adapted_ligand <dbl>, activity <dbl>, rank <dbl>, activity_zscore <dbl>, scaled_activity <dbl>, lfc_pval_receptor <dbl>,
## # p_val_adapted_receptor <dbl>, scaled_lfc_receptor <dbl>, scaled_p_val_receptor <dbl>, scaled_lfc_pval_receptor <dbl>, scaled_p_val_adapted_receptor <dbl>, scaled_avg_exprs_ligand <dbl>,
## # scaled_avg_exprs_receptor <dbl>, lfc_ligand_group <dbl>, p_val_ligand_group <dbl>, lfc_pval_ligand_group <dbl>, p_val_adapted_ligand_group <dbl>, scaled_lfc_ligand_group <dbl>,
## # scaled_p_val_ligand_group <dbl>, scaled_lfc_pval_ligand_group <dbl>, scaled_p_val_adapted_ligand_group <dbl>, lfc_receptor_group <dbl>, p_val_receptor_group <dbl>,
## # lfc_pval_receptor_group <dbl>, p_val_adapted_receptor_group <dbl>, scaled_lfc_receptor_group <dbl>, scaled_p_val_receptor_group <dbl>, scaled_lfc_pval_receptor_group <dbl>,
## # scaled_p_val_adapted_receptor_group <dbl>, prioritization_score <dbl>, prioritization_rank <dbl>As you can see, the resulting table now show the rankings for
ligand-receptor interactions of a sender-receiver cell type pair,
instead of just the prioritized ligands. Cxcl10 now went up in the
rankings due to both the high expression of its potential receptor Dpp4
and its high celltype specificity (scaled_lfc_ligand). You can also
see this in the visualizations further below.
We included all columns here, but if you just want relevant columns that were used to calculate the ranking:
prior_table %>% select(c('sender', 'receiver', 'ligand', 'receptor', 'scaled_p_val_adapted_ligand', 'scaled_p_val_adapted_receptor', 'scaled_avg_exprs_ligand', 'scaled_avg_exprs_receptor', 'scaled_p_val_adapted_ligand_group', 'scaled_p_val_adapted_receptor_group', 'scaled_activity'))
## # A tibble: 1,212 × 11
## sender receiver ligand receptor scaled_p_val_adapted_ligand scaled_p_val_adapted_re…¹ scaled_avg_exprs_lig…² scaled_avg_exprs_rec…³ scaled_p_val_adapted…⁴ scaled_p_val_adapted…⁵ scaled_activity
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NK CD8 T Ptprc Dpp4 0.871 0.846 1.00 1.00 0.844 0.833 0.660
## 2 Mono CD8 T Ptprc Dpp4 0.841 0.846 0.867 1.00 0.844 0.833 0.660
## 3 Treg CD8 T Ptprc Dpp4 0.754 0.846 0.741 1.00 0.844 0.833 0.660
## 4 Mono CD8 T Cxcl10 Dpp4 0.958 0.846 1.00 1.00 0.926 0.833 0.309
## 5 B CD8 T Ptprc Dpp4 0.747 0.846 0.666 1.00 0.844 0.833 0.660
## 6 Mono CD8 T Ebi3 Il6st 0.937 0.692 1.00 0.635 0.680 0.515 1.00
## 7 Mono CD8 T Cxcl9 Dpp4 0.974 0.846 1.00 1.00 0.779 0.833 0.263
## 8 DC CD8 T Icam1 Il2rg 0.878 0.723 1.00 0.995 0.713 0.985 0.273
## 9 DC CD8 T B2m Tap2 0.862 0.631 1.00 0.781 1 0.864 0.254
## 10 Mono CD8 T Cxcl11 Dpp4 0.972 0.846 1.00 1.00 0.721 0.833 0.273
## # ℹ 1,202 more rows
## # ℹ abbreviated names: ¹scaled_p_val_adapted_receptor, ²scaled_avg_exprs_ligand, ³scaled_avg_exprs_receptor, ⁴scaled_p_val_adapted_ligand_group, ⁵scaled_p_val_adapted_receptor_groupNote that we appended the suffix ’_group’ to columns that refer to
differential expression between conditions, e.g., lfc_ligand_group and
lfc_receptor_group.
generate_info_tables is a wrapper function that calculates all the
information needed for the prioritization. However, in some cases you
may need more flexibility on how these values are calculated (but note
that you can pass extra arguments to generate_info_tables that will
get passed on to FindMarkers, FindAllMarkers, and
AverageExpression). Below, we show how we use helper functions
calculate_de and get_exprs_avg to calculate the DE and get the
average expression used for cell type specificity. process_table_to_ic
transforms these different dataframes so they are compatible with the
generate_prioritization_tables function.
# Only calculate DE for LCMV condition, with genes that are in the ligand-receptor network
DE_table <- FindAllMarkers(subset(seuratObj, subset = aggregate == "LCMV"),
min.pct = 0, logfc.threshold = 0, return.thresh = 1,
features = unique(unlist(lr_network_filtered)))
# Average expression information - only for LCMV condition
expression_info <- get_exprs_avg(seuratObj, "celltype", condition_colname = "aggregate", condition_oi = condition_oi,
features = unique(unlist(lr_network_filtered)))
# Calculate condition specificity - only for datasets with two conditions!
condition_markers <- FindMarkers(object = seuratObj, ident.1 = condition_oi, ident.2 = condition_reference,
group.by = "aggregate", min.pct = 0, logfc.threshold = 0,
features = unique(unlist(lr_network_filtered))) %>% rownames_to_column("gene")
# Combine DE of senders and receivers -> used for prioritization
processed_DE_table <- process_table_to_ic(DE_table, table_type = "celltype_DE", lr_network_filtered,
senders_oi = sender_celltypes, receivers_oi = receiver)
processed_expr_table <- process_table_to_ic(expression_info, table_type = "expression", lr_network_filtered)
processed_condition_markers <- process_table_to_ic(condition_markers, table_type = "group_DE", lr_network_filtered)And here is how you can define custom weights:
prioritizing_weights = c("de_ligand" = 1,
"de_receptor" = 1,
"activity_scaled" = 1,
"exprs_ligand" = 1,
"exprs_receptor" = 1,
"ligand_condition_specificity" = 1,
"receptor_condition_specificity" = 1)prior_table <- generate_prioritization_tables(processed_expr_table,
processed_DE_table,
ligand_activities,
processed_condition_markers,
prioritizing_weights)
prior_table %>% head
## # A tibble: 6 × 52
## sender receiver ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor pct_expressed_sender pct_expressed_receiver avg_ligand
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NK CD8 T Ptprc Dpp4 0.642 0.299 0.471 2.18e- 7 2.96e- 3 0.000663 1 0.894 0.148 16.6
## 2 Mono CD8 T Ptprc Dpp4 0.474 0.299 0.386 3.52e- 5 4.77e- 1 0.000663 1 0.867 0.148 14.9
## 3 Treg CD8 T Ptprc Dpp4 0.307 0.299 0.303 1.44e- 2 1 e+ 0 0.000663 1 0.685 0.148 13.2
## 4 Mono CD8 T Cxcl10 Dpp4 4.86 0.299 2.58 2.53e-79 3.43e-75 0.000663 1 0.867 0.148 54.8
## 5 B CD8 T Ptprc Dpp4 0.201 0.299 0.250 2.08e- 2 1 e+ 0 0.000663 1 0.669 0.148 12.3
## 6 Mono CD8 T Ebi3 Il6st 4.00 0.141 2.07 9.77e-49 1.32e-44 0.0147 1 0.147 0.09 0.546
## # ℹ 38 more variables: avg_receptor <dbl>, ligand_receptor_prod <dbl>, lfc_pval_ligand <dbl>, p_val_adapted_ligand <dbl>, scaled_lfc_ligand <dbl>, scaled_p_val_ligand <dbl>,
## # scaled_lfc_pval_ligand <dbl>, scaled_p_val_adapted_ligand <dbl>, activity <dbl>, rank <dbl>, activity_zscore <dbl>, scaled_activity <dbl>, lfc_pval_receptor <dbl>,
## # p_val_adapted_receptor <dbl>, scaled_lfc_receptor <dbl>, scaled_p_val_receptor <dbl>, scaled_lfc_pval_receptor <dbl>, scaled_p_val_adapted_receptor <dbl>, scaled_avg_exprs_ligand <dbl>,
## # scaled_avg_exprs_receptor <dbl>, lfc_ligand_group <dbl>, p_val_ligand_group <dbl>, lfc_pval_ligand_group <dbl>, p_val_adapted_ligand_group <dbl>, scaled_lfc_ligand_group <dbl>,
## # scaled_p_val_ligand_group <dbl>, scaled_lfc_pval_ligand_group <dbl>, scaled_p_val_adapted_ligand_group <dbl>, lfc_receptor_group <dbl>, p_val_receptor_group <dbl>,
## # lfc_pval_receptor_group <dbl>, p_val_adapted_receptor_group <dbl>, scaled_lfc_receptor_group <dbl>, scaled_p_val_receptor_group <dbl>, scaled_lfc_pval_receptor_group <dbl>,
## # scaled_p_val_adapted_receptor_group <dbl>, prioritization_score <dbl>, prioritization_rank <dbl>As NicheNet is a receiver-based pipeline, to prioritize ligand-receptor pairs across multiple receivers, we need to perform the NicheNet analysis for each receiver separately. Let’s suppose we want to prioritize ligand-receptor pairs across all T cells (CD4, CD8, and Tregs). The CD8 T analysis has already been performed above. We will use the wrapper function to perform a basic NicheNet analysis on the other two:
nichenet_outputs <- lapply(c("CD8 T", "CD4 T", "Treg"), function(receiver_ct){
output <- nichenet_seuratobj_aggregate(receiver = receiver_ct,
seurat_obj = seuratObj,
condition_colname = "aggregate",
condition_oi = condition_oi,
condition_reference = condition_reference,
sender = sender_celltypes,
ligand_target_matrix = ligand_target_matrix,
lr_network = lr_network,
weighted_networks = weighted_networks,
expression_pct = 0.05)
# Add receiver cell type in ligand activity table
output$ligand_activities$receiver <- receiver_ct
return(output)
})
## [1] "The RNA assay will be used for the analysis."
## [1] "Read in and process NicheNet's networks"
## [1] "Define expressed ligands and receptors in receiver and sender cells"
## [1] "Perform DE analysis in receiver cell"
## [1] "Perform NicheNet ligand activity analysis"
## [1] "Infer active target genes of the prioritized ligands"
## [1] "Infer receptors of the prioritized ligands"
## [1] "Perform DE analysis in sender cells"
## [1] "The RNA assay will be used for the analysis."
## [1] "Read in and process NicheNet's networks"
## [1] "Define expressed ligands and receptors in receiver and sender cells"
## [1] "Perform DE analysis in receiver cell"
## [1] "Perform NicheNet ligand activity analysis"
## [1] "Infer active target genes of the prioritized ligands"
## [1] "Infer receptors of the prioritized ligands"
## [1] "Perform DE analysis in sender cells"
## [1] "The RNA assay will be used for the analysis."
## [1] "Read in and process NicheNet's networks"
## [1] "Define expressed ligands and receptors in receiver and sender cells"
## [1] "Perform DE analysis in receiver cell"
## [1] "Perform NicheNet ligand activity analysis"
## [1] "Infer active target genes of the prioritized ligands"
## [1] "Infer receptors of the prioritized ligands"
## [1] "Perform DE analysis in sender cells"To generate the dataframes used for prioritization, we will simply
change the lr_network_filtered argument to only calculate DE and
expression values for ligand-receptor pairs of interest.
# Calculate prioritization criteria for each receiver cell type
info_tables <- lapply(nichenet_outputs, function(output) {
lr_network_filtered <- lr_network %>% select(from, to) %>%
filter(from %in% output$ligand_activities$test_ligand &
to %in% output$background_expressed_genes)
generate_info_tables(seuratObj,
celltype_colname = "celltype",
senders_oi = sender_celltypes,
receivers_oi = unique(output$ligand_activities$receiver),
lr_network_filtered = lr_network_filtered,
condition_colname = "aggregate",
condition_oi = condition_oi,
condition_reference = condition_reference,
scenario = "case_control")
})We can then combine the results from generate_info_tables using
bind_rows, which will concatenate the rows together. Note that for the
average expression table (sender_receiver_info) and condition
specificity (lr_condition_de), we need to remove duplicate rows.
# bind rows of each element of info_tables using pmap
info_tables_combined <- purrr::pmap(info_tables, bind_rows)
ligand_activities_combined <- purrr::map_dfr(nichenet_outputs, "ligand_activities")
prior_table_combined <- generate_prioritization_tables(
sender_receiver_info = info_tables_combined$sender_receiver_info %>% distinct,
sender_receiver_de = info_tables_combined$sender_receiver_de,
ligand_activities = ligand_activities_combined,
lr_condition_de = info_tables_combined$lr_condition_de %>% distinct,
scenario = "case_control")
head(prior_table_combined)
## # A tibble: 6 × 52
## sender receiver ligand receptor lfc_ligand lfc_receptor ligand_receptor_lfc_avg p_val_ligand p_adj_ligand p_val_receptor p_adj_receptor pct_expressed_sender pct_expressed_receiver avg_ligand
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NK CD8 T Ptprc Dpp4 0.642 0.299 0.471 0.000000218 0.00296 6.63e- 4 1 e+ 0 0.894 0.148 16.6
## 2 NK CD4 T Ptprc Cd4 0.642 1.71 1.18 0.000000218 0.00296 2.63e-34 3.56e-30 0.894 0.226 16.6
## 3 Mono CD8 T Ptprc Dpp4 0.474 0.299 0.386 0.0000352 0.477 6.63e- 4 1 e+ 0 0.867 0.148 14.9
## 4 B CD4 T H2-Eb1 Cd4 5.00 1.71 3.36 0 0 2.63e-34 3.56e-30 0.93 0.226 31.0
## 5 Mono CD4 T Ptprc Cd4 0.474 1.71 1.09 0.0000352 0.477 2.63e-34 3.56e-30 0.867 0.226 14.9
## 6 NK CD4 T Ptprc Cd247 0.642 0.599 0.620 0.000000218 0.00296 5.61e- 4 1 e+ 0 0.894 0.309 16.6
## # ℹ 38 more variables: avg_receptor <dbl>, ligand_receptor_prod <dbl>, lfc_pval_ligand <dbl>, p_val_adapted_ligand <dbl>, scaled_lfc_ligand <dbl>, scaled_p_val_ligand <dbl>,
## # scaled_lfc_pval_ligand <dbl>, scaled_p_val_adapted_ligand <dbl>, activity <dbl>, rank <dbl>, activity_zscore <dbl>, scaled_activity <dbl>, lfc_pval_receptor <dbl>,
## # p_val_adapted_receptor <dbl>, scaled_lfc_receptor <dbl>, scaled_p_val_receptor <dbl>, scaled_lfc_pval_receptor <dbl>, scaled_p_val_adapted_receptor <dbl>, scaled_avg_exprs_ligand <dbl>,
## # scaled_avg_exprs_receptor <dbl>, lfc_ligand_group <dbl>, p_val_ligand_group <dbl>, lfc_pval_ligand_group <dbl>, p_val_adapted_ligand_group <dbl>, scaled_lfc_ligand_group <dbl>,
## # scaled_p_val_ligand_group <dbl>, scaled_lfc_pval_ligand_group <dbl>, scaled_p_val_adapted_ligand_group <dbl>, lfc_receptor_group <dbl>, p_val_receptor_group <dbl>,
## # lfc_pval_receptor_group <dbl>, p_val_adapted_receptor_group <dbl>, scaled_lfc_receptor_group <dbl>, scaled_p_val_receptor_group <dbl>, scaled_lfc_pval_receptor_group <dbl>,
## # scaled_p_val_adapted_receptor_group <dbl>, prioritization_score <dbl>, prioritization_rank <dbl>In addition to the usual heatmap visualizations, we provide a function
make_circos_lr to visualize the ligand-receptor pairs in a circos
plot. This was originally written for the (now deprecated) Differential
NicheNet vignettes. The function takes in a prioritization table and a
named vector for the color of senders and receivers. We first specify
the number of top ligand-receptor pairs to show with n.
# Get top n ligand-receptor pairs
prior_table_oi <- prior_table_combined %>% slice_max(prioritization_score, n = 50)
# Define colors for senders and receivers
senders_receivers <- prior_table_oi %>% select(sender, receiver) %>% unlist %>% unique %>% sort
celltype_colors <- RColorBrewer::brewer.pal(length(senders_receivers), name = 'Set3') %>%
magrittr::set_names(senders_receivers)
circos_plot <- make_circos_lr(prior_table_oi,
colors_sender = celltype_colors, colors_receiver = celltype_colors)circos_plot
Furthermore, we provide the function make_mushroom_plot which allows
you to display expression of ligand-receptor pairs in a specific
receiver. By default, the fill gradient shows the LFC between cell
types, while the size of the semicircle corresponds to the scaled mean
expression. You can also choose to show the rankings of each
ligand-receptor-sender pair with show_rankings, as well as show all
data points for context (show_all_datapoints).
true_color_range = TRUE will adjust the limits of the color gradient
to the min-max of the values, instead of the limit being from 0 to 1.
Note that the numbers displayed here are the rankings across all
receiver cell types (in case of multiple receivers), and by default the
top_n ligand-receptor pairs are shown despite the absolute ranking. To
show only pairs that have an absolute ranking within top_n across all
receivers, set use_absolute_rank = TRUE.
receiver_oi <- "CD8 T"
legend_adjust <- c(0.7, 0.7)
make_mushroom_plot(prior_table_combined %>% filter(receiver == receiver_oi),
top_n = 30,
true_color_range = TRUE,
show_rankings = TRUE,
show_all_datapoints = TRUE) +
theme(legend.justification = legend_adjust,
axis.title.x = element_text(hjust = 0.25))
Furthermore, you can change the “size” and “fill” values to certain
columns from the prioritization table (those with the _ligand or
_receptor suffix).
print(paste0("Column names that you can use are: ", paste0(prior_table %>% select(ends_with(c("_ligand", "_receptor", "_sender", "_receiver"))) %>% colnames() %>%
str_remove("_ligand|_receptor|_sender|_receiver") %>% unique, collapse = ", ")))
## [1] "Column names that you can use are: lfc, p_val, p_adj, avg, lfc_pval, p_val_adapted, scaled_lfc, scaled_p_val, scaled_lfc_pval, scaled_p_val_adapted, scaled_avg_exprs, pct_expressed"# Change size and color columns
make_mushroom_plot(prior_table, top_n = 30, size = "pct_expressed", color = "scaled_avg_exprs") +
theme(legend.justification = legend_adjust,
axis.title.x = element_text(hjust = 0.25))
sessionInfo()
## R version 4.3.3 (2024-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Stream 8
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so; LAPACK version 3.9.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8
## [8] LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Brussels
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] Seurat_5.1.0 SeuratObject_5.0.2 sp_2.1-4 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2 readr_2.1.5 tidyr_1.3.1
## [11] tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0 nichenetr_2.2.0
##
## loaded via a namespace (and not attached):
## [1] RcppAnnoy_0.0.22 splines_4.3.3 later_1.3.2 bitops_1.0-7 polyclip_1.10-6 hardhat_1.3.1 pROC_1.18.5 rpart_4.1.23
## [9] fastDummies_1.7.3 lifecycle_1.0.4 rstatix_0.7.2 doParallel_1.0.17 globals_0.16.3 lattice_0.22-5 MASS_7.3-60.0.1 backports_1.4.1
## [17] magrittr_2.0.3 limma_3.58.1 rmarkdown_2.27 Hmisc_5.1-2 plotly_4.10.4 yaml_2.3.8 httpuv_1.6.15 sctransform_0.4.1
## [25] spam_2.10-0 spatstat.sparse_3.0-3 reticulate_1.37.0 cowplot_1.1.3 pbapply_1.7-2 RColorBrewer_1.1-3 abind_1.4-5 Rtsne_0.17
## [33] presto_1.0.0 BiocGenerics_0.48.1 nnet_7.3-19 tweenr_2.0.3 ipred_0.9-14 circlize_0.4.16 lava_1.8.0 IRanges_2.36.0
## [41] S4Vectors_0.40.2 ggrepel_0.9.5 irlba_2.3.5.1 listenv_0.9.1 spatstat.utils_3.0-4 goftest_1.2-3 RSpectra_0.16-1 spatstat.random_3.2-3
## [49] fitdistrplus_1.1-11 parallelly_1.37.1 leiden_0.4.3.1 codetools_0.2-19 ggforce_0.4.2 tidyselect_1.2.1 shape_1.4.6.1 farver_2.1.2
## [57] matrixStats_1.3.0 stats4_4.3.3 base64enc_0.1-3 spatstat.explore_3.2-7 jsonlite_1.8.8 caret_6.0-94 GetoptLong_1.0.5 e1071_1.7-14
## [65] progressr_0.14.0 Formula_1.2-5 ggridges_0.5.6 survival_3.5-8 iterators_1.0.14 foreach_1.5.2 tools_4.3.3 ggnewscale_0.4.10
## [73] ica_1.0-3 Rcpp_1.0.12 glue_1.7.0 prodlim_2023.08.28 gridExtra_2.3 xfun_0.44 withr_3.0.0 fastmap_1.2.0
## [81] fansi_1.0.6 caTools_1.18.2 digest_0.6.35 gridGraphics_0.5-1 timechange_0.3.0 R6_2.5.1 mime_0.12 colorspace_2.1-0
## [89] scattermore_1.2 tensor_1.5 spatstat.data_3.0-4 DiagrammeR_1.0.11 utf8_1.2.4 generics_0.1.3 data.table_1.15.4 recipes_1.0.10
## [97] class_7.3-22 httr_1.4.7 htmlwidgets_1.6.4 uwot_0.2.2 ModelMetrics_1.2.2.2 pkgconfig_2.0.3 gtable_0.3.5 timeDate_4032.109
## [105] ComplexHeatmap_2.18.0 lmtest_0.9-40 shadowtext_0.1.3 htmltools_0.5.8.1 carData_3.0-5 dotCall64_1.1-1 clue_0.3-65 scales_1.3.0
## [113] png_0.1-8 gower_1.0.1 knitr_1.46 rstudioapi_0.16.0 tzdb_0.4.0 reshape2_1.4.4 rjson_0.2.21 checkmate_2.3.1
## [121] visNetwork_2.1.2 nlme_3.1-164 proxy_0.4-27 zoo_1.8-12 GlobalOptions_0.1.2 KernSmooth_2.23-22 parallel_4.3.3 miniUI_0.1.1.1
## [129] foreign_0.8-86 pillar_1.9.0 grid_4.3.3 vctrs_0.6.5 RANN_2.6.1 ggpubr_0.6.0 randomForest_4.7-1.1 promises_1.3.0
## [137] car_3.1-2 xtable_1.8-4 cluster_2.1.6 htmlTable_2.4.2 evaluate_0.23 cli_3.6.2 compiler_4.3.3 rlang_1.1.3
## [145] crayon_1.5.2 ggsignif_0.6.4 future.apply_1.11.2 labeling_0.4.3 fdrtool_1.2.17 plyr_1.8.9 stringi_1.8.4 viridisLite_0.4.2
## [153] deldir_2.0-4 munsell_0.5.1 lazyeval_0.2.2 spatstat.geom_3.2-9 Matrix_1.6-5 RcppHNSW_0.6.0 hms_1.1.3 patchwork_1.2.0
## [161] future_1.33.2 statmod_1.5.0 shiny_1.8.1.1 highr_0.10 ROCR_1.0-11 broom_1.0.5 igraph_2.0.3Medaglia, Chiara, Amir Giladi, Liat Stoler-Barak, Marco De Giovanni, Tomer Meir Salame, Adi Biram, Eyal David, et al. 2017. “Spatial Reconstruction of Immune Niches by Combining Photoactivatable Reporters and scRNA-Seq.” Science, December, eaao4277. https://doi.org/10.1126/science.aao4277.