Single-Cell Spatial Transcriptomics: Cell Communication Analysis (based on NicheNet)

Author: SeekGene

Time: 27 min

Words: 5.4k words

Updated: 2026-06-12

Reads: 0 times

Spatial-seq Analysis Guide Notebooks Cell Communication Analysis

1. Module Introduction

This module is based on the NicheNet algorithm and is used for cell communication analysis on single-cell spatial transcriptomics data.

NicheNet's core advantage is that it not only infers ligand-receptor interactions, but also links the ligands expressed by sender cells to changes in target gene expression in receiver cells. It utilizes prior knowledge integrating ligand-receptor binding, signal transduction, and gene regulatory networks, enabling the prediction of which ligands are most likely to cause the observed downstream gene expression changes. It achieves the following core analysis goals:

• Identify active ligands: Evaluate the extent to which ligands expressed by signal-sending cells (Sender) can explain the expression changes of a specific gene set in signal-receiving cells (Receiver).
• Infer downstream target genes: Predict the signaling pathways activated within receiver cells after ligand-receptor binding and the target gene network they ultimately regulate.

Note: The input for this pipeline is a Seurat object and defined Sender / Receiver cell pairs (Niches). The analysis results will output ligand activity rankings, ligand-receptor-target gene relationship prioritization tables, as well as visualizations like heatmaps, circos plots, and dot plots.

Brief Algorithm Principle: NicheNet calculates the "Regulatory Potential" from specific ligands to target genes by constructing a comprehensive network containing ligand-receptors, signal transduction proteins, and transcription factors. In real data, it evaluates which ligands' regulatory potential scores best predict the expression of these target genes by comparing the differentially expressed genes (or highly expressed gene sets) of receiver cells under different conditions (or different microenvironments), thereby identifying truly active communication events.

2. Input File Preparation

This module requires the following input files or data:

1) Seurat object (.rds): Contains single-cell spatial transcriptomics data that has undergone preprocessing steps such as normalization, dimensionality reduction, and clustering, with completed cell type annotations.

2) Define Microenvironment (Niche):

• Niche: Microenvironment name.
• Sender: The name of the cell type sending signals in the Niche (must match the annotation in Seurat).
• Receiver: The name of the cell type receiving signals in the Niche (must match the annotation in Seurat).

3) NicheNet Prior Network Model Files: Contains the pre-built ligand-receptor network, ligand-target gene regulatory potential matrix, etc. (e.g., ligand_target_matrix.rds, lr_network.rds, which will be loaded automatically by the pipeline).

3. Data Loading and Preprocessing

library(nichenetr)
library(RColorBrewer)
library(tidyverse)
library(Seurat)
library(DT)
library(ggplot2)
library(base64enc)
library(figpatch)
library(circlize)
library(gridExtra)
library(cowplot)

output

Registered S3 method overwritten by 'pROC':
method from
plot.roc spatstat.explore

Warning message:
“replacing previous import ‘e1071::element’ by ‘ggplot2::element’ when loading ‘nichenetr’”
Warning message:
“package ‘RColorBrewer’ was built under R version 4.3.3”
Warning message:
“package ‘tibble’ was built under R version 4.3.3”
Warning message:
“package ‘tidyr’ was built under R version 4.3.3”
Warning message:
“package ‘dplyr’ was built under R version 4.3.3”
Warning message:
“package ‘stringr’ was built under R version 4.3.3”
Warning message:
“package ‘lubridate’ was built under R version 4.3.3”
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──

✓ dplyr 1.1.4

✓ readr 2.1.5

✓ forcats 1.0.0

✓ stringr 1.5.1

✓ ggplot2 4.0.0

✓ tibble 3.2.1

✓ lubridate 1.9.4

✓ tidyr 1.3.1

✓ purrr 1.0.2

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✗ dplyr::filter() masks stats::filter()
✗ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package () to force all conflicts to become errors
Attaching SeuratObject


Attaching package: ‘DT’


The following object is masked from ‘package:SeuratObject’:

JS


The following object is masked from ‘package:Seurat’:

JS


Warning message:
“package ‘circlize’ was built under R version 4.3.3”
========================================
circlize version 0.4.16
CRAN page: https://cran.r-project.org/package=circlize
Github page: https://github.com/jokergoo/circlize
Documentation: https://jokergoo.github.io/circlize_book/book/

If you use it in published research, please cite:
Gu, Z. circlize implements and enhances circular visualization
in R. Bioinformatics 2014.

