单细胞空间转录组：细胞通讯分析（基于 NicheNet）

作者: SeekGene

时长: 29 分钟

字数: 6.4k 字

更新: 2026-06-12

阅读: 0 次

空间转录组 分析指南 Notebooks 细胞通讯分析

1. 模块简介

本模块基于 NicheNet 算法，用于对单细胞空间转录组数据进行细胞通讯分析。

NicheNet 的核心优势在于它不仅能推断配体-受体互作，还能将发送细胞表达的配体与接收细胞中靶基因表达的变化联系起来。它利用先验知识整合了配体-受体结合、信号传导和基因调控网络，从而能够预测哪些配体最有可能引起了观察到的下游基因表达变化。实现以下核心分析目标：

• 识别活跃配体：评估信号发送细胞 (Sender) 表达的配体在多大程度上能够解释信号接收细胞 (Receiver) 中特定基因集的表达变化。
• 推断下游靶基因：预测配体结合受体后，在接收细胞内部激活的信号通路及其最终调控的靶基因网络。

注意： 本流程的输入为 Seurat 对象及定义的 Sender / Receiver 细胞对 (Niches)。分析结果将输出配体活性排名、配体-受体-靶基因关系优先级表格，以及热图、和弦图、气泡图等可视化结果。

算法原理简述： NicheNet 通过构建一个包含配体-受体、信号转导蛋白和转录因子的综合网络，计算出从特定配体到目标基因的“调控潜力 (Regulatory Potential)”。在实际数据中，它通过比较接收细胞在不同条件（或不同微环境）下的差异表达基因（或高表达基因集），评估哪些配体的调控潜力得分能够最好地预测这些目标基因的表达，从而识别出真正活跃的通讯事件。

2. 输入文件准备

本模块需要提供以下输入文件或数据：

1) Seurat 对象 (.rds)：包含经过标准化、降维和聚类等预处理步骤的单细胞空间转录组数据，且已完成细胞类型注释。

2) 定义微环境 (Niche)：

• Niche：微环境名称。
• Sender：分析 Niche 中的发送信号的细胞类型名称（需与 Seurat 中的注释一致）。
• Receiver：分析 Niche 中的接收信号的细胞类型名称（需与 Seurat 中的注释一致）。

3) NicheNet 先验网络模型文件：包含预先构建的配体-受体网络、配体-靶基因调控潜力矩阵等（如 ligand_target_matrix.rds、lr_network.rds，流程将自动加载）。

3. 数据加载与预处理

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：Seurat 对象的 rds 文件路径，包含表达矩阵及细胞注释
rds_path="./rds/25020230_nao_CS.rds"
## meta_path：meta表的路径
meta_path="./rds/meta_micro.tsv"
## col_sam：选取meta中分析样本的列名
col_sam="Sample"
## Samples：选取分析的样本名，存在于参数col_sam列下的样本
sample_name="25020230_nao_CS_expression"
## species：物种，可选"human"，"mouse"
species="mouse"
## col_celltype：Seurat 对象 metadata 中存储细胞类型信息的列名
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

#读取互作受配体库以及靶向配体矩阵
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))]
}

#定义niche中发送信号以及接受信号的细胞类型
niches = list(
  "niche_1" = list(
    "sender" = c("c0", "c10"),
    "receiver" = c("c11")),
  "niche_2" = list(
    "sender" = c("c1","c13"),
    "receiver" = c("c14"))
  )

#计算受体与受体细胞，配体与配体细胞的差异表达基因
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"

# 过滤低表达基因 ：利用 expression_pct 参数，过滤掉那些在目标细胞群中表达细胞比例过低的基因
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

# 计算接收细胞的差异表达基因，通过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) #niche1中受体的target
length(geneset_niche2) #niche2中受体的target

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

28

143

#获得显著高表达target基因列表
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
    )
}

#通过受体细胞差异表达target基因，计算配体基因活性得分
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. 关键基因在各细胞群中的表达气泡图

图注说明：该气泡图展示了细胞通讯分析中涉及的关键基因在所选生态位 (Niches) 不同细胞类型中的表达情况。

• X轴（横坐标）：感兴趣的基因 (Features)，包含了 NicheNet 分析推断出的关键高表达配体、受体以及被调控的靶基因。
• Y轴（纵坐标）：参与细胞通讯的不同细胞类型（发送细胞与接收细胞）。
• 气泡大小 (Dot Size)：代表该基因在对应细胞亚群中表达的细胞比例 (Percent Expressed)。气泡越大，表示该亚群中有越多的细胞表达该基因。
• 气泡颜色 (Dot Color)：代表该基因在对应细胞亚群中的平均表达水平（Average Expression，已进行标准化缩放 Scaled）。
• 生物学意义：通过该图可以直观地验证 NicheNet 预测的通讯分子（如发送细胞的配体、接收细胞的受体和靶基因）是否在预期的细胞亚群中存在真实的、特异性的高表达支持。

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))

配体-受体及靶基因优先级表格 (Prioritization Tables)

get_prioritization_tables 是 NicheNet 分析的一个关键汇总步骤。它将前面所有步骤计算出的多维度证据融合在一起，最终为每一对“配体-受体”和“配体-靶基因”打出一个综合的优先级得分 (Prioritization Score)。

• 输入信息 (output 列表)：
  • DE_sender_receiver：Sender 和 Receiver 中配体与受体的差异表达共定位数据。
  • ligand_activities_targets：推断出的配体活性得分及预测的靶基因。
  • exprs_tbl_* 系列：配体、受体、靶基因在对应细胞群体中的绝对表达量和表达比例 (Fraction)。
• 权重配置 (prioritizing_weights)：这是一个向量，定义了在最终打分时，赋予各个维度指标多大的权重（例如：配体活性权重占多少，Sender 表达特异性占多少，Receiver 表达特异性占多少）。
• 输出与生物学意义：
  • 返回一个列表，包含两个核心表格：prioritization_tbl_ligand_receptor（配体-受体优先级表）和 prioritization_tbl_ligand_target（配体-靶基因优先级表）。
  • 这两张表是后续所有可视化绘图（如和弦图、气泡图、热图）的数据基础。它帮研究者从成百上千可能的分子交互中，定量地筛选出最 Top、最值得进行后续实验验证的关键通讯链路。

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) 面板：展示了经过 NicheNet 优先级排序后的配体在不同发送细胞 (Sender cells) 中的表达量变化 (logFC)。
• 受体 (Receptor) 面板：展示了与这些配体相对应的受体在目标接收细胞 (Receiver cells) 中的表达量变化 (logFC)。
• 生物学意义：用于直观评估预测的配体-受体相互作用是否有显著的表达量变化支持。如果一个配体在发送细胞中显著上调，且其受体在接收细胞中也显著上调，说明该通讯连接在当前条件下（如疾病与健康对比）非常活跃且具有较高的可信度。

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. 配体-受体互作和弦图

图注说明：该和弦图直观展示了发送细胞 (Sender) 与接收细胞 (Receiver) 之间高优先级的配体-受体 (Ligand-Receptor) 相互作用关系。

• 圆环的区块 (Sectors)：
  • 一部分区块代表 Sender 细胞分泌的配体 (Ligands)，其颜色通常与特定的 Sender 细胞类型对应。
  • 另一部分区块代表 Receiver 细胞表达的受体 (Receptors)。
• 内部连线 (Links/Chords)：连接特定的配体和其对应的受体。连线的存在代表 NicheNet 推断出该配体-受体对在当前的细胞通讯中起关键作用。
• 连线的宽度 (Width)：反映了该配体-受体对的优先级得分 (Prioritization Score) 或推断的互作权重。连线越宽，表示该通讯链路的优先级越高、生物学上可能越重要。
• 生物学意义：帮助研究者快速识别特定微环境中，哪种细胞类型通过分泌哪些核心配体，主要作用于目标细胞的哪些受体，从而勾勒出细胞间信号传导的骨架。

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. 配体-活性-靶基因表达综合关联图

图注说明：这是一个综合性的关联展示图，从左至右包含 5 个并排的子面板，整合了细胞通讯中配体表达、调控活性以及下游靶基因的多个维度信息：

1. 配体表达 (Ligand Expression)：展示了高优先级配体在不同 Niche（生态位）中不同细胞类型的表达情况。
2. 标准化配体活性 (Scaled Ligand Activity)：展示了配体在 Receiver（接收）细胞中推测的调控活性得分（经过 Scaled 标准化处理），用于直观比较不同配体间相对调控能力的强弱。
3. 配体活性 (Ligand Activity)：展示了配体在 Receiver 细胞中真实的调控潜力原始得分，反映配体驱动下游靶基因表达变化的绝对能力。
4. 预测的靶基因 (Predicted Target Genes)：展示了配体与其预测的下游目标基因之间的调控关系。矩阵中的颜色反映了该配体对特定靶基因的调控潜力。
5. 靶基因表达 (Target Expression)：展示了这些被调控的下游关键靶基因在 Receiver 细胞（或不同组别）中的实际表达水平。

生物学意义：这张图提供了一个从“源头产生信号（配体表达）”到“受体接收信号产生整体响应（配体活性）”，再到“具体作用靶点（预测靶基因）”和“最终输出分子表型变化（靶基因表达）”的完整证据链，全面且直观地论证了 NicheNet 预测结果的可靠性。

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)
}