This message can be suppressed by:
suppressPackageStartupMessages(library(circlize))
========================================


Warning message:
“package ‘gridExtra’ was built under R version 4.3.3”

Attaching package: ‘gridExtra’


The following object is masked from ‘package:dplyr’:

combine


Warning message:
“package ‘cowplot’ was built under R version 4.3.3”

Attaching package: ‘cowplot’


The following object is masked from ‘package:lubridate’:

stamp

## rds_path: Path to the Seurat object rds file, containing expression matrix and cell annotations
rds_path="./rds/25020230_nao_CS.rds"
## meta_path: Path to the metadata table
meta_path="./rds/meta_micro.tsv"
## col_sam: Column name in metadata to select samples for analysis
col_sam="Sample"
## Samples: Sample names to be analyzed, which exist in the col_sam column
sample_name="25020230_nao_CS_expression"
## species: Species, options are "human", "mouse"
species="mouse"
## col_celltype: Column name in the Seurat object metadata storing cell type information
col_celltype="RNA_snn_res.0.2"

data = readRDS(rds_path)
data_subset = subset(data,cells = rownames(data@meta.data[data@meta.data[[col_sam]] %in% sample_name,]))
meta = read.delim(meta_path)
meta_subset = meta[meta[[col_sam]] %in% sample_name,]
data_subset@meta.data = meta_subset

# Read interacting ligand-receptor database and ligand-target matrix
ligand_target_matrix = readRDS(url("https://zenodo.org/record/3260758/files/ligand_target_matrix.rds"))
lr_network = readRDS(url("https://zenodo.org/record/3260758/files/lr_network.rds"))
lr_network = lr_network %>% mutate(bonafide = ! database %in% c("ppi_prediction","ppi_prediction_go"))
lr_network = lr_network %>% dplyr::rename(ligand = from, receptor = to) %>% distinct(ligand, receptor, bonafide)

if(species == "mouse"){
  lr_network = lr_network %>% mutate(ligand = convert_human_to_mouse_symbols(ligand), receptor = convert_human_to_mouse_symbols(receptor)) %>% drop_na()

  colnames(ligand_target_matrix) = ligand_target_matrix %>% colnames() %>% convert_human_to_mouse_symbols()
  rownames(ligand_target_matrix) = ligand_target_matrix %>% rownames() %>% convert_human_to_mouse_symbols()
  ligand_target_matrix = ligand_target_matrix %>% .[!is.na(rownames(ligand_target_matrix)), !is.na(colnames(ligand_target_matrix))]
}

# Define the cell types sending and receiving signals in the niche
niches = list(
  "niche_1" = list(
    "sender" = c("c0", "c10"),
    "receiver" = c("c11")),
  "niche_2" = list(
    "sender" = c("c1","c13"),
    "receiver" = c("c14"))
  )

# Calculate differentially expressed genes of receptors in receiver cells and ligands in sender cells
data_subset = SetIdent(data_subset, value = data_subset[[col_celltype]])
assay_oi = "RNA" 
DE_sender = calculate_niche_de(seurat_obj = data_subset %>% subset(features = lr_network$ligand %>% unique()), niches = niches, type = "sender", assay_oi = assay_oi)
DE_receiver = calculate_niche_de(seurat_obj = data_subset %>% subset(features = lr_network$receptor %>% unique()), niches = niches, type = "receiver", assay_oi = assay_oi) # only receptors now, later on: DE analysis to find targets
DE_sender = DE_sender %>% mutate(avg_log2FC = ifelse(avg_log2FC == Inf, max(avg_log2FC[is.finite(avg_log2FC)]), ifelse(avg_log2FC == -Inf, min(avg_log2FC[is.finite(avg_log2FC)]), avg_log2FC)))
DE_receiver = DE_receiver %>% mutate(avg_log2FC = ifelse(avg_log2FC == Inf, max(avg_log2FC[is.finite(avg_log2FC)]), ifelse(avg_log2FC == -Inf, min(avg_log2FC[is.finite(avg_log2FC)]), avg_log2FC)))

output

[1] "Calculate Sender DE between: c0 and c1"
[2] "Calculate Sender DE between: c0 and c13"
[1] "Calculate Sender DE between: c10 and c1"
[2] "Calculate Sender DE between: c10 and c13"
[1] "Calculate Sender DE between: c1 and c0"
[2] "Calculate Sender DE between: c1 and c10"
[1] "Calculate Sender DE between: c13 and c0"
[2] "Calculate Sender DE between: c13 and c10"
# A tibble: 1 × 2
receiver receiver_other_niche

1 c11 c14
[1] "Calculate receiver DE between: c11 and c14"
[1] "Calculate receiver DE between: c14 and c11"

# Filter low expression genes: use the expression_pct parameter to filter out genes with low expression proportions in the target cell population
expression_pct = 0.75
DE_sender_processed = process_niche_de(DE_table = DE_sender, niches = niches, expression_pct = expression_pct, type = "sender")
DE_receiver_processed = process_niche_de(DE_table = DE_receiver, niches = niches, expression_pct = expression_pct, type = "receiver")

specificity_score_LR_pairs = "min_lfc"
DE_sender_receiver = combine_sender_receiver_de(DE_sender_processed, DE_receiver_processed, lr_network, specificity_score = specificity_score_LR_pairs)
table(DE_sender_receiver$sender)
table(DE_sender_receiver$receiver)

output

c0 c1 c10 c13
11435 11435 11435 11435




c11 c14
22870 22870

# Calculate differentially expressed genes in receiver cells, obtaining highly expressed genes by filtering with lfc_cutoff
lfc_cutoff = 0.15  
specificity_score_targets = "min_lfc"

DE_receiver_targets = calculate_niche_de_targets(seurat_obj = data_subset, niches = niches, lfc_cutoff = lfc_cutoff, expression_pct = expression_pct, assay_oi = assay_oi) 
DE_receiver_processed_targets = process_receiver_target_de(DE_receiver_targets = DE_receiver_targets, niches = niches, expression_pct = expression_pct, specificity_score = specificity_score_targets)

output

[1] "Calculate receiver DE between: c11 and c14"
[1] "Calculate receiver DE between: c14 and c11"

geneset_niche1 = DE_receiver_processed_targets %>% filter(receiver == niches[[1]]$receiver & target_score >= lfc_cutoff & target_significant == 1 & target_present == 1) %>% pull(target) %>% unique()
geneset_niche2 = DE_receiver_processed_targets %>% filter(receiver == niches[[2]]$receiver & target_score >= lfc_cutoff & target_significant == 1 & target_present == 1) %>% pull(target) %>% unique()
  
geneset_niche1 %>% setdiff(rownames(ligand_target_matrix))  
geneset_niche2 %>% setdiff(rownames(ligand_target_matrix))  

length(geneset_niche1) # targets of receiver in niche1
length(geneset_niche2) # targets of receiver in niche2

1. 'Calm3'
2. 'Ubb'
3. 'Cdr1os'
4. 'Ptms'
5. 'AY036118'
6. 'Gm42418'
7. 'Gm19951'
8. 'mt-Co1'

28

143

# Get the list of significantly highly expressed target genes
geneset_niche_list = list()
background = DE_receiver_processed_targets  %>% pull(target) %>% unique()
for(i in 1:length(names(niches))){
    geneset_niche_list[[names(niches)[i]]] = DE_receiver_processed_targets %>% filter(receiver == niches[[i]]$receiver & target_score >= lfc_cutoff & target_significant == 1 & target_present == 1) %>% pull(target) %>% unique()
}

top_n_target = 250
niche_geneset_list = list()
for(i in 1:length(names(niches))){
    niche_geneset_list[[names(niches)[i]]] = list(
        "receiver" = niches[[i]]$receiver,
        "geneset" = geneset_niche_list[[names(niches)[i]]],
        "background" = background
    )
}

# Calculate ligand gene activity scores through differentially expressed target genes in receiver cells
ligand_activities_targets = get_ligand_activities_targets(niche_geneset_list = niche_geneset_list, ligand_target_matrix = ligand_target_matrix, top_n_target = top_n_target)

output

[1] "Calculate Ligand activities for: c11"


Warning message in evaluate_target_prediction(setting, ligand_target_matrix, ligands_position):
“all target gene probability score predictions have same value”
Warning message in cor(prediction, response):
“the standard deviation is zero”
Warning message in cor(prediction, response, method = "s"):
“the standard deviation is zero”


[1] "Calculate Ligand activities for: c14"


Warning message in evaluate_target_prediction(setting, ligand_target_matrix, ligands_position):
“all target gene probability score predictions have same value”
Warning message in cor(prediction, response):
“the standard deviation is zero”
Warning message in cor(prediction, response, method = "s"):
“the standard deviation is zero”

features_oi = union(lr_network$ligand, lr_network$receptor) %>% union(ligand_activities_targets$target) %>% setdiff(NA)

dotplot = suppressWarnings(Seurat::DotPlot(data_subset %>% subset(idents = niches %>% unlist() %>% unique()), features = features_oi, assay = assay_oi, cols = "RdYlBu") + RotatedAxis())
exprs_tbl = dotplot$data %>% as_tibble()
exprs_tbl = exprs_tbl %>% rename(celltype = id, gene = features.plot, expression = avg.exp, expression_scaled = avg.exp.scaled, fraction = pct.exp) %>%
    mutate(fraction = fraction/100) %>% as_tibble() %>% select(celltype, gene, expression, expression_scaled, fraction) %>% distinct() %>% arrange(gene) %>% mutate(gene = as.character(gene))
  
exprs_tbl_ligand = exprs_tbl %>% filter(gene %in% lr_network$ligand) %>% rename(sender = celltype, ligand = gene, ligand_expression = expression, ligand_expression_scaled = expression_scaled, ligand_fraction = fraction) 
exprs_tbl_receptor = exprs_tbl %>% filter(gene %in% lr_network$receptor) %>% rename(receiver = celltype, receptor = gene, receptor_expression = expression, receptor_expression_scaled = expression_scaled, receptor_fraction = fraction)
exprs_tbl_target = exprs_tbl %>% filter(gene %in% ligand_activities_targets$target) %>% rename(receiver = celltype, target = gene, target_expression = expression, target_expression_scaled = expression_scaled, target_fraction = fraction)

4. Expression Dot Plot of Key Genes Across Cell Populations

Figure Legend: This dot plot displays the expression of key genes involved in cell communication analysis across different cell types in the selected Niches.

• X-axis: Features of interest, including key highly expressed ligands, receptors, and regulated target genes inferred by the NicheNet analysis.
• Y-axis: Different cell types participating in cell communication (Sender and Receiver cells).
• Dot Size: Represents the Percent Expressed of the gene in the corresponding cell subpopulation. A larger dot indicates that more cells in that subpopulation express the gene.
• Dot Color: Represents the Average Expression level of the gene in the corresponding cell subpopulation (Scaled).
• Biological Significance: This plot provides visual validation of whether the communication molecules predicted by NicheNet (such as sender cell ligands, receiver cell receptors, and target genes) have genuine, specific high-expression support in the expected cell subpopulations.

options(repr.plot.height=4, repr.plot.width=16)
suppressWarnings(Seurat::DotPlot(data_subset %>% subset(idents = niches %>% unlist() %>% unique()), features = features_oi[1:30], assay = assay_oi, cols = "RdYlBu") + RotatedAxis())
#dotplot

exprs_tbl_ligand = exprs_tbl_ligand %>%  mutate(scaled_ligand_expression_scaled = scale_quantile_adapted(ligand_expression_scaled)) %>% mutate(ligand_fraction_adapted = ligand_fraction) %>% mutate_cond(ligand_fraction >= expression_pct, ligand_fraction_adapted = expression_pct)  %>% mutate(scaled_ligand_fraction_adapted = scale_quantile_adapted(ligand_fraction_adapted))

exprs_tbl_receptor = exprs_tbl_receptor %>% mutate(scaled_receptor_expression_scaled = scale_quantile_adapted(receptor_expression_scaled))  %>% mutate(receptor_fraction_adapted = receptor_fraction) %>% mutate_cond(receptor_fraction >= expression_pct, receptor_fraction_adapted = expression_pct)  %>% mutate(scaled_receptor_fraction_adapted = scale_quantile_adapted(receptor_fraction_adapted))

exprs_sender_receiver = lr_network %>% 
  inner_join(exprs_tbl_ligand, by = c("ligand")) %>% 
  inner_join(exprs_tbl_receptor, by = c("receptor")) %>% inner_join(DE_sender_receiver %>% distinct(niche, sender, receiver))
  
ligand_scaled_receptor_expression_fraction_df = exprs_sender_receiver %>% group_by(ligand, receiver) %>% mutate(rank_receptor_expression = dense_rank(receptor_expression), rank_receptor_fraction  = dense_rank(receptor_fraction)) %>% mutate(ligand_scaled_receptor_expression_fraction = 0.5*( (rank_receptor_fraction / max(rank_receptor_fraction)) + ((rank_receptor_expression / max(rank_receptor_expression))) ) )  %>% distinct(ligand, receptor, receiver, ligand_scaled_receptor_expression_fraction, bonafide) %>% distinct() %>% ungroup()

output

Warning message in inner_join(., exprs_tbl_ligand, by = c("ligand")):
“Detected an unexpected many-to-many relationship between \`x\` and \`y\`.
ℹ Row 1 of \`x\` matches multiple rows in \`y\`.
ℹ Row 109 of \`y\` matches multiple rows in \`x\`.
ℹ If a many-to-many relationship is expected, set \`relationship =
"many-to-many"\` to silence this warning.”
Warning message in inner_join(., exprs_tbl_receptor, by = c("receptor")):
“Detected an unexpected many-to-many relationship between \`x\` and \`y\`.
ℹ Row 1 of \`x\` matches multiple rows in \`y\`.
ℹ Row 673 of \`y\` matches multiple rows in \`x\`.
ℹ If a many-to-many relationship is expected, set \`relationship =
"many-to-many"\` to silence this warning.”
Joining with \`by = join_by(sender, receiver)\`

prioritizing_weights = c("scaled_ligand_score" = 5,
                         "scaled_ligand_expression_scaled" = 1,
                         "ligand_fraction" = 1,
                         "scaled_ligand_score_spatial" = 2, 
                         "scaled_receptor_score" = 0.5,
                         "scaled_receptor_expression_scaled" = 0.5,
                          "receptor_fraction" = 1, 
                         "ligand_scaled_receptor_expression_fraction" = 1,
                         "scaled_receptor_score_spatial" = 0,
                         "scaled_activity" = 0,
                         "scaled_activity_normalized" = 1,
                         "bona_fide" = 1)

include_spatial_info_sender = FALSE # if not spatial info to include: put this to false # user adaptation required on own dataset
include_spatial_info_receiver = FALSE # if spatial info to include: put this to true # user adaptation required on own dataset
if(include_spatial_info_sender == FALSE & include_spatial_info_receiver == FALSE){
    spatial_info = tibble(celltype_region_oi = NA, celltype_other_region = NA) %>% mutate(niche =  niches %>% names() %>% head(1), celltype_type = "sender")
}

sender_spatial_DE_processed = get_non_spatial_de(niches = niches, spatial_info = spatial_info, type = "sender", lr_network = lr_network)
sender_spatial_DE_processed = sender_spatial_DE_processed %>% mutate(scaled_ligand_score_spatial = scale_quantile_adapted(ligand_score_spatial))  
receiver_spatial_DE_processed = get_non_spatial_de(niches = niches, spatial_info = spatial_info, type = "receiver", lr_network = lr_network)
receiver_spatial_DE_processed = receiver_spatial_DE_processed %>% mutate(scaled_receptor_score_spatial = scale_quantile_adapted(receptor_score_spatial))

Ligand-Receptor and Target Gene Prioritization Tables

get_prioritization_tables is a key summarization step in the NicheNet analysis. It integrates the multi-dimensional evidence calculated in all previous steps and ultimately assigns a comprehensive Prioritization Score to each "ligand-receptor" and "ligand-target gene" pair.

• Input Information (output list):
  • DE_sender_receiver: Differential expression colocalization data of ligands and receptors in the Sender and Receiver.
  • ligand_activities_targets: Inferred ligand activity scores and predicted target genes.
  • exprs_tbl_* series: Absolute expression levels and Fraction of ligands, receptors, and target genes in the corresponding cell populations.
• Weight Configuration (prioritizing_weights): This is a vector defining how much weight is assigned to each dimensional metric during the final scoring (e.g., the proportion of ligand activity weight, Sender expression specificity, Receiver expression specificity).
• Output and Biological Significance:
  • Returns a list containing two core tables: prioritization_tbl_ligand_receptor and prioritization_tbl_ligand_target.
  • These two tables form the data foundation for all subsequent visualizations (e.g., circos plots, dot plots, heatmaps). They help researchers quantitatively screen out the top, most worthwhile key communication links for subsequent experimental validation from hundreds or thousands of possible molecular interactions.

output = list(DE_sender_receiver = DE_sender_receiver, ligand_scaled_receptor_expression_fraction_df = ligand_scaled_receptor_expression_fraction_df, sender_spatial_DE_processed = sender_spatial_DE_processed, receiver_spatial_DE_processed = receiver_spatial_DE_processed,
         ligand_activities_targets = ligand_activities_targets, DE_receiver_processed_targets = DE_receiver_processed_targets, exprs_tbl_ligand = exprs_tbl_ligand,  exprs_tbl_receptor = exprs_tbl_receptor, exprs_tbl_target = exprs_tbl_target)
prioritization_tables = get_prioritization_tables(output, prioritizing_weights)

top_ligand_niche_df = prioritization_tables$prioritization_tbl_ligand_receptor %>% select(niche, sender, receiver, ligand, receptor, prioritization_score) %>% group_by(ligand) %>% top_n(1, prioritization_score) %>% ungroup() %>% select(ligand, receptor, niche) %>% rename(top_niche = niche)
top_ligand_receptor_niche_df = prioritization_tables$prioritization_tbl_ligand_receptor %>% select(niche, sender, receiver, ligand, receptor, prioritization_score) %>% group_by(ligand, receptor) %>% top_n(1, prioritization_score) %>% ungroup() %>% select(ligand, receptor, niche) %>% rename(top_niche = niche)

ligand_prioritized_tbl_oi = prioritization_tables$prioritization_tbl_ligand_receptor %>% select(niche, sender, receiver, ligand, prioritization_score) %>% group_by(ligand, niche) %>% top_n(1, prioritization_score) %>% ungroup() %>% distinct() %>% inner_join(top_ligand_niche_df) %>% filter(niche == top_niche) %>% group_by(niche) %>% top_n(50, prioritization_score) %>% ungroup() # get the top50 ligands per niche
table(ligand_prioritized_tbl_oi$receiver)

output

Joining with \`by = join_by(ligand)\`
Warning message in inner_join(., top_ligand_niche_df):
“Detected an unexpected many-to-many relationship between \`x\` and \`y\`.
ℹ Row 273 of \`x\` matches multiple rows in \`y\`.
ℹ Row 39 of \`y\` matches multiple rows in \`x\`.
ℹ If a many-to-many relationship is expected, set \`relationship =
"many-to-many"\` to silence this warning.”




c11 c14
50 50

5. Ligand-Receptor Expression Difference Heatmap

Figure Legend: Ligand-receptor expression difference heatmap.

• Ligand Panel: Displays the expression changes (logFC) of ligands in different Sender cells after NicheNet prioritization.
• Receptor Panel: Displays the expression changes (logFC) of the corresponding receptors in target Receiver cells.
• Biological Significance: Used to visually assess whether the predicted ligand-receptor interactions are supported by significant expression changes. If a ligand is significantly upregulated in sender cells and its receptor is also significantly upregulated in receiver cells, it indicates that this communication link is highly active and highly credible under the current conditions (e.g., disease vs. health comparison).

niche_receiver = c()
for(i in 1:length(names(niches))){
    niche_receiver1 = paste0(names(niches)[i],"_",niches[[i]]$receiver)
    niche_receiver = c(niche_receiver,niche_receiver1)
}


for(i in 1:length(names(niches))){
    receiver_oi = niches[[i]]$receiver
    filtered_ligands = ligand_prioritized_tbl_oi %>% filter(receiver == receiver_oi) %>% pull(ligand) %>% unique()
    prioritized_tbl_oi = prioritization_tables$prioritization_tbl_ligand_receptor %>% filter(ligand %in% filtered_ligands) %>% select(niche, sender, receiver, ligand,  receptor, ligand_receptor, prioritization_score) %>% distinct() %>% inner_join(top_ligand_receptor_niche_df) %>% group_by(ligand) %>% filter(receiver == receiver_oi) %>% top_n(1, prioritization_score) %>% ungroup() 
    lfc_plot = make_ligand_receptor_lfc_plot(receiver_oi, prioritized_tbl_oi, prioritization_tables$prioritization_tbl_ligand_receptor, plot_legend = FALSE, heights = NULL, widths = NULL)
    options(repr.plot.height=8, repr.plot.width=16)
    print(lfc_plot)
    
}

output

Joining with \`by = join_by(ligand, receptor)\`
Warning message in inner_join(., top_ligand_receptor_niche_df):
“Detected an unexpected many-to-many relationship between \`x\` and \`y\`.
ℹ Row 280 of \`x\` matches multiple rows in \`y\`.
ℹ Row 102 of \`y\` matches multiple rows in \`x\`.
ℹ If a many-to-many relationship is expected, set \`relationship =
"many-to-many"\` to silence this warning.”
Joining with \`by = join_by(ligand_receptor, prioritization_score)\`
Joining with \`by = join_by(ligand)\`
Joining with \`by = join_by(prioritization_score, niche, sender, ligand)\`
Joining with \`by = join_by(ligand_receptor, ligand, receptor)\`
Warning message:
“The \`size\` argument of \`element_rect()\` is deprecated as of ggplot2 3.4.0.
ℹ Please use the \`linewidth\` argument instead.
ℹ The deprecated feature was likely used in the nichenetr package.
Please report the issue at .”
Joining with \`by = join_by(ligand, receptor)\`
Joining with \`by = join_by(ligand_receptor, prioritization_score)\`
Joining with \`by = join_by(ligand)\`
Joining with \`by = join_by(prioritization_score, niche, sender, ligand)\`
Joining with \`by = join_by(ligand_receptor, ligand, receptor)\`

6. Ligand-Receptor Interaction Circos Plot

Figure Legend: This circos plot visually displays the high-priority Ligand-Receptor interaction relationships between Sender cells and Receiver cells.

• Sectors:
  • One part of the sectors represents Ligands secreted by Sender cells, whose colors typically correspond to specific Sender cell types.
  • Another part of the sectors represents Receptors expressed by Receiver cells.
• Links/Chords: Connect specific ligands to their corresponding receptors. The presence of a link indicates that NicheNet infers this ligand-receptor pair plays a key role in the current cell communication.
• Width: Reflects the Prioritization Score or inferred interaction weight of the ligand-receptor pair. A wider link indicates a higher priority and potentially greater biological importance for that communication link.
• Biological Significance: Helps researchers quickly identify which core ligands are secreted by which cell types and act on which receptors of target cells within a specific microenvironment, thereby outlining the backbone of intercellular signal transduction.

circos_lr_list = list()

niche_receiver = c()
for(i in 1:length(names(niches))){
    niche_receiver1 = paste0(names(niches)[i],"_",niches[[i]]$receiver)
    niche_receiver = c(niche_receiver,niche_receiver1)
}
n=1
for(i in 1:length(names(niches))){
    
    circos.clear()
    receiver_oi = niches[[i]]$receiver
    print(receiver_oi)

    filtered_ligands = ligand_prioritized_tbl_oi %>% filter(receiver == receiver_oi) %>% top_n(15, prioritization_score) %>% pull(ligand) %>% unique()

    prioritized_tbl_oi = prioritization_tables$prioritization_tbl_ligand_receptor %>% filter(ligand %in% filtered_ligands) %>% select(niche, sender, receiver, ligand,  receptor, ligand_receptor, prioritization_score) %>% distinct() %>% inner_join(top_ligand_receptor_niche_df) %>% group_by(ligand) %>% filter(receiver == receiver_oi) %>% top_n(2, prioritization_score) %>% ungroup() 

    colors_sender = brewer.pal(n = prioritized_tbl_oi$sender %>% unique() %>% sort() %>% length(), name = 'Spectral') %>% magrittr::set_names(prioritized_tbl_oi$sender %>% unique() %>% sort())
    colors_receiver = c("lavender")  %>% magrittr::set_names(prioritized_tbl_oi$receiver %>% unique() %>% sort())

    circos_output = make_circos_lr(prioritized_tbl_oi, colors_sender, colors_receiver,separate_legend = FALSE)
    circos_lr_list[[paste0(names(niches)[i],"_circos")]] = circos_output
    if(n>1){
        graphics.off()
    }
    n=n+1
    
}

circos_output

7. Ligand-Activity-Target Gene Expression Comprehensive Association Plot

Figure Legend: This is a comprehensive association display plot containing 5 side-by-side sub-panels from left to right, integrating multiple dimensions of information including ligand expression, regulatory activity, and downstream target genes in cell communication:

1. Ligand Expression: Displays the expression of high-priority ligands across different cell types in different Niches.
2. Scaled Ligand Activity: Displays the inferred regulatory activity scores (after Scaled normalization) of ligands in Receiver cells, used for intuitive comparison of relative regulatory strengths among different ligands.
3. Ligand Activity: Displays the true raw regulatory potential scores of ligands in Receiver cells, reflecting the absolute ability of ligands to drive downstream target gene expression changes.
4. Predicted Target Genes: Displays the regulatory relationships between ligands and their predicted downstream target genes. The colors in the matrix reflect the regulatory potential of the ligand for specific target genes.
5. Target Expression: Displays the actual expression levels of these regulated key downstream target genes in Receiver cells (or across different groups).

Biological Significance: This plot provides a complete chain of evidence from "signal source generation (Ligand Expression)" to "receptor receiving signal and generating overall response (Ligand Activity)", then to "specific action targets (Predicted Target Genes)" and "final output molecular phenotype changes (Target Expression)", comprehensively and intuitively demonstrating the reliability of NicheNet's prediction results.

receiver_oi = niches[[1]]$receiver
filtered_ligands = ligand_prioritized_tbl_oi %>% filter(receiver == receiver_oi) %>% pull(ligand) %>% unique()
prioritized_tbl_oi = prioritization_tables$prioritization_tbl_ligand_receptor %>% filter(ligand %in% filtered_ligands) %>% select(niche, sender, receiver, ligand,  receptor, ligand_receptor, prioritization_score) %>% distinct() %>% inner_join(top_ligand_receptor_niche_df) %>% group_by(ligand) %>% filter(receiver == receiver_oi) %>% top_n(2, prioritization_score) %>% ungroup() 
exprs_activity_target_plot = make_ligand_activity_target_exprs_plot(receiver_oi, prioritized_tbl_oi,  prioritization_tables$prioritization_tbl_ligand_receptor,  prioritization_tables$prioritization_tbl_ligand_target, output$exprs_tbl_ligand,  output$exprs_tbl_target, lfc_cutoff, ligand_target_matrix, plot_legend = FALSE, heights = NULL, widths = NULL)

output

Joining with \`by = join_by(ligand, receptor)\`
Warning message in inner_join(., top_ligand_receptor_niche_df):
“Detected an unexpected many-to-many relationship between \`x\` and \`y\`.
ℹ Row 280 of \`x\` matches multiple rows in \`y\`.
ℹ Row 102 of \`y\` matches multiple rows in \`x\`.
ℹ If a many-to-many relationship is expected, set \`relationship =
"many-to-many"\` to silence this warning.”
Joining with \`by = join_by(ligand, ligand_score)\`
Joining with \`by = join_by(ligand)\`
Joining with \`by = join_by(sender)\`
Joining with \`by = join_by(target)\`
Warning message:
“Using \`size\` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use \`linewidth\` instead.
ℹ The deprecated feature was likely used in the nichenetr package.
Please report the issue at .”
Warning message:
“The \`size\` argument of \`element_line()\` is deprecated as of ggplot2 3.4.0.
ℹ Please use the \`linewidth\` argument instead.
ℹ The deprecated feature was likely used in the nichenetr package.
Please report the issue at .”
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

exprs_activity_target_plot$combined_plot

receiver_oi = niches[[2]]$receiver
filtered_ligands = ligand_prioritized_tbl_oi %>% filter(receiver == receiver_oi) %>% pull(ligand) %>% unique()
prioritized_tbl_oi = prioritization_tables$prioritization_tbl_ligand_receptor %>% filter(ligand %in% filtered_ligands) %>% select(niche, sender, receiver, ligand,  receptor, ligand_receptor, prioritization_score) %>% distinct() %>% inner_join(top_ligand_receptor_niche_df) %>% group_by(ligand) %>% filter(receiver == receiver_oi) %>% top_n(2, prioritization_score) %>% ungroup() 
exprs_activity_target_plot = make_ligand_activity_target_exprs_plot(receiver_oi, prioritized_tbl_oi,  prioritization_tables$prioritization_tbl_ligand_receptor,  prioritization_tables$prioritization_tbl_ligand_target, output$exprs_tbl_ligand,  output$exprs_tbl_target, lfc_cutoff, ligand_target_matrix, plot_legend = FALSE, heights = NULL, widths = NULL)

output

Joining with \`by = join_by(ligand, receptor)\`
Joining with \`by = join_by(ligand, ligand_score)\`
Joining with \`by = join_by(ligand)\`
Joining with \`by = join_by(sender)\`
Joining with \`by = join_by(target)\`
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

exprs_activity_target_plot$combined_plot

for(i in 1:length(names(niches))){
    receiver_oi = niches[[i]]$receiver
    print(receiver_oi)

    filtered_ligands = ligand_prioritized_tbl_oi %>% filter(receiver == receiver_oi) %>% top_n(15, prioritization_score) %>% pull(ligand) %>% unique()

    prioritized_tbl_oi = prioritization_tables$prioritization_tbl_ligand_receptor %>% filter(ligand %in% filtered_ligands) %>% select(niche, sender, receiver, ligand,  receptor, ligand_receptor, prioritization_score) %>% distinct() %>% inner_join(top_ligand_receptor_niche_df) %>% group_by(ligand) %>% filter(receiver == receiver_oi) %>% top_n(2, prioritization_score) %>% ungroup() 
    write.table(prioritized_tbl_oi,paste0(names(niches)[i],"_ligand_activities.txt"),sep = "\t",col.names = T,row.names = F,quote = F)
}