本教程基于Linux中的Rstudio环境演示，计算资源不足的同学可参考：生信分析为什么要使用服务器？足够支持你完成硕博生涯的生信环境 配置一个心仪的工作站(硬件+环境配置) 独享服务器，生信分析不求人 为实验室准备一份生物信息学不动产

二、可视化集锦 三、绪论 1、定义 细胞活动的协调对于多细胞存在至关重要，它取决于整个生物体中各种细胞类型和组织之间的细胞间相互作用。细胞间通讯是一个重要过程，它深刻影响生物体的体内平衡、发育和疾病过程。通常，细胞通讯涉及与分泌的配体和质膜受体的相互作用，但它也包括分泌酶、细胞外基质蛋白、转运蛋白和直接的细胞间接触机制。不同的细胞使用不同的细胞通讯来确保生物发育、稳态和组织修复。

例如，在免疫反应期间，细胞通讯使免疫细胞能够识别和对抗病原体。在生长发育过程中，细胞通讯调节细胞增殖和分化，促进器官和组织的正常发育。当细胞无法正确相互作用或误解分子信息时，通常会表现出疾病。细胞通讯的复杂性已被公认为发育生物学、致癌作用和器官功能障碍分子机制的一部分。探索不同条件下的细胞通讯动态变化，可以更深入地了解不同生物过程的潜在机制，并有助于阐明疾病发生和发展背后的机制。

目前，单细胞RNA测序（scRNA-Seq）已在多个研究领域得到广泛应用，以研究配体-受体动力学在细胞间通讯中的关键作用。单细胞RNA测序等技术使研究人员能够探索多细胞生物体内不同细胞类型之间错综复杂的通讯模式，为细胞通讯机制、细胞功能和细胞群组织提供新的视角。细胞间通讯分析有助于了解细胞之间的相互作用，剖析通讯网络，揭示发育过程中的各种细胞相互作用，探索肿瘤免疫微环境，并确定疾病的潜在治疗靶点。

因此，识别和量化细胞间信号通路已成为不同学科的标准做法。通过配体-受体相互作用（LRI）激活特定的细胞信号通路构成了细胞通讯的基本模式，并与各种退行性过程和疾病有着错综复杂的联系。从历史上看，细胞通讯研究主要局限于在体外进行的实验，涉及一种或两种类型的细胞和一组有限的基因。随着科学技术的进步，单细胞水平的数据采集能够检测低丰度基因，并为细胞通讯研究提供了坚实的基础。近年来，许多研究工作致力于制定基于这些相互作用构建细胞通讯网络的策略。

利用这些技术，许多实验室开发了各种用于细胞通讯研究的算法和软件。2、不同细胞通讯软件的比较 这里我们参考了一篇2024年发表在Nature Reviews Genetics杂志上面的文献，文献的标题为“The diversification of methods for studying cell–cell interactions and communication”。

2.1 互作预测过程 文中列举的核心工具可以利用已知的配体-受体相互作用（LRIs）来推断样本中的细胞间相互作用（CCIs）。这些核心方法通过分析配体和受体的表达水平来阐明细胞间如何进行交流，从而提出新的生物学假说。其预测步骤大致可以分为三步：（1）对基因表达矩阵进行筛选，仅保留配体和受体。（2）将每种基因的表达水平在特定细胞类型的所有单细胞中进行汇总（例如，通过计算单细胞的平均表达水平）。

（3）第三，对于样本中每一对细胞类型，评估每个候选LRI，考虑发送细胞类型的配体表达水平和接收细胞类型的受体表达水平。最后针对每一对细胞类型中的每个LRI分别计算一个交流得分（例如，这可以是发送 - 接收细胞类型对中配体和受体表达水平的乘积、平均值或几何平均值）。为了进一步识别基于假设的细胞间相互作用（例如识别特定细胞类型的长程相互作用），可以进一步进行包括置换检验、参数检验和非参数检验在内的统计分析，以确定显著的相互作用。

因此，计算工具能够识别重要的细胞间相互作用，并生成可进行实验验证的生物学假设。2.2 主要研究成果 下图展示了推断细胞间相互作用的计算工具的系统发育树。用于推断细胞间相互作用（CCI）的计算工具从核心工具的“根”部发展而来。从这个根部出发，方法变得更加专业化以应对特定的场景（中心的灰色箭头）。从中心生长出的主要分支代表工具的主要特征（“核心工具”、“更精细”、“更深入”、“更广泛”、“更局部化”或“其他改进”）。

彩色方框表示每个工具命运或子组的次要特征，共展示了105种工具。核心工具依赖于核心或类似的评分函数，是CCI分析的一般框架。例如我们在本教程给大家带来的CellChat、icellnet以及依据python语法开发的CellPhoneDB工具为代表的核心工具，位于系统发育树的核心工具的分支中，供研究界从转录组学中推断细胞间相互作用。2.2.1 基于规则的工具（如icellnet）纳入了关于CCI行为的假设或先验知识。

它们利用与配体和受体数量相关的原理（例如对配体和受体表达进行阈值处理，或使用表达水平作为描述相互作用模式的连续核心函数的输入）来建模相互作用。基于规则的工具通常能产生一致的结果，这是因为它们依赖于基于基因表达的公式。但是，基于规则的方法可能难以应对更高的CCI复杂性、噪声数据和未考虑的变量。2.2.2 基于数据的工具主要使用统计测试或机器学习（如差异分析、标签置换、回归模型、因子分解方法和深度学习）来解释基因表达。

这些方法即使在潜在机制尚不明确的情况下（例如涉及非配体-受体相互作用（非LRI）的情况），也能揭示大型数据集中的意外相关性和隐藏模式。例如，DIALOGUE、MISTy、scITD采用不同类型的分解方法来提取细胞间通讯的特性，尽管它们需要大量的数据。相比之下，数据驱动的工具在处理相同的数据集时可能会产生不同的输出。这是由于统计测试（如置换）的内在随机性以及机器学习算法（如基于梯度的方法）的初始化所导致的。

2.2.3 混合模型通过结合基于规则和数据驱动的策略来发挥各自的优势。例如，CellChat和CellPhoneDB首先通过基于表达式的公式推断细胞间相互作用以确保一致性，然后使用统计测试提取显著的长程相互作用。2.2.4 其他分支 更精细 分支能够以单细胞分辨率全面获取见解。前面我们提到的一些计算工具主要是在伪总体水平上应用的，即把单个细胞聚集成簇或细胞类型。

这种方法有效地解决了scRNA-seq中每个细胞中测量到的转录本通常较为稀疏的问题。然而，最近的方法现在能够以真正的单细胞分辨率处理这些数据，推断出成对单个细胞之间的相互作用。如SPRUCE和DeepCOLOR将单个细胞的基因表达投影到潜在空间，并利用这些信息推断单细胞对之间的细胞间相互作用。更深入 分支可以窥探细胞内活动。细胞通讯接口涉及多种分子类型，包括离子、小分子化合物、肽和蛋白质。

这些分子在细胞外（例如对于长程细胞间通讯）和细胞内（例如对于信号传导、基因表达调控和配体生物合成）都很重要。然而，基于转录组学的工具通常仅限于推断蛋白质配体，其丰度与基因表达的相关性高于其他类型的配体。要估算涉及非蛋白质配体（如小分子）的细胞间通讯，新一代工具可以分析产生或消耗代谢物配体的酶的表达情况，从而捕捉其通讯潜力。MEBOCOST通过整合代谢物受体和代谢物-转运体相互作用的精选数据库来推断基于代谢物的细胞间通讯。

另外，细胞内信号通路还包括转录因子及其下游靶点，也可以被考虑在内。正如NicheNet则可以实现配体→受体→转录因子→靶基因的预测。更局部化 分支从空间上对细胞进行定位。细胞位置会影响细胞间相互作用以及相关的基因表达。从产生它们的细胞中扩散出去，形成浓度梯度，并在接收细胞中触发不同的信号传导程序。这个分支上的工具更广泛地促进了空间数据的分析。

例如，Giotto和Squidpy帮助用户可视化空间数据，并在细胞或点之间进行分析，如聚类和邻域检测，但它们也提供了用于推断CCI的信息。更广泛的 分支可以考虑多种情况。细胞间相互作用在很大程度上取决于其发生的生物学环境。不同的遗传、细胞、细胞外或时间状态等生物学环境会影响单个细胞和整个组织的行为。这些环境在不同表型条件下存在差异，这使得在单细胞图谱中研究细胞间相互作用变得困难，尤其是在样本数量庞大的情况下。

随着样本中实验变量数量的增加，分析的复杂性呈指数级增长。这时候我们可以借助一些工具进行不同条件之间的差异分析帮助我们快速定位不同条件下的细胞通讯差异。如CellChat可以进行细胞类型间、样本间或不同组间两两比较来寻找差异CCI。3、实验追踪细胞间相互作用 除了基于软件分析的方法推断细胞间的细胞通讯，我们其实还可以借助实验的方法去追踪细胞之间的相互作用。例如荧光原位杂交、免疫染色等方法。

然而，这些方法难以用于验证大规模单细胞RNA测序研究得出的CCI预测，因为它们在给定实验中仅能检测少数细胞和长程相互作用（LRIs）。新一代实验方法也正在涌现，提高了CCI测量的通量。如ProximID、PIC-seq这些方法可以同时测量许多LRIs和细胞对。4、总结 近年来，用于研究细胞间通讯（CCIs）的工具经历了显著的多样化发展。计算工具已从基于核心基因表达的方法演进到结合细胞其他生物学特性的下一代策略。

同样，更先进的实验方法提高了通量能力，能够同时分析多个通讯途径，从而为细胞间通讯提供了更全面且具有生物学意义的见解。计算方法和实验方法具有很强的互补性和协同性，这将扩大其在生物医学和个性化医疗等领域的应用潜力。因此，把握这些新机遇对于进一步深化对细胞间通讯的理解至关重要。

四、CellChat 1、简介 视频教程：请点击 Cellchat 基础分析教程视频链接。

2、安装并加载R包 suppressMessages ( if ( ! require (CellChat))devtools :: install_github ( "sqjin/CellChat" )) suppressMessages ( if ( ! require (ggplot2)) install.packages ( "ggplot2" )) suppressMessages ( if ( ! require (patchwork)) install.packages ( "patchwork" ) ) suppressMessages ( if ( ! require (ggalluvial)) install.packages ( "ggalluvial" )) suppressMessages ( if ( ! require (igraph)) install.packages ( "igraph" )) suppressMessages ( if ( ! require (dplyr)) install.packages ( "dplyr" )) 3、下载、读取、创建cellchat对象 Sys.time # [1] "2025-05-15 11:06:22 CST" # 示例数据：来源于人类皮肤 # 下载地址：https://ndownloader.figshare.com/files/25950872 rm ( list = ls ) gc # used (Mb) gc trigger (Mb) max used (Mb) # Ncells 4378697 233.9 8513810 454.7 5569683 297.5 # Vcells 7313817 55.8 12458379 95.1 10290480 78.6 setwd ( "/data/07_scRNA-seq_cell_communication/" ) load ( "data/data_humanSkin_CellChat.rda" ) # 查看数据类型 class (data_humanSkin) # [1] "list" # 获取表达矩阵 data <- data_humanSkin $ data[ 1 : 4 , 1 : 4 ] # 4 x 4 sparse Matrix of class "dgCMatrix" # S1_AACTCCCAGAGCTGCA S1_CAACCAATCCTCATTA S1_CGCTATCTCCTAGTGA # A1BG 0.5774867 1.121317 . # A1BG-AS1 0.5774867 . . # A2M . . . # A2M-AS1 . . . # S1_ATTTCTGCAGGACGTA # A1BG 0.8604636 # A1BG-AS1 . # A2M . # A2M-AS1 . dim (data) # [1] 17328 7563 # 获取细胞描述信息，行名为细胞barcode名称 meta <- data_humanSkin $ meta head (meta) # patient.id condition labels # S1_AACTCCCAGAGCTGCA Patient1 LS Inflam. FIB # S1_CAACCAATCCTCATTA Patient1 LS FBN1+ FIB # S1_CGCTATCTCCTAGTGA Patient1 LS Inflam. FIB # S1_ATTTCTGCAGGACGTA Patient1 LS Inflam. FIB # S1_TGAGCCGAGCTGGAAC Patient1 LS Inflam. FIB # S1_CAGGTGCAGCCCAACC Patient1 LS Inflam. FIB dim (meta) # [1] 7563 3 # 查看分组名称 unique (meta $ condition) # [1] "LS" "NL" # 按指定的变量提取细胞baocode，这里提取出LS组的细胞baocode cell.use <- rownames (meta)[meta $ condition == "LS" ] head (cell.use) # [1] "S1_AACTCCCAGAGCTGCA" "S1_CAACCAATCCTCATTA" "S1_CGCTATCTCCTAGTGA" # [4] "S1_ATTTCTGCAGGACGTA" "S1_TGAGCCGAGCTGGAAC" "S1_CAGGTGCAGCCCAACC" # 取出对应细胞表达矩阵 data_LS <- data[, cell.use] dim (data_LS) # [1] 17328 5011 # 取出对应细胞的meta信息 meta_LS <- meta[cell.use, ] dim (meta_LS) # [1] 5011 3 # 手动创建meta_LS数据框 # meta_test <- data.frame(patient.id = meta[cell.use, "patient.id"], # condition = meta[cell.use, "condition"], # labels = meta[cell.use,"labels"], # row.names = colnames(data_LS)) # dim(meta_test) cellchat <- createCellChat ( object = data_LS, meta = meta_LS, group.by = "labels" ) # [1] "Create a CellChat object from a data matrix" # Set cell identities for the new CellChat object # The cell groups used for CellChat analysis are APOE+ FIB FBN1+ FIB COL11A1+ FIB Inflam. FIB cDC1 cDC2 LC Inflam. DC TC Inflam. TC CD40LG+ TC NKT cellchat # An object of class CellChat created from a single dataset # 17328 genes. # 5011 cells. # CellChat analysis of single cell RNA-seq data! # 如果上述没有指定meta=meta，可以如下手动加入 # cellchat_test <- createCellChat(object = data) # cellchat_test <- addMeta(cellchat, meta = meta) # 用meta中的注释作为分组依据 # cellchat_test <- setIdent(cellchat, ident.use = "labels") # cellchat_test # 细胞类型种类 levels (cellchat @ idents) # show factor levels of the cell labels # [1] "APOE+ FIB" "FBN1+ FIB" "COL11A1+ FIB" "Inflam. FIB" "cDC1" # [6] "cDC2" "LC" "Inflam. DC" "TC" "Inflam. TC" # [11] "CD40LG+ TC" "NKT" # 每种细胞类型数目 table (cellchat @ idents) # # APOE+ FIB FBN1+ FIB COL11A1+ FIB Inflam. FIB cDC1 cDC2 # 1228 813 181 484 121 294 # LC Inflam. DC TC Inflam. TC CD40LG+ TC NKT # 67 81 765 266 630 81 4、载入数据库并开始计算 # 物种为人，设置数据库为CellChatDB.human，如果是小鼠的则换成CellChatDB.mouse CellChatDB <- CellChatDB.human # 展示 CellChat 数据库（CellChatDB）中的配体 - 受体相互作用数据的分类信息，“Secreted Signaling”（分泌信号）、“ECM - Receptor”（细胞外基质 - 受体）和 “Cell - Cell Contact”（细胞 - 细胞接触） showDatabaseCategory (CellChatDB) # dplyr::glimpse 函数用于快速查看数据框或者数据框类似对象的结构和内容，它会以一种紧凑的格式展示数据的基本信息，包括列名、列的数据类型以及每列的前几个值 dplyr :: glimpse (CellChatDB $ interaction) # Rows: 1,939 # Columns: 11 # $ interaction_name <chr> "TGFB1_TGFBR1_TGFBR2", "TGFB2_TGFBR1_TGFBR2", "TGFB… # $ pathway_name <chr> "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TG… # $ ligand <chr> "TGFB1", "TGFB2", "TGFB3", "TGFB1", "TGFB1", "TGFB2… # $ receptor <chr> "TGFbR1_R2", "TGFbR1_R2", "TGFbR1_R2", "ACVR1B_TGFb… # $ agonist <chr> "TGFb agonist", "TGFb agonist", "TGFb agonist", "TG… # $ antagonist <chr> "TGFb antagonist", "TGFb antagonist", "TGFb antagon… # $ co_A_receptor <chr> "", "", "", "", "", "", "", "", "", "", "", "", "",… # $ co_I_receptor <chr> "TGFb inhibition receptor", "TGFb inhibition recept… # $ evidence <chr> "KEGG: hsa04350", "KEGG: hsa04350", "KEGG: hsa04350… # $ annotation <chr> "Secreted Signaling", "Secreted Signaling", "Secret… # $ interaction_name_2 <chr> "TGFB1 - (TGFBR1+TGFBR2)", "TGFB2 - (TGFBR1+TGFBR2)… # 使用所有的分类用作分析数据库 # CellChatDB.use <- CellChatDB # 取出部分分类用作分析数据库 CellChatDB.use <- subsetDB (CellChatDB, search = "Secreted Signaling" ) # 将数据库内容载入cellchat对象中 cellchat @ DB <- CellChatDB.use # 表达量预处理 # 在 CellChatDB.use 中对信号基因的表达数据进行子集处理 cellchat <- subsetData (cellchat, features = NULL ) # 寻找高表达的基因 cellchat <- identifyOverExpressedGenes (cellchat) # identifyOverExpressedInteractions函数用于寻找高表达的细胞间相互作用（配体 - 受体相互作用）。

细胞间的通讯是通过配体与受体的结合来实现的，高表达的相互作用可能代表着更活跃的细胞通讯通路。函数会基于基因表达数据以及配体 - 受体相互作用信息，找出那些表达水平较高的相互作用，结果同样会存储在cellchat对象中 cellchat <- identifyOverExpressedInteractions (cellchat) # 将cellchat对象中的数据投影到人类蛋白质 - 蛋白质相互作用网络（PPI.human）上。

蛋白质 - 蛋白质相互作用网络能够提供额外的信息，帮助进一步理解细胞间通讯的分子机制。通过投影操作，可以将细胞间通讯分析与蛋白质相互作用网络相结合，挖掘出更多潜在的生物学意义 cellchat <- projectData (cellchat, PPI.human) # 计算细胞间通讯的概率。细胞间通讯并非是绝对的，而是存在一定的概率。

computeCommunProb函数会根据基因表达数据、配体 - 受体相互作用信息等，采用默认的计算方式（这里 type = "truncatedMean"）来计算不同细胞群体之间通讯的概率。raw.use = T 表示使用原始的基因表达数据进行计算。

计算得到的通讯概率会存储在 cellchat 对象中，为后续分析细胞间通讯的强度和模式提供依据 # 默认cutoff的值为20%，即表达比例在25%以下的基因会被认为是0，trim = 0.1可以调整比例阈值 cellchat <- computeCommunProb (cellchat, raw.use = T) # triMean is used for calculating the averagene expression per cell group. # [1] ">>> Run CellChat on sc/snRNA-seq data <<< [2025-05-15 11:10:57.088421]" # [1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2025-05-15 11:14:10.908108]" # 对细胞间通讯结果进行过滤。

min.cells = 10 表示只保留那些涉及至少 10 个细胞的通讯相互作用。在实际分析中，一些细胞间的通讯可能只涉及极少数细胞，这些通讯可能是噪声或者不具有普遍意义。

通过设置最小细胞数的阈值，可以过滤掉这些可能的噪声数据，使得分析结果更加可靠 cellchat <- filterCommunication (cellchat, min.cells = 10 ) # 从 cellchat 对象中提取细胞间通讯的相关信息，并将其存储在一个数据框 df.net 中。这个数据框包含了经过前面一系列分析和过滤后得到的细胞间通讯的关键信息，例如细胞群体之间的通讯概率、涉及的配体 - 受体相互作用等。

通过这个数据框，你可以方便地进行后续的数据可视化、统计分析等操作 df.net <- subsetCommunication (cellchat) # DT 包是一个用于创建交互式表格的 R 包，它基于 JavaScript 库 DataTables。

datatable 函数的主要作用是将数据框转换为一个交互式的 HTML 表格，这个表格具有丰富的交互功能，例如排序、筛选、分页等 DT :: datatable (df.net) write.csv (df.net, "./result_cellchat/01.df.net.csv" ) # computeCommunProbPathway 函数的主要功能是计算细胞间通讯在信号通路层面的概率。

在细胞间通讯分析里，仅仅知道细胞群体之间存在通讯是不够的，还需要了解这些通讯是通过哪些信号通路来实现的，以及每条信号通路的通讯概率是多少。该函数会基于之前计算得到的细胞间通讯概率（例如通过 computeCommunProb 函数计算），进一步分析每条信号通路在不同细胞群体之间的通讯概率。

# 每对配受体的预测结果存在net中，每条通路的预测结果存在netp中 cellchat <- computeCommunProbPathway (cellchat) # 计算联路数与通讯概率，可用sources.use and targets.use指定来源与去向 cellchat <- aggregateNet (cellchat) 5、可视化 # 统计每类细胞的数目 groupSize <- as.numeric ( table (cellchat @ idents)) library (patchwork) # 表示将绘图区域划分为 1 行 3 列的网格布局，并且允许图形元素绘制在绘图区域之外，这样在绘制一些带有标签、注释等可能超出绘图框的图形时，这些超出部分也能正常显示 par ( mfrow = c ( 1 , 3 ), xpd= TRUE ) netVisual_circle (cellchat @ net $ count, vertex.weight = groupSize, # 指定节点（细胞类型）权重，节点大小会根据这个权重进行调整，这里依据细胞数目 weight.scale = T, label.edge= F, title.name = "Number of interactions" ) netVisual_circle (cellchat @ net $ weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength" ) netVisual_circle (cellchat @ net $ weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength" , targets.use = "cDC2" ) cDC2 细胞群体相关的通讯相互作用 mat <- cellchat @ net $ count par ( mfrow = c ( 3 , 4 ), xpd= TRUE ) for (i in 1 : nrow (mat)) { # mat2 <- matrix ( 0 , nrow = nrow (mat), ncol = ncol (mat), dimnames = dimnames (mat)) mat2[i, ] <- mat[i, ] netVisual_circle (mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max (mat), title.name = rownames (mat)[i]) } mat <- cellchat @ net $ weight par ( mfrow = c ( 3 , 4 ), xpd= TRUE ) for (i in 1 : nrow (mat)) { mat2 <- matrix ( 0 , nrow = nrow (mat), ncol = ncol (mat), dimnames = dimnames (mat)) mat2[i, ] <- mat[i, ] netVisual_circle (mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max (mat), title.name = rownames (mat)[i]) } 6、进阶可视化 pathways.show <- df.net $ pathway_name pathways.show[ 1 ] # [1] "FGF" # 层次图（Hierarchy plot） par ( mfrow = c ( 1 , 2 ), xpd= TRUE ) vertex.receiver = seq ( 1 , 4 ) netVisual_aggregate (cellchat, signaling = pathways.show[ 1 ], vertex.receiver = vertex.receiver, layout = "hierarchy" ) netVisual_aggregate (cellchat, signaling = pathways.show[ 1 ], layout = "circle" ) # 采用弦图（chord diagram）的布局方式来展示细胞间的通讯关系 netVisual_aggregate (cellchat, signaling = pathways.show[ 1 ], layout = "chord" ) # 分组展示弦图 group.cellType <- c ( rep ( "FIB" , 4 ), rep ( "DC" , 4 ), rep ( "TC" , 4 )) names (group.cellType) <- levels (cellchat @ idents) group.cellType # APOE+ FIB FBN1+ FIB COL11A1+ FIB Inflam. FIB cDC1 cDC2 # "FIB" "FIB" "FIB" "FIB" "DC" "DC" # LC Inflam. DC TC Inflam. TC CD40LG+ TC NKT # "DC" "DC" "TC" "TC" "TC" "TC" par ( mfrow = c ( 1 , 2 ), xpd= TRUE ) netVisual_chord_cell (cellchat, signaling = pathways.show[ 1 ], group = group.cellType, title.name = paste0 (pathways.show[ 1 ], " signaling network" )) netVisual_chord_cell (cellchat, signaling = pathways.show[ 1 ], title.name = paste0 (pathways.show[ 1 ], " signaling network" )) 7、配受体展示 # 展现对特定通路的贡献程度 p1 <- netAnalysis_contribution (cellchat, signaling = "MIF" , title = "MIF" ) # 计算每对配体-受体对整体信号通路的贡献，并可视化由单一配体-受体对介导的细胞间通信。

p2 <- netAnalysis_contribution (cellchat, signaling = pathways.show) cowplot :: plot_grid (p1, p2, align = "h" , ncol= 2 ) pairLR.CXCL <- extractEnrichedLR (cellchat, signaling = pathways.show[ 1 ], geneLR.return = FALSE ) # 提取第一个显著的配受体对 LR.show <- pairLR.CXCL[ 1 ,] # 对应前四个可以作为受体的细胞类型 # 如果出现Error in plot.new : figure margins too large报错，请把plot窗口拉大一些 vertex.receiver = seq ( 1 , 4 ) p1 <- netVisual_individual (cellchat, signaling = pathways.show, pairLR.use = LR.show, vertex.receiver = vertex.receiver, layout = "hierarchy" ) vertex.receiver = seq ( 1 , 7 ) p2 <- netVisual_individual (cellchat, signaling = pathways.show, pairLR.use = LR.show, vertex.receiver = vertex.receiver, layout = "hierarchy" ) # cowplot::plot_grid(p1, p2 ,align = "h",ncol=1) LR.show # [1] "FGF7_FGFR1" # 单个互作对圈图 netVisual_individual (cellchat, signaling = pathways.show[ 1 ], pairLR.use = LR.show, layout = "circle" ) # [[1]] # 单个互作对弦图 netVisual_individual (cellchat, signaling = pathways.show[ 1 ], pairLR.use = LR.show, layout = "chord" ) # [[1]] 8、通路展示 pathways.show.all <- cellchat @ netP $ pathways.show.all # [1] "MIF" "GALECTIN" "CXCL" "COMPLEMENT" "FGF" # [6] "TNF" "CCL" "GAS" "IL4" "CD40" # [11] "LIGHT" "CSF" "VEGF" for (i in 1 : length (pathways.show.all)) { gg <- netAnalysis_contribution (cellchat, signaling = pathways.show.all[i]) ggsave ( filename= paste0 ( "result_cellchat/" ,pathways.show.all[i], "_L-R_contribution.pdf" ), plot= gg, width = 7 , height = 7 , units = "in" , dpi = 300 ) } netAnalysis_contribution (cellchat, signaling = pathways.show.all) 9、多个细胞通讯通路气泡图可视化 levels (cellchat @ idents) # [1] "APOE+ FIB" "FBN1+ FIB" "COL11A1+ FIB" "Inflam. FIB" "cDC1" # [6] "cDC2" "LC" "Inflam. DC" "TC" "Inflam. TC" # [11] "CD40LG+ TC" "NKT" par ( mfrow = c ( 1 , 4 ), xpd= TRUE ) p1 <- netVisual_bubble (cellchat, sources.use = 4 , targets.use = c ( 5 : 11 ), remove.isolate = FALSE ) + # 将 x 轴的刻度标签旋转 45 度，hjust = 1：标签右对齐（与坐标轴的刻度线对齐），vjust = 1：标签与 x 轴顶部对齐，给 x 轴的刻度标签上方添加 1 个单位的空间（即距离顶部增加） theme ( axis.text.x = element_text ( angle = 45 , hjust = 1 , vjust = 1 , margin = margin ( t = 1 ))) p2 <- netVisual_bubble (cellchat, sources.use = 4 , targets.use = c ( "cDC1" , "cDC2" , "LC" , "Inflam. DC" , "TC" , "Inflam. TC" , "CD40LG+ TC" ), remove.isolate = FALSE ) + theme ( axis.text.x = element_text ( angle = 45 , hjust = 1 , vjust = 1 , margin = margin ( t = 1 ))) p3 <- netVisual_bubble (cellchat, sources.use = 4 , targets.use = c ( "cDC1" , "cDC2" ), remove.isolate = FALSE ) + theme ( axis.text.x = element_text ( angle = 45 , hjust = 1 , vjust = 1 , margin = margin ( t = 1 ))) p4 <- netVisual_bubble (cellchat, sources.use = 4 , targets.use = c ( 5 : 11 ), signaling = c ( "CCL" , "CXCL" ), remove.isolate = FALSE ) + theme ( axis.text.x = element_text ( angle = 45 , hjust = 1 , vjust = 1 , margin = margin ( t = 1 ))) cowplot :: plot_grid (p1, p2,p3,p4, align = "h" , ncol= 4 ) par ( mfrow = c ( 1 , 3 ), xpd= TRUE ) pairLR.use <- extractEnrichedLR (cellchat, signaling = c ( "CCL" , "CXCL" , "MIF" )) pairLR.use # interaction_name # 1 CCL19_CCR7 # 2 CXCL12_CXCR4 # 3 CXCL12_ACKR3 # 4 MIF_CD74_CXCR4 # 5 MIF_CD74_CD44 # 6 MIF_ACKR3 p1 <- netVisual_bubble (cellchat, sources.use = c ( 3 , 4 ), targets.use = c ( 5 : 8 ), pairLR.use = pairLR.use, remove.isolate = TRUE ) + # 将 x 轴的刻度标签旋转 45 度，hjust = 1：标签右对齐（与坐标轴的刻度线对齐），vjust = 1：标签与 x 轴顶部对齐，给 x 轴的刻度标签上方添加 1 个单位的空间（即距离顶部增加） theme ( axis.text.x = element_text ( angle = 45 , hjust = 1 , vjust = 1 , margin = margin ( t = 1 ))) p2 <- netVisual_bubble (cellchat, sources.use = 4 , targets.use = c ( 5 : 11 ), signaling = c ( "CCL" , "CXCL" , "MIF" ), remove.isolate = FALSE ) + theme ( axis.text.x = element_text ( angle = 45 , hjust = 1 , vjust = 1 , margin = margin ( t = 1 ))) # 孤立节点：在细胞通讯网络分析里，孤立节点指的是那些没有与其他节点（细胞类型）产生信号通讯的节点。

也就是该细胞类型既不向其他细胞类型发送信号，也不接收其他细胞类型发出的信号 # remove.isolate 设置成 TRUE 时，netVisual_bubble函数会在绘图之前将所有孤立的节点从网络中移除。

这样做的好处是能够让可视化的结果更加简洁，把重点聚焦在那些存在信号通讯的细胞类型上，避免孤立节点对整体网络结构和信号通讯模式的干扰 # p3 <- netVisual_bubble (cellchat, sources.use = 4 , targets.use = c ( 5 : 11 ), signaling = c ( "CCL" , "CXCL" , "MIF" ), remove.isolate = T) + theme ( axis.text.x = element_text ( angle = 45 , hjust = 1 , vjust = 1 , margin = margin ( t = 1 ))) cowplot :: plot_grid (p1, p2,p3, align = "h" , ncol= 3 ) 10、多组别细胞通讯 视频教程：请点击 CellChat 多组别通讯视频链接 这部分的内容主要分析目的为：1、细胞通讯在不同组别中是否发生变化 2、细胞通讯在不同细胞类型中是否发生变化 3、细胞通讯的“发起者”与接收者是否随组别发生变化 10.1 准备工作 library (CellChat) library (patchwork) library (cowplot) dir.create ( "./comparison" ) setwd ( "./comparison" ) # 下载地址：cellchat.NL <- readRDS ( "data/cellchat_humanSkin_NL.rds" ) cellchat.LS <- readRDS ( "data/cellchat_humanSkin_LS.rds" ) # 如果出现报错：no slot of name "images" for this object of class "CellChat"，是因为cellchat版本不同导致的，可以先使用updateCellChat函数将之前对象更新后再执行mergeCellChat合并的操作 cellchat.NL <- updateCellChat (cellchat.NL) cellchat.LS <- updateCellChat (cellchat.LS) object.list <- list ( NL = cellchat.NL, LS = cellchat.LS) cellchat <- mergeCellChat (object.list, add.names = names (object.list)) cellchat # An object of class CellChat created from a merged object with multiple datasets # 555 signaling genes. # 7563 cells. # CellChat analysis of single cell RNA-seq data! 10.2 可视化 最简单的展示，查看细胞互作的数量在不同条件下是否有差异 # 比较cellchat对象里group1和group2之间的细胞间相互作用 # measure：该参数指定了用于比较相互作用的衡量指标。

默认值为"count"，也就是相互作用的次数 # 在第二行代码中，measure = "weight" 表明使用相互作用的权重来进行比较。gg1 <- compareInteractions (cellchat, show.legend = F, group = c ( 1 , 2 )) # 权重反映了特定细胞群之间信号传递的相对强度。

例如，一个高权重的相互作用可能意味着发送细胞群能够高效地产生并释放信号分子，而接收细胞群对这些信号分子有高度的响应能力。这种相互作用在细胞的生理过程（如细胞增殖、分化、免疫反应等）中可能起着更为关键的作用。

影响因素包括了信号分子的表达水平、受体的表达水平、信号通路的活性 gg2 <- compareInteractions (cellchat, show.legend = F, group = c ( 1 , 2 ), measure = "weight" ) gg1 + gg2 # 查看细胞通路在两组间的富集程度 # 比较每个信号传导途径的信息流来识别保守的和特定于情境的信号传导途径，信息流由推断网络中所有细胞对之间的通信概率之和 (即网络中的总权重) 定义。

# 条形图可以用堆叠模式绘制，也可以不用。重要的信号传导途径根据NL和LS组之间推断网络内整体信息流的差异进行排序。顶部红色的信号传导途径在NL中富集，而这些绿色的信号传导途径在LS中富集。

gg1 <- rankNet (cellchat, mode = "comparison" , stacked = T, do.stat = TRUE ) gg2 <- rankNet (cellchat, mode = "comparison" , stacked = F, do.stat = TRUE ) gg1 + gg2 以 circle plot 的形式展示第二个组别中相较于第一个组别细胞通讯发生的变化，红色为上调蓝色为下调 par ( mfrow = c ( 1 , 2 ), xpd= TRUE ) # 不同细胞群体之间的交互次数或交互强度差异，可以使用圆形图可视化两个数据集之间细胞-细胞通信网络中的交互次数或交互强度差异 netVisual_diffInteraction (cellchat, weight.scale = T) netVisual_diffInteraction (cellchat, weight.scale = T, measure = "weight" ) 以上是直接展示二组“相减”的结果，当然你也可以直接将两组分开展示: # getMaxWeight 是用来提取在 object.list 中每个 CellChat 对象的最大权重值。

这个函数根据指定的 attribute 参数计算最大的权重值。attribute = c("idents", "count") 表示从每个对象中的 idents（即细胞标识符）和 count（即互动计数）属性中提取最大权重值。

weight.max <- getMaxWeight (object.list, attribute = c ( "idents" , "count" )) # par(mfrow = c(1,2)) 表示将绘图区域分为1行2列，用来并排显示两个图，xpd=TRUE允许文本标签在绘图区域外显示 par ( mfrow = c ( 1 , 2 ), xpd= TRUE ) # 遍历object.list中的每一个CellChat对象细胞通讯网络的计数矩阵，其中每个元素代表细胞之间的交互数量 # 将权重缩放到一个标准范围，以便使不同的网络图具有一致的显示比例 # 使用之前计算的最大权重来设定边的最大权重，以确保所有的网络图在权重上保持一致性 # 设置边的最大宽度，用来调节视觉效果，确保较强的互动路径更为突出 # 设置每个子图的标题 # 线条越粗则互作对数越多 for (i in 1 : length (object.list)) { netVisual_circle (object.list[[i]] @ net $ count, weight.scale = T, label.edge= F, edge.weight.max = weight.max[ 2 ], edge.width.max = 12 , title.name = paste0 ( "Number of interactions - " , names (object.list)[i])) } 同理，用 heatmap 也可以进行展示 # 以heatmap的形式展示第二个组别中相较于第一个组别细胞通讯发生的变化，红色为上调蓝色为下调 # 左图count 反映的是细胞间通讯的频率或次数 # 右图weight反映的是细胞间通讯的互作强度 gg1 <- netVisual_heatmap (cellchat) gg2 <- netVisual_heatmap (cellchat, measure = "weight" ) gg1 + gg2 展示特定通路展示圈图、热图和弦图 # 指定要分析的细胞通讯通路是 "CXCL"，一种常见的趋化因子家族，通常涉及免疫系统中的细胞迁移 # 圈图 pathways.show <- c ( "CXCL" ) # 这行代码获取指定通路 ("CXCL") 在object.list中所有CellChat对象的最大权重 # slot.name = "netP" 表示从每个CellChat对象中获取名为netP的slot，这个slot包含了细胞间通讯的网络信息 weight.max <- getMaxWeight (object.list, slot.name = c ( "netP" ), attribute = pathways.show) par ( mfrow = c ( 1 , 2 ), xpd= TRUE ) # 将最大权重设置为先前计算的最大值，以确保不同数据集之间的权重尺度一致 # 设置边的最大宽度为10，影响图中边的粗细，帮助突出重要的通讯路径 # signaling.name设置图的标题,显示的是当前通路和数据集的名称 for (i in 1 : length (object.list)) { netVisual_aggregate (object.list[[i]], signaling = pathways.show, layout = "circle" , edge.weight.max = weight.max[ 1 ], edge.width.max = 10 , signaling.name = paste (pathways.show, names (object.list)[i])) } # 热图 pathways.show <- c ( "CXCL" ) par ( mfrow = c ( 1 , 2 ), xpd= TRUE ) ht <- list # 热图展示 "CXCL" 通路在不同细胞类型之间的通讯强度 # 设置热图颜色方案为红色渐变，用于显示通讯强度的高低 # 为每个热图设置标题，显示通路名称和数据集名称 for (i in 1 : length (object.list)) { ht[[i]] <- netVisual_heatmap (object.list[[i]], signaling = pathways.show, color.heatmap = "Reds" , title.name = paste (pathways.show, "signaling " , names (object.list)[i])) } # 保存在ht列表中的热图绘制在同一面板上，ht_gap = unit(0.5, "cm") 设置热图之间的间距为 0.5 厘米 ComplexHeatmap :: draw (ht[[ 1 ]] + ht[[ 2 ]], ht_gap = unit ( 0.5 , "cm" )) # 弦图 # 更粗的箭头 表示该细胞对之间的通讯更强 pathways.show <- c ( "CXCL" ) par ( mfrow = c ( 1 , 2 ), xpd= TRUE ) for (i in 1 : length (object.list)) { netVisual_aggregate (object.list[[i]], signaling = pathways.show, layout = "chord" , signaling.name = paste (pathways.show, names (object.list)[i])) } 基因表达情况 # 设定图形中组别展示顺序，NL在前，LS在后 cellchat @ meta $ datasets = factor (cellchat @ meta $ datasets, levels = c ( "NL" , "LS" )) # 指定了要可视化的信号通路为CXCL # 划分不同的数据集进行对比展示 # 在可视化中使用 ggplot2 的默认配色方案 plotGeneExpression (cellchat, signaling = "CXCL" , split.by = "datasets" , colors.ggplot = T) 研究细胞类型之间的通讯差异 # 检查这两个数据集的细胞类型（idents）是否一致 levels (object.list[[ 1 ]] @ idents) # [1] "APOE+ FIB" "FBN1+ FIB" "COL11A1+ FIB" "Inflam. FIB" "cDC1" # [6] "cDC2" "LC" "Inflam. DC" "TC" "Inflam. TC" # [11] "CD40LG+ TC" "NKT" levels (object.list[[ 2 ]] @ idents) # [1] "APOE+ FIB" "FBN1+ FIB" "COL11A1+ FIB" "Inflam. FIB" "cDC1" # [6] "cDC2" "LC" "Inflam. DC" "TC" "Inflam. TC" # [11] "CD40LG+ TC" "NKT" # 定义了一个包含细胞类型信息的向量。

这里假设有三种细胞类型（"FIB"、"DC"、"TC"），并且每种类型有4种细胞 group.cellType <- c ( rep ( "FIB" , 4 ), rep ( "DC" , 4 ), rep ( "TC" , 4 )) # 将 group.cellType 转换为因子，并指定因子的顺序为"FIB" -> "DC" -> "TC" group.cellType <- factor (group.cellType, levels = c ( "FIB" , "DC" , "TC" )) # 将group.cellType作为新的细胞类型信息合并到每个CellChat对象中。

这会更新object.list中的每个 CellChat 对象，确保每个对象的细胞类型标签被正确地合并 object.list <- lapply (object.list, function (x) { mergeInteractions (x, group.cellType)}) # 将object.list中的所有CellChat对象合并成一个单一的CellChat对象。

这一步是将不同数据集（例如，不同实验条件下的数据）进行合并，便于进行跨数据集的比较 cellchat <- mergeCellChat (object.list, add.names = names (object.list)) # 直接展示相减的结果 # 展示基于合并后的细胞通讯计数（count.merged）的差异 par ( mfrow = c ( 1 , 2 ), xpd= TRUE ) netVisual_diffInteraction (cellchat, weight.scale = T, measure = "count.merged" , label.edge = T) # 展示基于合并后的细胞通讯强度（weight.merged）的差异 netVisual_diffInteraction (cellchat, weight.scale = T, measure = "weight.merged" , label.edge = T) # 两组分别展示 # getMaxWeight 函数用于获取object.list中指定槽位和属性的最大权重，这里分别考虑了细胞身份、相互作用计数和合并后的相互作用计数 weight.max <- getMaxWeight (object.list, slot.name = c ( "idents" , "net" , "net" ), attribute = c ( "idents" , "count" , "count.merged" )) par ( mfrow = c ( 1 , 2 ), xpd= TRUE ) # 使用for循环遍历 object.list中的每个对象，netVisual_circle函数以圆形布局可视化细胞间相互作用网络 # 使用合并后的相互作用计数count.merged作为边的权重 # weight.scale=T表示对权重进行缩放 # label.edge=T表示显示边的标签 # edge.weight.max 设置边的最大权重为之前获取的最大权重 # edge.width.max设置边的最大宽度为12 # title.name 为每个图添加标题，显示相互作用数量和对象名称 for (i in 1 : length (object.list)) { netVisual_circle (object.list[[i]] @ net $ count.merged, weight.scale = T, label.edge= T, edge.weight.max = weight.max[ 3 ], edge.width.max = 12 , title.name = paste0 ( "Number of interactions - " , names (object.list)[i])) } 展示两组间各类细胞 incoming 与 outcoming 通讯的强度 # x@net$count为CellChat对象中存储细胞间通讯计数的矩阵。

矩阵中的每个元素表示一个细胞对之间的通讯计数。# rowSums(x@net$count)：计算每行（细胞源）上的总通讯数量，即每个细胞源的所有目标细胞的通讯总数。# colSums(x@net$count)：计算每列（细胞目标）上的总通讯数量，即每个目标细胞的所有源细胞的通讯总数。

# diag(x@net$count)取矩阵的对角线元素，这些通常表示细胞本身与自己的通讯，需要将这些对角线元素减去，因为它们不代表细胞间的通讯 # num.link：这会输出一个包含每个CellChat对象中每个细胞类型的通讯总数 num.link <- sapply (object.list, function (x) { rowSums (x @ net $ count) + colSums (x @ net $ count) - diag (x @ net $ count)}) num.link # weight.MinMax 存储了所有对象中链接数量的最小值和最大值，用于控制不同数据集散点图中散点的大小 weight.MinMax <- c ( min (num.link), max (num.link)) gg <- list # 生成每个CellChat对象不同细胞类型通讯总数的散点图 for (i in 1 : length (object.list)) { gg[[i]] <- netAnalysis_signalingRole_scatter (object.list[[i]], title = names (object.list)[i], weight.MinMax = weight.MinMax) } p1 <- patchwork :: wrap_plots ( plots = gg) png ( "/data/07_scRNA-seq_cell_communication/img/cellchat.png" , width = 720 , height = 360 ) print (p1) dev.off 展示两组间的差异 # 对 CellChat 对象中的两种不同细胞类型（"Inflam. DC" 和 "cDC1"）进行信号传导变化的散点图分析，同时排除 "MIF" 信号通路的影响，最后将两个散点图组合在一起，以便直观地比较这两种细胞类型在信号传导方面的差异 gg1 <- netAnalysis_signalingChanges_scatter (cellchat, idents.use = "Inflam. DC" , signaling.exclude = "MIF" ) gg2 <- netAnalysis_signalingChanges_scatter (cellchat, idents.use = "cDC1" , signaling.exclude = c ( "MIF" )) patchwork :: wrap_plots ( plots = list (gg1,gg2)) 继续探索outgoing与incoming的组间模式差异 suppressMessages ( library (ComplexHeatmap)) i = 1 # outgoing 传出信号，作为发送者角色即配体 # 使用union函数将这两个CellChat对象的信号通路合并，得到这两个数据集中所有信号通路的并集 pathway.union <- union (object.list[[i]] @ netP $ pathways, object.list[[i + 1 ]] @ netP $ pathways) # 生成信号通路角色的热图，将展示在每个数据集中，细胞在信号通路中的“outgoing”角色（即发送信号的角色） ht1 = netAnalysis_signalingRole_heatmap (object.list[[i]], pattern = "outgoing" , signaling = pathway.union, title = names (object.list)[i], width = 5 , height = 6 ) ht2 = netAnalysis_signalingRole_heatmap (object.list[[i + 1 ]], pattern = "outgoing" , signaling = pathway.union, title = names (object.list)[i + 1 ], width = 5 , height = 6 ) draw (ht1 + ht2, ht_gap = unit ( 0.5 , "cm" )) # incoming 传入信号，作为接收者即受体 ht1 = netAnalysis_signalingRole_heatmap (object.list[[i]], pattern = "incoming" , signaling = pathway.union, title = names (object.list)[i], width = 5 , height = 6 , color.heatmap = "GnBu" ) ht2 = netAnalysis_signalingRole_heatmap (object.list[[i + 1 ]], pattern = "incoming" , signaling = pathway.union, title = names (object.list)[i + 1 ], width = 5 , height = 6 , color.heatmap = "GnBu" ) draw (ht1 + ht2, ht_gap = unit ( 0.5 , "cm" )) # overall # 这次展示的是 overall 角色，即展示细胞在信号通路中的所有角色（既包括发送信号的角色，又包括接收信号的角色） ht1 = netAnalysis_signalingRole_heatmap (object.list[[i]], pattern = "all" , signaling = pathway.union, title = names (object.list)[i], width = 5 , height = 6 , color.heatmap = "OrRd" ) ht2 = netAnalysis_signalingRole_heatmap (object.list[[i + 1 ]], pattern = "all" , signaling = pathway.union, title = names (object.list)[i + 1 ], width = 5 , height = 6 , color.heatmap = "OrRd" ) draw (ht1 + ht2, ht_gap = unit ( 0.5 , "cm" )) 从配受体的角度画一画 气泡图 # sources.use指定了源细胞类型的编号，4表示选择编号为4的细胞类型作为信号源。

在气泡图中，源细胞是负责发送信号的细胞类型 # targets.use指定了目标细胞类型的编号，c(5:11) 表示选择编号为5到11的细胞类型作为信号的接收者 # comparison指定需要比较的两个数据集或实验条件。c(1,2)表示将会比较数据集或实验组1和实验组2中细胞通讯的差异 # angle.x 控制 X 轴标签的旋转角度。

在这里，设置为45度，目的是使X轴标签的文字更容易阅读，尤其是当标签较长时 netVisual_bubble (cellchat, sources.use = 4 , targets.use = c ( 5 : 11 ), comparison = c ( 1 , 2 ), angle.x = 45 ) # 我们可以识别一个数据集中与另一个数据集相比上调 (增加) 和下调 (减少) 的信号配体 - 受体对。

这可以通过在 netVisual_bubble函数中指定max.dataset和min.dataset来实现。信号增加意味着这些信号在一个数据集中比另一个数据集具有更高的通信概率 (强度)。

gg1 <- netVisual_bubble (cellchat, sources.use = 4 , targets.use = c ( 5 : 11 ), comparison = c ( 1 , 2 ), max.dataset = 2 , title.name = "Increased signaling in LS" , angle.x = 45 , remove.isolate = T) gg2 <- netVisual_bubble (cellchat, sources.use = 4 , targets.use = c ( 5 : 11 ), comparison = c ( 1 , 2 ), max.dataset = 1 , title.name = "Decreased signaling in LS" , angle.x = 45 , remove.isolate = T) gg1 + gg2 以上的计算依赖的都是细胞通讯的可能性(强度)，接下来我们将通过配受体对基因表达的角度进一步研究 pos.dataset = "LS" features.name = pos.dataset # 差异计算 # identifyOverExpressedGenes函数用于识别在LS组中过表达的基因 # group.dataset = "datasets"：指定用于比较的分组（如实验组、对照组等）。

这个值指的是 cellchat@meta$datasets 中的分组信息 # only.pos = FALSE：表示不仅仅考虑上调基因，还包括下调基因。默认是 FALSE，表示同时考虑上调和下调基因 # thresh.pc = 0.1：指定基因在该组中需要过表达的最小百分比。

如果某个基因在指定的实验组中过表达的细胞占比超过10%，则认为该基因被过表达 # thresh.fc = 0.1：指定基因在该组中的最小对数变化倍数（log fold change）。

如果基因的对数倍数变化大于 0.1，则认为该基因过表达 # thresh.p = 1：p 值的阈值，设置为 1 意味着不对 p 值进行过滤 cellchat <- identifyOverExpressedGenes (cellchat, group.dataset = "datasets" , pos.dataset = pos.dataset, features.name = features.name, only.pos = FALSE , thresh.pc = 0.1 , thresh.fc = 0.1 , thresh.p = 1 ) # 提取LS中上调的配受体对 # netMappingDEG 函数用于将差异表达基因（Differentially Expressed Genes, DEG）映射到细胞通讯网络上。

它会分析哪些差异表达基因参与了细胞间的信号通讯 net <- netMappingDEG (cellchat, features.name = features.name) # 提取LS组上调的配体 # ligand.logFC表示筛选出配体基因的对数倍数变化（logFC）大于等于 0.2 的信号通讯对，即筛选出上调的配体基因参与的信号通讯 net.up <- subsetCommunication (cellchat, net = net, datasets = "LS" , ligand.logFC = 0.2 , receptor.logFC = NULL ) # 提取LS组下调的配体基因和下调受体基因参与的信号通讯 net.down <- subsetCommunication (cellchat, net = net, datasets = "LS" , ligand.logFC = - 0.1 , receptor.logFC = - 0.1 ) # extractGeneSubsetFromPair 函数用于从细胞通讯对中提取相关的基因子集 gene.up <- extractGeneSubsetFromPair (net.up, cellchat) gene.down <- extractGeneSubsetFromPair (net.down, cellchat) # 可视化：# 气泡图 # 从 net.up 中提取与上调信号相关的配体-受体（ligand-receptor）对的名称（interaction_name） pairLR.use.up = net.up[, "interaction_name" , drop = F] gg1 <- netVisual_bubble (cellchat, pairLR.use = pairLR.use.up, sources.use = 4 , targets.use = c ( 5 : 11 ), comparison = c ( 1 , 2 ), angle.x = 45 , remove.isolate = T, title.name = paste0 ( "Up-regulated signaling in " , names (object.list)[ 2 ])) # 从 net.down 中提取下调的配体-受体对的名称（interaction_name） pairLR.use.down = net.down[, "interaction_name" , drop = F] gg2 <- netVisual_bubble (cellchat, pairLR.use = pairLR.use.down, sources.use = 4 , targets.use = c ( 5 : 11 ), comparison = c ( 1 , 2 ), angle.x = 45 , remove.isolate = T, title.name = paste0 ( "Down-regulated signaling in " , names (object.list)[ 2 ])) gg1 + gg2 # 弦图 # lab.cex = 0.8：标签字体大小 # small.gap = 3.5：每两个细胞群体之间的间隙大小，数值越大，间隙越明显 par ( mfrow = c ( 1 , 2 ), xpd= TRUE ) netVisual_chord_gene (object.list[[ 1 ]], sources.use = 4 , targets.use = c ( 5 : 11 ), slot.name = "net" , net = net.down, lab.cex = 0.8 , small.gap = 3.5 , title.name = paste0 ( "Down-regulated signaling in " , names (object.list)[ 1 ])) netVisual_chord_gene (object.list[[ 2 ]], sources.use = 4 , targets.use = c ( 5 : 11 ), slot.name = "net" , net = net.up, lab.cex = 0.8 , small.gap = 3.5 , title.name = paste0 ( "Up-regulated signaling in " , names (object.list)[ 2 ])) Sys.time # [1] "2025-05-15 11:15:10 CST" save.image ( "result_cellchat/cellchat.rdata" )

五、nichenetr 测试数据来源于: Medaglia et al. 等人提供的NICHE-seq数据，用于研究淋巴细胞性脑膜脊髓炎病毒( lymphocytic choriomeningitis virus, LCVM)注射后72小时在腹股沟淋巴结中是否会对T Cell的细胞通讯产生影响。NicheNet 会选择一定的基因集用于分析，因此在做 NicheNet 分析前需要确定一些受体细胞的特定基因集(通常是一些差异基因)。

1、准备工作 1.1 安装并加载R包 # 安装R包 # install.packages("mlrMBO") # if(!require(nichenetr))devtools::install_github("saeyslab/nichenetr") # 加载R包 library (nichenetr) library (Seurat) library (tidyverse) # 查看R包版本 packageVersion ( "Seurat" ) # [1] '5.0.3' 1.2 表达矩阵获取 # 读入数据 # 这个网址可能打不开 # seuratObj = readRDS(url("https://zenodo.org/record/3531889/files/seuratObj.rds")) # 我把他下载到了本地：Sys.time # [1] "2025-05-15 11:16:14 CST" rm ( list = ls ) gc # used (Mb) gc trigger (Mb) max used (Mb) # Ncells 6995062 373.6 13029127 695.9 13029127 695.9 # Vcells 13767826 105.1 84017254 641.1 205115830 1565.0 setwd ( "/data/07_scRNA-seq_cell_communication/" ) seuratObj = readRDS ( "data/seuratObj.rds" ) # 如遇报错"no slot of name "images" for this object of class "Seurat""，可以使用UpdateSeuratObject更细seurat对象 seuratObj <- UpdateSeuratObject (seuratObj) seuratObj # An object of class Seurat # 13541 features across 5027 samples within 1 assay # Active assay: RNA (13541 features, 1575 variable features) # 3 layers present: counts, data, scale.data # 4 dimensional reductions calculated: cca, cca.aligned, tsne, pca # 查看数据基本信息 # 注释信息 head (seuratObj @ meta.data) # nGene nUMI orig.ident aggregate res.0.6 celltype nCount_RNA # W380370 880 1611 LN_SS 1 CD8 T 1607 # W380372 541 891 LN_SS 0 CD4 T 885 # W3803742 1229 LN_SS 0 CD4 T 1223 # W380378 847 1546 LN_SS 1 CD8 T 1537 # W380379 839 1606 LN_SS 0 CD4 T 1603 # W380381 517 844 LN_SS 0 CD4 T 840 # nFeature_RNA # W380370 876 # W380372 536 # W380374 737 # W380378 838 # W380379 836 # W380381 513 # celltype数量 table (seuratObj @ meta.data $ celltype) # # B CD4 T CD8 T DC Mono NK Treg # 382 2562 1645 18 90 131 199 # 降维图 DimPlot (seuratObj, reduction = "tsne" ) # 查看分组信息，这里分组变量为"aggregate" seuratObj @ meta.data $ aggregate %>% table # . # LCMV SS # 3886 1141 table (seuratObj @ meta.data $ celltype, seuratObj @ meta.data $ aggregate) # # LCMV SS # B 344 38 # CD4 T 1961 601 # CD8 T 1252 393 # DC 14 4 # Mono 75 15 # NK 94 37 # Treg 146 53 DimPlot (seuratObj, reduction = "tsne" , group.by = "aggregate" ) 1.3 下载先验模型 目前支持的物种为小鼠 mouse 与人类 human # 定义分析物种 organism = "mouse" # 在线读取库文件，可能连接比较困难 if (F){ if (organism == "human" ){ lr_network = readRDS ( url ( "https://zenodo.org/record/7074291/files/lr_network_human_21122021.rds" )) ligand_target_matrix = readRDS ( url ( "https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final.rds" )) weighted_networks = readRDS ( url ( "https://zenodo.org/record/7074291/files/weighted_networks_nsga2r_final.rds" )) } else if (organism == "mouse" ){ 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)、配体-靶基因关联矩阵 (ligand_target_matrix) 和加权网络 (weighted_networks)，这些是 NicheNet 核心分析所需的数据库 if (organism == "human" ){ lr_network = readRDS ( "data/lr_network_human_21122021.rds" ) ligand_target_matrix = readRDS ( "data/ligand_target_matrix_nsga2r_final.rds" ) weighted_networks = readRDS ( "data/weighted_networks_nsga2r_final.rds" ) } else if (organism == "mouse" ){ lr_network = readRDS ( "data/lr_network_mouse_21122021.rds" ) ligand_target_matrix = readRDS ( "data/ligand_target_matrix_nsga2r_final_mouse.rds" ) weighted_networks = readRDS ( "data/weighted_networks_nsga2r_final_mouse.rds" ) } # 查看三个文件前几行 head (lr_network) # # A tibble: 6 × 4 # from to database source # <chr> <chr> <chr> <chr> # 1 2300002M23Rik Ddr1 omnipath # 2 2610528A11Rik Gpr15 omnipath # 3 9530003J23Rik Itgal omnipath # 4 a Atrn omnipath # 5 a F11r omnipath # 6 a Mc1r omnipath head (ligand_target_matrix[, 1 : 6 ]) # 2300002M23Rik 2610528A11Rik 9530003J23Rik a # 0610005C13Rik 0.000000e+00 0.000000e+00 1.311297e-05 0.000000e+00 # 0610009B22Rik 0.000000e+00 0.000000e+00 1.269301e-05 0.000000e+00 # 0610009L18Rik 8.872902e-05 4.977197e-05 2.581909e-04 7.570125e-05 # 0610010F05Rik 2.194046e-03 1.111556e-03 3.142374e-03 1.631658e-03 # 0610010K14Rik 2.271606e-03 9.360769e-04 3.546140e-03 1.697713e-03 # 0610012G03Rik 6.619036e-04 2.958259e-04 9.295219e-04 4.211287e-04 # A2m Aanat # 0610005C13Rik 1.390053e-05 0.000000e+00 # 0610009B22Rik 1.345536e-05 0.000000e+00 # 0610009L18Rik 9.802264e-05 6.938897e-05 # 0610010F05Rik 2.585820e-03 1.829664e-03 # 0610010K14Rik 2.632082e-03 1.658516e-03 # 0610012G03Rik 6.271871e-04 4.077356e-04 head (weighted_networks) # $lr_sig # # A tibble: 3,865,137 × 3 # from to weight # <chr> <chr> <dbl> # 1 0610010F05Rik App 0.110 # 2 0610010F05Rik Cat 0.0673 # 3 0610010F05Rik H1f2 0.0660 # 4 0610010F05Rik Lrrc49 0.0829 # 5 0610010F05Rik Nicn1 0.0864 # 6 0610010F05Rik Srpk1 0.123 # 7 0610010F05Rik Ttll1 0.0855 # 8 0610010K14Rik Akap8 0.131 # 9 0610010K14Rik Ash2l 0.312 # 10 0610010K14Rik Bptf 0.146 # # ℹ 3,865,127 more rows # # $gr # # A tibble: 4,364,411 × 3 # from to weight # <chr> <chr> <dbl> # 1 0610010K14Rik 0610010K14Rik 0.121 # 2 0610010K14Rik 2510039O18Rik 0.121 # 3 0610010K14Rik 2610021A01Rik 0.0256 # 4 0610010K14Rik 9130401M01Rik 0.0263 # 5 0610010K14Rik Alg1 0.127 # 6 0610010K14Rik Alox12 0.128 # 7 0610010K14Rik Ambra1 0.122 # 8 0610010K14Rik Ankib1 0.127 # 9 0610010K14Rik Anxa1 0.125 # 10 0610010K14Rik Arhgap1 0.125 # # ℹ 4,364,401 more rows # 预处理 # 去除 lr_network 数据框中 from 和 to 列组合重复的行，保留唯一的 from 和 to 组合 lr_network = lr_network %>% distinct (from, to) head (lr_network) # # A tibble: 6 × 2 # from to # <chr> <chr> # 1 2300002M23Rik Ddr1 # 2 2610528A11Rik Gpr15 # 3 9530003J23Rik Itgal # 4 a Atrn # 5 a F11r # 6 a Mc1r weighted_networks_lr = weighted_networks $ lr_sig %>% inner_join (lr_network, by = c ( "from" , "to" )) # 1. lr_sig（配体-受体信号网络）：包含了每个配体-受体对的权重，这些权重通常反映了配体通过其受体对靶基因表达的影响程度。

# from: 配体（ligand） # to: 受体（receptor） # weight: 配体-受体相互作用的权重，通常是信号强度或潜在调控能力 head (weighted_networks $ lr_sig) # # A tibble: 6 × 3 # from to weight # <chr> <chr> <dbl> # 1 0610010F05Rik App 0.110 # 2 0610010F05Rik Cat 0.0673 # 3 0610010F05Rik H1f2 0.0660 # 4 0610010F05Rik Lrrc49 0.0829 # 5 0610010F05Rik Nicn1 0.0864 # 6 0610010F05Rik Srpk1 0.123 # 2. gr（受体-基因调控网络）：描述的是受体与靶基因之间的调控关系，即受体激活后，如何通过转录因子或其他调控机制影响靶基因的表达。

它揭示了受体如何调控基因表达，是配体通过其受体影响靶基因的下游环节。通过 gr 网络，NicheNet 能够进一步推测从配体到基因的调控路径。

# from: 受体（receptor） # to: 靶基因（target gene） # weight: 受体对靶基因表达的影响程度或调控强度 head (weighted_networks $ gr) # # A tibble: 6 × 3 # from to weight # <chr> <chr> <dbl> # 1 0610010K14Rik 0610010K14Rik 0.121 # 2 0610010K14Rik 2510039O18Rik 0.121 # 3 0610010K14Rik 2610021A01Rik 0.0256 # 4 0610010K14Rik 9130401M01Rik 0.0263 # 5 0610010K14Rik Alg1 0.127 # 6 0610010K14Rik Alox12 0.128 2、细胞通讯分析 2.1 定义sender与receiver 这里作者的研究目的中，CD8 T 作为 receiver，CD4 T , Treg , Mono , NK , B 与 DC 作为 sender。

并且在表达矩阵中，只有至少在一种细胞内表达比例至少达到10%的基因才会参与下游分析。

# receiver = "CD8 T" # 筛选在 seuratObj 这个 Seurat对象中 CD8 T细胞 群中，至少10%的细胞表达的基因 expressed_genes_receiver = get_expressed_genes (receiver, seuratObj, pct = 0.10 ) # 从 expressed_genes_receiver 中，进一步筛选出在 ligand_target_matrix 中也存在的基因 background_expressed_genes = expressed_genes_receiver %>% .[. %in% rownames (ligand_target_matrix)] # sender_celltypes = c ( "CD4 T" , "Treg" , "Mono" , "NK" , "B" , "DC" ) # 由于这里的sender细胞种类大于1，因此通过lapply获取每种细胞类型中表达比例大于10%的基因 list_expressed_genes_sender = sender_celltypes %>% unique %>% lapply (get_expressed_genes, seuratObj, 0.10 ) expressed_genes_sender = list_expressed_genes_sender %>% unlist %>% unique 2.2 定义一组感兴趣的基因集 receiver 细胞类型中可能存在经细胞互作影响表达的基因，这些基因可帮助后续的数据分析更加的有目的性。

由于这里作者需要对比LCVM与steady-state(SS)间 CD8 T 的差别，因此这里的基因集便为T Cell在组间的差异基因。

# 取出受体细胞的Seurat对象 seurat_obj_receiver = subset (seuratObj, idents = receiver) # 将分组信息aggregate作为Idents (seurat_obj_receiver) <- "aggregate" # 实验组：condition_oi <- "LCMV" # 对照组 condition_reference <- "SS" # 计算组间差异基因 # rownames_to_column("gene") 将基因名从行名转换为新的一列 gene DE_table_receiver <- FindMarkers ( object = seurat_obj_receiver, ident.1 = condition_oi, ident.2 = condition_reference, min.pct = 0.10 ) %>% rownames_to_column ( "gene" ) head (DE_table_receiver) # gene p_val avg_log2FC pct.1 pct.2 p_val_adj # 1 Ifi27l2b 2.164082e-150 2.922767 0.973 0.387 2.930383e-146 # 2 Irf7 5.425060e-129 3.707013 0.870 0.165 7.346074e-125 # 3 Ly6a 8.023340e-124 3.358513 0.879 0.239 1.086440e-119 # 4 Ifi27l2a 2.882324e-115 3.116906 0.861 0.191 3.902955e-111 # 5 Stat1 2.742811e-92 2.210800 0.888 0.389 3.714040e-88 # 6 Ly6c2 3.421692e-89 1.907967 0.940 0.534 4.633313e-85 # 过滤并获得基因名称 # 过滤标准：多重校正后的 p 值 ≤ 0.05（显著性）和 log2FoldChange 的绝对值 ≥ 0.25 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)] # 最后查看一下前五个基因 geneset_oi[ 1 : 5 ] # [1] "Ifi27l2b" "Irf7" "Ly6a" "Stat1" "Ly6c2" 2.3 获取配受体基因表达矩阵 结合前面获得的 sender 与 receiver 的表达矩阵与库文件，分别取出受体基因在 receiver 的表达矩阵与配体基因在 sender 中对应的表达矩阵：# 从库文件中取出配体基因 ligands = lr_network %>% pull (from) %>% unique ligands[ 1 : 5 ] # [1] "2300002M23Rik" "2610528A11Rik" "9530003J23Rik" "a" # [5] "A2m" # 从库文件中取出受体基因 receptors = lr_network %>% pull (to) %>% unique receptors[ 1 : 5 ] # [1] "Ddr1" "Gpr15" "Itgal" "Atrn" "F11r" # 在`sender`发送细胞中表达的配体基因 expressed_ligands = intersect (ligands,expressed_genes_sender) expressed_ligands[ 1 : 5 ] # [1] "Ace" "Adam10" "Adam17" "Adam23" "Aimp1" # 在`receiver`接收细胞表达的受体基因 expressed_receptors = intersect (receptors,expressed_genes_receiver) expressed_receptors[ 1 : 5 ] # [1] "Itgal" "Notch1" "Itga4" "Itgb7" "Treml2" # 并且这个配体-受体对本身在 lr_network 网络中是存在的 # 如果存在则取出真正的配体基因作为潜在配体 potential_ligands = lr_network %>% filter (from %in% expressed_ligands & to %in% expressed_receptors) %>% pull (from) %>% unique potential_ligands[ 1 : 5 ] # [1] "Adam10" "Adam17" "App" "B2m" "Btla" 2.4 配体活性分析 结合配体的靶向基因在作者感兴趣的基因集和背景基因的表达情况，对潜在配体进行评级。

# 核心逻辑是：对每一个 potential_ligand，看它根据 ligand_target_matrix 最可能调控 geneset_oi 中的基因。

计算每个配体的一个打分指标，衡量它能否解释目标基因集的表达变化 # geneset：是你在前面筛选出来的目标基因集合（"CD8 T"组件差异表达且在ligand_target_matrix中存在的基因），即配体可能调控的基因 # background_expressed_genes：是在接收细胞中表达的背景基因，用于设定合理的比较基线 # ligand_target_matrix：是预先训练好的配体-靶基因打分矩阵（每个配体对每个基因的调控潜力打分） # potential_ligands：是之前筛选出的、在发送细胞表达且对应受体在接收细胞表达的潜在配体 ligand_activities = predict_ligand_activities ( geneset = geneset_oi, background_expressed_genes = background_expressed_genes, ligand_target_matrix = ligand_target_matrix, potential_ligands = potential_ligands ) # 按照 aupr_corrected（打分指标,adjusted Area Under the Precision-Recall curve，一种综合考虑准确率和召回率的打分）从高到低排序 # mutate为每个配体生成一个排名,并存储到ligand_activities的rank列 ligand_activities = ligand_activities %>% arrange ( - aupr_corrected) %>% mutate ( rank = rank ( desc (aupr_corrected))) head (ligand_activities) # # A tibble: 6 × 6 # test_ligand auroc aupr_corrected pearson rank # <chr> <dbl> <dbl> <dbl> <dbl> <dbl> # 1 Ebi3 0.655 0.380 0.234 0.291 1 # 2 Ptprc 0.640 0.304 0.158 0.160 2 # 3 H2-M3 0.609 0.287 0.141 0.181 3 # 4 H2-M2 0.613 0.272 0.125 0.146 5 # 5 H2-T10 0.613 0.272 0.125 0.146 5 # 6 H2-T22 0.613 0.272 0.125 0.146 5 基于area under the precision-recall curve(AUPR)选择Top30的ligand用户下游分析 # 把筛选出来的30个配体再按 aupr_corrected 分数从高到低排一下顺序 # 从 ligand_activities 中提取出列名叫 test_ligand 的这一列（即配体的名字），不带其他列 best_upstream_ligands = ligand_activities %>% top_n ( 30 , aupr_corrected) %>% arrange ( - aupr_corrected) %>% pull (test_ligand) %>% unique best_upstream_ligands # [1] "Ebi3" "Ptprc" "H2-M3" "H2-M2" "H2-T10" "H2-T22" "H2-T23" "H2-K1" # [9] "H2-Q4" "H2-Q6" "H2-Q7" "H2-D1" "Sirpa" "Cd48" "App" "Tgfb1" # [17] "Ccl22" "Selplg" "Btla" "Cxcl10" "Adam17" "Cxcl11" "Icam1" "Tgm2" # [25] "Cxcl9" "B2m" "Cd72" "C3" "Hp" "Itgb2" # 查看配体在各细胞中的表达情况 DotPlot (seuratObj, features = best_upstream_ligands %>% rev , cols = "RdYlBu" ) + RotatedAxis 2.5 推断受体以及配体潜在调控靶基因与配体的调控活性 2.5.1 配体与Target Gene的调控关系 # get_weighted_ligand_target_links获取配体对应的靶基因打分矩阵 # 对每个配体，从 ligand_target_matrix 中选出打分最高的 200 个目标基因（如果这些目标基因在你的 geneset_oi 中） active_ligand_target_links_df = best_upstream_ligands %>% lapply (get_weighted_ligand_target_links, geneset = geneset_oi, ligand_target_matrix = ligand_target_matrix, n = 200 ) %>% bind_rows %>% drop_na # 得到的 active_ligand_target_links_df 是一个数据框，其中包含每个配体对多个靶基因的“调控潜力得分” head (active_ligand_target_links_df) # # A tibble: 6 × 3 # ligand target weight # <chr> <chr> <dbl> # 1 Ebi3 Bst2 0.0500 # 2 Ebi3 Cd274 0.0504 # 3 Ebi3 Cxcl10 0.0570 # 4 Ebi3 Cxcr4 0.0430 # 5 Ebi3 Ddit4 0.0485 # 6 Ebi3 Ddx58 0.0402 # 将上一步的 dataframe 转换成适合绘制热图的矩阵形式。

# cutoff = 0.33 表示过滤掉得分低于这个值的配体-靶基因对，降低噪音，突出高潜力的调控关系 active_ligand_target_links = prepare_ligand_target_visualization ( ligand_target_df = active_ligand_target_links_df, ligand_target_matrix = ligand_target_matrix, cutoff = 0.33 ) # 保留best_upstream_ligands和colnames(active_ligand_target_links)都存在的配体 # rev翻转配体基因排列顺序 # make.names 把名字里带 -、.、空格等特殊字符的基因名转换成 R 合法变量名（避免热图绘图时报错） order_ligands = intersect (best_upstream_ligands, colnames (active_ligand_target_links)) %>% rev %>% make.names # 保留active_ligand_target_links_df$target和rownames(active_ligand_target_links)都存在的受体 # 同样处理基因可能存在的特殊字符 order_targets = active_ligand_target_links_df $ target %>% unique %>% intersect ( rownames (active_ligand_target_links)) %>% make.names # 处理active_ligand_target_links行列配体和靶基因名称可能存在的特殊字符，如 rownames (active_ligand_target_links) = rownames (active_ligand_target_links) %>% make.names colnames (active_ligand_target_links) = colnames (active_ligand_target_links) %>% make.names # 得到最终的绘图矩阵 vis_ligand_target = active_ligand_target_links[order_targets,order_ligands] %>% t # 可视化：p_ligand_target_network = vis_ligand_target %>% make_heatmap_ggplot ( "Prioritized ligands" , # y轴标题为配体 "Predicted target genes" , # x轴标题为配体的靶基因 color = "purple" , legend_position = "top" , # 图例放在图形顶部 x_axis_position = "top" , # x轴标签放在图形顶部 legend_title = "Regulatory potential" ) + # 为图例设置标题，表明图例代表的是配体对靶基因的调控潜力 theme ( axis.text.x = element_text ( face = "italic" )) + # 将x轴上的文本字体设置为斜体 scale_fill_gradient2 ( low = "whitesmoke" , high = "purple" ) p_ligand_target_network # 热图中展示的结果经过严格的过滤，不包含target gene或调控分值不足0.33者均被忽略 2.5.2 配体与受体间的活性分析 # 筛选配体在 best_upstream_ligands中（预测最活跃的前30个配体）和受体在expressed_receptors中的高可信配受体对 lr_network_top = lr_network %>% filter (from %in% best_upstream_ligands & to %in% expressed_receptors) %>% distinct (from,to) # 得到最相关的受体基因名 best_upstream_receptors = lr_network_top %>% pull (to) %>% unique # 从 weighted network 中获取配受体打分 lr_network_top_df_large = weighted_networks_lr %>% filter (from %in% best_upstream_ligands & to %in% best_upstream_receptors) # 转成矩阵用于绘图（热图前的准备） # spread函数把长数据变成宽数据，lr_network_top_df_large数据框中的 "from"列的值转换为新的列名，将 "weight" 列的值填充到相应的新列中，若某个组合不存在值，则用 0 进行填充 lr_network_top_df = lr_network_top_df_large %>% spread ( "from" , "weight" , fill = 0 ) # 选取除 "to" 列之外的所有列 # 数据框转换成矩阵 # 从原始数据框 lr_network_top_df 中提取 "to" 列的值。

这些值会被用作转换后矩阵的行名 lr_network_top_matrix = lr_network_top_df %>% select ( - to) %>% as.matrix %>% magrittr :: set_rownames (lr_network_top_df $ to) # dist函数计算 lr_network_top_matrix 中受体（矩阵的行）之间的二元距离 # hclust函数对计算得到的距离矩阵进行层次聚类，采用Ward方法（ward.D2），该方法试图最小化聚类内部的方差 dist_receptors = dist (lr_network_top_matrix, method = "binary" ) hclust_receptors = hclust (dist_receptors, method = "ward.D2" ) # 从聚类结果中提取排序后的受体名称，存储在 order_receptors 中 order_receptors = hclust_receptors $ labels[hclust_receptors $ order] # 对矩阵进行转置，使得配体变为矩阵的行，受体变为矩阵的列 # 同样计算配体之间的二元距离，进行层次聚类，并提取排序后的配体名称存储在 order_ligands 中 dist_ligands = dist (lr_network_top_matrix %>% t , method = "binary" ) hclust_ligands = hclust (dist_ligands, method = "ward.D2" ) order_ligands = hclust_ligands $ labels[hclust_ligands $ order] # 使用 intersect 函数确保 order_receptors 和 order_ligands 中的名称与原始矩阵的行名和列名一致，避免出现不匹配的情况 order_receptors = order_receptors %>% intersect ( rownames (lr_network_top_matrix)) order_ligands = order_ligands %>% intersect ( colnames (lr_network_top_matrix)) # 根据排序后的受体和配体名称对原始矩阵进行重排 # make.names 函数用于将名称转换为合法的 R 名称，避免出现非法字符 vis_ligand_receptor_network = lr_network_top_matrix[order_receptors, order_ligands] rownames (vis_ligand_receptor_network) = order_receptors %>% make.names colnames (vis_ligand_receptor_network) = order_ligands %>% make.names # 可视化 # 配体在 x 轴，受体在 y 轴 p_ligand_receptor_network = vis_ligand_receptor_network %>% t %>% make_heatmap_ggplot ( "Ligands" , # y轴标题为配体 "Receptors" , # x轴标题为受体 color = "mediumvioletred" , x_axis_position = "top" , legend_title = "Prior interaction potential" ) p_ligand_receptor_network 2.6 sender 中配体的组间差异log fold change值 # 对每个发送细胞类型调用 get_lfc_celltype，进行LCMV vs SS组间差异分析 DE_table_all = Idents (seuratObj) %>% levels %>% intersect (sender_celltypes) %>% lapply (get_lfc_celltype, seurat_obj = seuratObj, condition_colname = "aggregate" , condition_oi = condition_oi, condition_reference = condition_reference, # 注意修改自己的对照组与实验组名称 expression_pct = 0.10 , # 表达比例大于 10% 才考虑 celltype_col = NULL ) %>% reduce (full_join) # 将所有发送细胞类型的结果合并为一个大表（按基因名对齐） # 将非差异基因者数值转换为0 DE_table_all[ is.na (DE_table_all)] = 0 head (DE_table_all) # # A tibble: 6 × 7 # gene `CD4 T` Treg B NK Mono DC # <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> # 1 Ifi27l2b 2.55 1.64 1.70 3.36 2.54 6.43 # 2 Irf7 4.16 4.30 3.54 3.47 4.65 6.16 # 3 Ly6a 2.67 2.94 1.47 2.81 1.59 4.33 # 4 Ifi27l2a 2.68 2.24 1.70 3.73 2.35 5.35 # 5 Ifit3 4.69 4.47 3.41 4.55 5.81 3.74 # 6 Stat1 2.27 2.04 1.88 2.00 1.75 6.10 # 合并配体活动信息和差异表达信息 # 从 ligand_activities 中选出 test_ligand 和 pearson（预测的调控强度），重命名为 ligand # 将其和刚才的 DE_table_all 按 ligand 字段做左连接 ligand_activities_de = ligand_activities %>% select (test_ligand, pearson) %>% rename ( ligand = test_ligand) %>% left_join (DE_table_all %>% rename ( ligand = gene)) # 再次将 NA 转为 0（说明某些配体没出现在差异表达结果中） ligand_activities_de[ is.na (ligand_activities_de)] = 0 head (ligand_activities_de) # # A tibble: 6 × 8 # ligand pearson `CD4 T` Treg B NK Mono DC # <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> # 1 Ebi3 0.291 0 4.23 0.958 0.663 3.07 0 # 2 Ptprc 0.160 0.314 0.519 0.448 0.616 0.253 3.87 # 3 H2-M3 0.181 0.944 0.676 1.23 0.178 0.990 2.54 # 4 H2-M2 0.146 0 0 0 0 0 -0.705 # 5 H2-T10 0.146 1.52 1.76 1.30 0 1.60 1.62 # 6 H2-T22 0.146 1.84 2.79 2.07 1.95 3.71 5.27 # 绘制热图 # 构建用于绘图的矩阵（logFC 矩阵） # 提取差异表达数据（去掉 ligand 名和 pearson 列） # 设置行名为配体名 lfc_matrix = ligand_activities_de %>% select ( - ligand, - pearson) %>% as.matrix %>% magrittr :: set_rownames (ligand_activities_de $ ligand) # make.names 防止非法字符导致绘图报错 rownames (lfc_matrix) = rownames (lfc_matrix) %>% make.names # 确保配体顺序与之前热图一致 order_ligands = order_ligands[order_ligands %in% rownames (lfc_matrix)] vis_ligand_lfc = lfc_matrix[order_ligands,] # 修正列名即细胞名称 colnames (vis_ligand_lfc) = vis_ligand_lfc %>% colnames %>% make.names p_ligand_lfc = vis_ligand_lfc %>% make_threecolor_heatmap_ggplot ( "Prioritized ligands" , "LFC in Sender" , low_color = "midnightblue" , mid_color = "white" , mid = median (vis_ligand_lfc), high_color = "red" , legend_position = "top" , x_axis_position = "top" , legend_title = "LFC" ) + theme ( axis.text.y = element_text ( face = "italic" )) p_ligand_lfc 2.7 可视化拼盘 # 配体活性热图 # 筛选出ligand_activities的aupr_corrected，数值表示某配体调控其靶基因集合的能力（越高代表越有效） # 将ligand_activities$test_ligand设置为结果的行名 ligand_aupr_matrix = ligand_activities %>% select (aupr_corrected) %>% as.matrix %>% magrittr :: set_rownames (ligand_activities $ test_ligand) head (ligand_aupr_matrix) # aupr_corrected # Ebi3 0.2335898 # Ptprc 0.1579384 # H2-M3 0.1406565 # H2-M2 0.1254201 # H2-T10 0.1254201 # H2-T22 0.1254201 # 确保配体名合法 rownames (ligand_aupr_matrix) = rownames (ligand_aupr_matrix) %>% make.names colnames (ligand_aupr_matrix) = colnames (ligand_aupr_matrix) %>% make.names # 设置配体展示顺序和前面热图一直 # 设置列名称为AUPR vis_ligand_aupr = ligand_aupr_matrix[order_ligands, ] %>% as.matrix ( ncol = 1 ) %>% magrittr :: set_colnames ( "AUPR" ) # 绘制热图 p_ligand_aupr = vis_ligand_aupr %>% make_heatmap_ggplot ( "Prioritized ligands" , "Ligand activity" , color = "darkorange" , legend_position = "top" , x_axis_position = "top" , legend_title = "AUPR \n (target gene prediction ability)" ) + theme ( legend.text = element_text ( size = 9 )) # 配体表达量气泡图 order_ligands_adapted = order_ligands # 因为前面用 make.names 把基因名中的点号转换成了 .，这里为了让 DotPlot 正常识别表达矩阵中真实的基因名，将部分基因名改回来（如 "H2.M3" → "H2-M3" order_ligands_adapted[order_ligands_adapted == "H2.M3" ] = "H2-M3" order_ligands_adapted[order_ligands_adapted == "H2.T23" ] = "H2-T23" # 绘制气泡图 rotated_dotplot = DotPlot (seuratObj %>% subset (celltype %in% sender_celltypes), features = order_ligands_adapted, cols = "RdYlBu" ) + coord_flip + theme ( legend.text = element_text ( size = 10 ), legend.title = element_text ( size = 12 )) # 将四张关键的配体-受体分析图整合成一个完整面板 figures_without_legend = cowplot :: plot_grid ( # p_ligand_aupr：配体活性热图（AUPR 值） p_ligand_aupr + theme ( legend.position = "none" , axis.ticks = element_blank ) + theme ( axis.title.x = element_text ), # rotated_dotplot：配体表达量气泡图（在 Sender 细胞中） rotated_dotplot + theme ( legend.position = "none" , axis.ticks = element_blank , axis.title.x = element_text ( size = 12 ), axis.text.y = element_text ( face = "italic" , size = 9 ), axis.text.x = element_text ( size = 9 , angle = 90 , hjust = 0 )) + ylab ( "Expression in Sender" ) + xlab ( "" ) + scale_y_discrete ( position = "right" ), # p_ligand_lfc：Sender 细胞中配体的差异表达（LogFC） p_ligand_lfc + theme ( legend.position = "none" , axis.ticks = element_blank ) + theme ( axis.title.x = element_text ) + ylab ( "" ), # p_ligand_target_network：配体-靶基因调控网络热图 p_ligand_target_network + theme ( legend.position = "none" , axis.ticks = element_blank ) + ylab ( "" ), align = "hv" , # 水平方向和垂直方向都对齐 nrow = 1 , # 按列宽比例调节四张图的视觉宽度 rel_widths = c ( ncol (vis_ligand_aupr) + 6 , ncol (vis_ligand_lfc) + 7 , ncol (vis_ligand_lfc) + 8 , ncol (vis_ligand_target)) ) # legends：抽出各个图的图例 legends = cowplot :: plot_grid ( ggpubr :: as_ggplot (ggpubr :: get_legend (p_ligand_aupr)), ggpubr :: as_ggplot (ggpubr :: get_legend (rotated_dotplot)), ggpubr :: as_ggplot (ggpubr :: get_legend (p_ligand_lfc)), ggpubr :: as_ggplot (ggpubr :: get_legend (p_ligand_target_network)), nrow = 1 , = 1） align = "h" , # 水平对齐 rel_widths = c ( 1.5 , 1 , 1 , 1 )) # 使用 rel_widths 控制每个图例所占的空间比例 # 最终组合：主图 + 图例 combined_plot = cowplot :: plot_grid ( figures_without_legend, legends, rel_heights = c ( 10 , 5 ), nrow = 2 , align = "hv" ) # 水平方向和垂直方向都对齐 combined_plot 3、差异细胞通讯分析 本部分教程的示例数据来自于一个肝脏 scRNA-Seq 数据集( https://www.sciencedirect.com/science/article/pii/S0092867421014811)的 一部分。

NicheNet 的差异分析目标是找到配受体对均存在差异表达且激活的互作通路。

3.1 准备工作 3.1.1 加载R包 library (nichenetr) library (RColorBrewer) library (tidyverse) library (Seurat) 3.1.2 测试数据 # 读取表达矩阵，依旧是Seurat对象的rds文件 # 这个网址大家可能连接不上 if (F){ seurat_obj = readRDS ( url ( "https://zenodo.org/record/4675430/files/seurat_obj_hnscc.rds" )) DimPlot (seurat_obj, group.by = "celltype" , label = TRUE ) } # 依旧是下在本地操作比较稳：seurat_obj = readRDS ( "data/seurat_obj_hnscc.rds" ) # 查看数据中包含的细胞类型：DimPlot (seurat_obj, group.by = "celltype" , label = TRUE ) seurat_obj = SetIdent (seurat_obj, value = "celltype" ) # 查看数据中的分组信息 # p-EMT(partial epithelial-to-mesenchymal transition) 部分上皮细胞到间质的转变 DimPlot (seurat_obj, group.by = "pEMT" ) table (seurat_obj @ meta.data $ celltype, seurat_obj @ meta.data $ pEMT) # # High Low # CAF 396 104 # Endothelial 105 53 # Malignant 1093 549 # Myeloid 92 7 # myofibroblast 382 61 # T.cell 689 3 # 类似于Seurat中的差异计算，这里也需要构建一个包含celltype信息与分组信息的idents：seurat_obj @ meta.data $ celltype_aggregate <- paste (seurat_obj @ meta.data $ celltype, seurat_obj @ meta.data $ pEMT, sep = "_" ) # user adaptation required on own dataset DimPlot (seurat_obj, group.by = "celltype_aggregate" ) seurat_obj @ meta.data $ celltype_aggregate %>% table %>% sort ( decreasing = TRUE ) # . # Malignant_High T.cell_High Malignant_Low CAF_High # 1093 689 549 396 # myofibroblast_High Endothelial_High CAF_Low Myeloid_High # 382 105 104 92 # myofibroblast_Low Endothelial_Low Myeloid_Low T.cell_Low # 61 53 7 3 celltype_id <- "celltype_aggregate" Idents (seurat_obj) <- celltype_id 3.1.3 库文件 # 读取库文件 organism <- "human" # 测试数据是一个人类数据 if (T){ if (organism == "human" ){ lr_network = readRDS ( "data/lr_network_human_21122021.rds" ) ligand_target_matrix = readRDS ( "data/ligand_target_matrix_nsga2r_final.rds" ) weighted_networks = readRDS ( "data/weighted_networks_nsga2r_final.rds" ) } else if (organism == "mouse" ){ lr_network = readRDS ( "data/lr_network_mouse_21122021.rds" ) ligand_target_matrix = readRDS ( "data/ligand_target_matrix_nsga2r_final_mouse.rds" ) weighted_networks = readRDS ( "data/weighted_networks_nsga2r_final_mouse.rds" ) } } # 预处理 lr_network = lr_network %>% dplyr :: rename ( ligand = from, receptor = to) %>% distinct (ligand, receptor) head (lr_network) # # A tibble: 6 × 2 # ligand receptor # <chr> <chr> # 1 A2M MMP2 # 2 A2M MMP9 # 3 A2M LRP1 # 4 A2M KLK3 # 5 AANAT MTNR1A # 6 AANAT MTNR1B ligand_target_matrix[ 1 : 5 , 1 : 5 ] # A2M AANAT ABCA1 ACE2 # A-GAMMA3'E 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.000000000 # A1BG 0.0018503922 0.0011108718 0.0014225077 0.0028594037 0.001139013 # A1BG-AS1 0.0007400797 0.0004677614 0.0005193137 0.0007836698 0.000375007 # A1CF 0.0024799266 0.0013026348 0.0020420890 0.0047921048 0.003273375 # A2M 0.0084693452 0.0040689323 0.0064256379 0.0105191365 0.005719199 3.2 定义 niches 即 sender 与 receiver 的细胞群体集合，在R语言中的本质为 list。

需要注意的是一个 niche 中 receiver 只能有一个，而 sender 可以包含很多。

# 例如这个包含两个niche的niches，第一个为pEMT_High_niche，第二个为pEMT_Low_niches = list ( "pEMT_High_niche" = list ( "sender" = c ( "myofibroblast_High" , "Endothelial_High" , "CAF_High" , "T.cell_High" , "Myeloid_High" ), "receiver" = c ( "Malignant_High" )), "pEMT_Low_niche" = list ( "sender" = c ( "myofibroblast_Low" , "Endothelial_Low" , "CAF_Low" ), "receiver" = c ( "Malignant_Low" )) ) niches # $pEMT_High_niche # $pEMT_High_niche$sender # [1] "myofibroblast_High" "Endothelial_High" "CAF_High" # [4] "T.cell_High" "Myeloid_High" # # $pEMT_High_niche$receiver # [1] "Malignant_High" # # # $pEMT_Low_niche # $pEMT_Low_niche$sender # [1] "myofibroblast_Low" "Endothelial_Low" "CAF_Low" # # $pEMT_Low_niche$receiver # [1] "Malignant_Low" 3.3 计算 niches 的差异基因 # 目标：将配体和受体在发送者（sender）和接收者（receiver）细胞群中的差异表达信息与配对的配体-受体网络（lr_network）整合起来，计算每对LR对的特异性评分（specificity score），以找出在不同生态位（niche）中更具特异性的细胞通讯对 # 设置后续的分析操作（例如差异表达分析）基于 RNA 表达矩阵来开展 assay_oi = "RNA" # lr_network$ligand %>% unique：从 lr_network 数据结构里提取所有独特的配体名称 # subset(features = ...)：从seurat_obj中选取仅包含这些配体基因的子集。

这样做的目的是把分析范围限定在配体相关的基因上，减少不必要的计算量 # niches：一个列表，定义了不同的细胞生态位（niche），其中每个生态位包含发送者细胞群和接收者细胞群的信息 # type = "sender"：指定要进行差异表达分析的细胞类型为发送者细胞。这意味着函数会比较不同生态位中发送者细胞群里配体基因的表达差异 # 计算过程: # (1) 分组：依据 niches 里定义的发送者细胞群，把细胞划分成不同的组。

例如，在 "pEMT_High_niche" 和 "pEMT_Low_niche" 中分别确定对应的发送者细胞群 # (2) 统计检验：对每个配体基因，在不同组的发送者细胞中进行差异表达检验。常见的统计检验方法有 Wilcoxon 秩和检验、t 检验等。

这些检验会比较不同组中基因表达水平的均值或分布，以判断基因是否存在差异表达 # (3) 结果输出：函数会返回一个结果对象 DE_sender，其中包含每个配体基因的差异表达统计信息，例如 p 值、调整后的 p 值、平均表达差异等。

这些信息能够帮助你识别在不同生态位中发送者细胞群里显著差异表达的配体基因 DE_sender <- calculate_niche_de ( seurat_obj = seurat_obj %>% subset ( features = lr_network $ ligand %>% unique ), niches = niches, type = "sender" , assay_oi = assay_oi) # [1] "Calculate Sender DE between: myofibroblast_High and myofibroblast_Low" # [2] "Calculate Sender DE between: myofibroblast_High and Endothelial_Low" # [3] "Calculate Sender DE between: myofibroblast_High and CAF_Low" # [1] "Calculate Sender DE between: Endothelial_High and myofibroblast_Low" # [2] "Calculate Sender DE between: Endothelial_High and Endothelial_Low" # [3] "Calculate Sender DE between: Endothelial_High and CAF_Low" # [1] "Calculate Sender DE between: CAF_High and myofibroblast_Low" # [2] "Calculate Sender DE between: CAF_High and Endothelial_Low" # [3] "Calculate Sender DE between: CAF_High and CAF_Low" # [1] "Calculate Sender DE between: T.cell_High and myofibroblast_Low" # [2] "Calculate Sender DE between: T.cell_High and Endothelial_Low" # [3] "Calculate Sender DE between: T.cell_High and CAF_Low" # [1] "Calculate Sender DE between: Myeloid_High and myofibroblast_Low" # [2] "Calculate Sender DE between: Myeloid_High and Endothelial_Low" # [3] "Calculate Sender DE between: Myeloid_High and CAF_Low" # [1] "Calculate Sender DE between: myofibroblast_Low and myofibroblast_High" # [2] "Calculate Sender DE between: myofibroblast_Low and Endothelial_High" # [3] "Calculate Sender DE between: myofibroblast_Low and CAF_High" # [4] "Calculate Sender DE between: myofibroblast_Low and T.cell_High" # [5] "Calculate Sender DE between: myofibroblast_Low and Myeloid_High" # [1] "Calculate Sender DE between: Endothelial_Low and myofibroblast_High" # [2] "Calculate Sender DE between: Endothelial_Low and Endothelial_High" # [3] "Calculate Sender DE between: Endothelial_Low and CAF_High" # [4] "Calculate Sender DE between: Endothelial_Low and T.cell_High" # [5] "Calculate Sender DE between: Endothelial_Low and Myeloid_High" # [1] "Calculate Sender DE between: CAF_Low and myofibroblast_High" # [2] "Calculate Sender DE between: CAF_Low and Endothelial_High" # [3] "Calculate Sender DE between: CAF_Low and CAF_High" # [4] "Calculate Sender DE between: CAF_Low and T.cell_High" # [5] "Calculate Sender DE between: CAF_Low and Myeloid_High" head (DE_sender) # # A tibble: 6 × 8 # gene p_val avg_log2FC pct.1 pct.2 p_val_adj sender_other_niche # <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr> # 1 THBS4 0.0612 9.63 0.055 0 1 myofibroblas… myofibroblast_Low # 2 APOO 0.104 9.29 0.042 0 1 myofibroblas… myofibroblast_Low # 3 LRIG2 0.111 9.03 0.068 0.016 1 myofibroblas… myofibroblast_Low # 4 MFAP2 0.256 8.92 0.021 0 1 myofibroblas… myofibroblast_Low # 5 LGI1 0.0937 8.83 0.045 0 1 myofibroblas… myofibroblast_Low # 6 ICAM2 0.256 8.73 0.021 0 1 myofibroblas… myofibroblast_Low # receiver数量较少，算起来也比较快：DE_receiver <- calculate_niche_de ( seurat_obj = seurat_obj %>% subset ( features = lr_network $ receptor %>% unique ), niches = niches, type = "receiver" , assay_oi = assay_oi) # # A tibble: 1 × 2 # receiver_other_niche # <chr> <chr> # 1 Malignant_High Malignant_Low # [1] "Calculate receiver DE between: Malignant_High and Malignant_Low" # [1] "Calculate receiver DE between: Malignant_Low and Malignant_High" head (DE_receiver) # # A tibble: 6 × 8 # gene p_val avg_log2FC pct.1 pct.2 p_val_adj receiver_other_niche # <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr> # 1 CD3D 1.38e- 2 11.8 0.011 0 1 e+ 0 Maligna… Malignant_Low # 2 CTLA4 2.81e- 1 11.4 0.005 0.002 1 e+ 0 Maligna… Malignant_Low # 3 MMP1 2.06e-26 10.6 0.22 0.022 2.03e-23 Maligna… Malignant_Low # 4 ASGR2 6.21e- 1 10.1 0.001 0.002 1 e+ 0 Maligna… Malignant_Low # 5 MMP9 2.82e- 8 9.81 0.08 0.013 2.78e- 5 Maligna… Malignant_Low # 6 CLEC1… 2.20e- 1 9.63 0.003 0 1 e+ 0 Maligna… Malignant_Low # 将无穷大于无穷小的avg_log2FC转换为极值：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))) 去除表达比例不足10%的基因，并计算最小 avg_log2FC 与平均 avg_log2FC : expression_pct <- 0.10 # process_niche_de:遍历DE_table中的每个基因，检查其在相应细胞群体中的表达比例。

如果某个基因的表达比例低于 expression_pct，则将其从结果中过滤掉。这有助于去除那些可能是噪声或低质量的基因表达信息 # 函数会根据 niches 中定义的细胞生态位信息，将差异表达结果进行分类。

例如，对于发送者细胞，会将结果按照不同的发送者细胞群和生态位进行分组，以便后续分析不同生态位中发送者细胞的基因表达差异 DE_sender_processed <- process_niche_de ( DE_table = DE_sender, niches = niches, expression_pct = expression_pct, type = "sender" ) head (DE_sender_processed) # # A tibble: 6 × 9 # niche sender gene mean_avg_log2FC min_avg_log2FC mean_significant # <chr> <chr> <chr> <dbl> <dbl> <dbl> # 1 pEMT_High_niche myofibr… COL9… 8.85 8.46 0 # 2 pEMT_High_niche myofibr… CALCA 7.45 7.06 0 # 3 pEMT_High_niche myofibr… S100B 7.86 6.97 0 # 4 pEMT_High_niche myofibr… COL2… 7.09 6.70 0 # 5 pEMT_High_niche myofibr… PTPR… 6.83 6.30 0 # 6 pEMT_High_niche myofibr… PRL 6.26 5.87 0 # # ℹ 3 more variables: mean_present <dbl>, mean_score <dbl>, min_score <dbl> DE_receiver_processed <- process_niche_de ( DE_table = DE_receiver, niches = niches, expression_pct = expression_pct, type = "receiver" ) head (DE_receiver_processed) # # A tibble: 6 × 9 # niche receiver gene mean_avg_log2FC min_avg_log2FC mean_significant # <chr> <chr> <chr> <dbl> <dbl> <dbl> # 1 pEMT_High_niche Maligna… CD3D 11.8 11.8 0 # 2 pEMT_High_niche Maligna… CTLA4 11.4 11.4 0 # 3 pEMT_High_niche Maligna… MMP1 10.6 10.6 1 # 4 pEMT_High_niche Maligna… ASGR2 10.1 10.1 0 # 5 pEMT_High_niche Maligna… MMP9 9.81 9.81 1 # 6 pEMT_High_niche Maligna… CLEC… 9.63 9.63 0 # # ℹ 3 more variables: mean_present <dbl>, mean_score <dbl>, min_score <dbl> # 将经过处理的发送者细胞差异表达结果（DE_sender_processed）和接收者细胞差异表达结果（DE_receiver_processed）进行整合，并结合配体 - 受体网络信息（lr_network），同时根据指定的特异性得分方法（specificity_score_LR_pairs）来计算相关得分，最终得到整合后的差异表达结果 DE_sender_receiver # 举例：配体在 sender 的 log2FC 是 1.2，受体在 receiver 是 0.3 → 这个对的得分就是 0.3，说明受体的表达是限制因素 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) 3.4 定义感兴趣的基因集 与2.2类似，在开展正式细胞通讯计算前需要预先为每个 niche 确定一个感兴趣的基因集，这里主要研究目的是组间差异，因此这个基因集自然就是Malignant_High与Malignant_Low的差异基因。

# 设定阈值 lfc_cutoff <- 0.15 # 利用min_lfc筛选target：specificity_score_targets <- "min_lfc" # 计算下游目标基因（target genes）在 receiver 中的差异表达 # 使用完整的 Seurat（不需提前 subset） # 每个目标基因的DE信息，结果用于target-based下游分析 DE_receiver_targets <- calculate_niche_de_targets ( seurat_obj = seurat_obj, niches = niches, lfc_cutoff = lfc_cutoff, expression_pct <- expression_pct, assay_oi = assay_oi) # [1] "Calculate receiver DE between: Malignant_High and Malignant_Low" # [1] "Calculate receiver DE between: Malignant_Low and Malignant_High" head (DE_receiver_targets) # # A tibble: 6 × 8 # receiver_other_niche gene p_val avg_log2FC pct.1 pct.2 # <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> # 1 Malignant_High Malignant_Low C9orf152 1 e+ 0 0 NA # 2 Malignant_High Malignant_Low RPS11 1.12e-15 -1.95 0.972 0.995 # 3 Malignant_High Malignant_Low ELMO2 5.23e- 2 -0.271 0.263 0.31 # 4 Malignant_High Malignant_Low CREB3L1 1 e+ 0 0 NA # 5 Malignant_High Malignant_Low PNMA1 1.01e- 9 -1.00 0.199 0.339 # 6 Malignant_High Malignant_Low MMP2 8.78e-30 2.87 0.369 0.106 # # ℹ 1 more variable: p_val_adj <dbl> # 过滤并计算target score # 依据在 receiver 群体中的表达比例（expression_pct）过滤低表达目标基因 # 根据选择的 specificity_score 方法（如 "min_lfc"）对 target gene 在各receiver群体中差异表达的 特异性进行打分 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) head (DE_receiver_processed_targets) # # A tibble: 6 × 6 # niche receiver target_score target_significant target_present # <chr> <chr> <chr> <dbl> <dbl> <dbl> # 1 pEMT_Low_niche Maligna… CT45A3 15.0 1 1 # 2 pEMT_Low_niche Maligna… CT45A1 14.9 1 1 # 3 pEMT_High_niche Maligna… MSLN 13.8 1 1 # 4 pEMT_Low_niche Maligna… NR0B1 13.3 1 1 # 5 pEMT_High_niche Maligna… SPINK6 13.1 1 1 # 6 pEMT_Low_niche Maligna… DSCR8 13.0 1 1 # 确定背景基因 # 合并所有 receiver 细胞群的 DE 结果 background <- DE_receiver_processed_targets %>% pull (target) %>% unique # niche1的target基因集 geneset_niche1 <- DE_receiver_processed_targets %>% filter (receiver == niches[[ 1 ]] $ receiver & target_score >= lfc_cutoff & target_significant == 1 & target_present == 1 ) %>% pull (target) %>% unique # niche2的target基因集 geneset_niche2 <- DE_receiver_processed_targets %>% filter (receiver == niches[[ 2 ]] $ receiver & target_score >= lfc_cutoff & target_significant == 1 & target_present == 1 ) %>% pull (target) %>% unique # 如果有很多基因在这个过程中被drop out，请检查自己的基因名版本 length (geneset_niche1 %>% setdiff ( rownames (ligand_target_matrix))) # [1] 97 length (geneset_niche2 %>% setdiff ( rownames (ligand_target_matrix))) # [1] 514 # 查看最终的target基因数目 length (geneset_niche1) # [1] 1366 length (geneset_niche2) # [1] 4992 # 查看最终的target数据类型 # 其实最终得到的就是字符串，所以你也可以结合前期研究基础填一些自己真正感兴趣的基因 class (geneset_niche1) # [1] "character" # 查看前五个target基因ID geneset_niche1[ 1 : 5 ] # [1] "MSLN" "SPINK6" "NEFM" "LINC00162" "GPR128" geneset_niche2[ 1 : 5 ] # [1] "CT45A3" "CT45A1" "NR0B1" "DSCR8" "CTAG2" geneSet 的大小通常建议在20~1000之间，背景基因推荐至少5000个。

如果你的DE数量过多，可以设置较高的阈值对基因集进行过滤。

如果你有两个以上的 receiver cells/niches，lfc_cutoff 推荐为 0.15，若只有2个 receiver cells/niches，可以考虑用更高一些的的阈值，例如 0.25，若你的数据为 Smart-seq2 这可高测序深度的数据，可以使用更高的阈值，例如这里，高阈值会让数据的可信度增加：# 这里重新调整阈值，执行和前面相同的步骤 lfc_cutoff <- 0.75 specificity_score_targets <- "min_lfc" 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) # 背景基因 background <- DE_receiver_processed_targets %>% pull (target) %>% unique # niche1的target基因 geneset_niche1 <- DE_receiver_processed_targets %>% filter (receiver == niches[[ 1 ]] $ receiver & target_score >= lfc_cutoff & target_significant == 1 & target_present == 1 ) %>% pull (target) %>% unique # niche2的target基因 geneset_niche2 <- DE_receiver_processed_targets %>% filter (receiver == niches[[ 2 ]] $ receiver & target_score >= lfc_cutoff & target_significant == 1 & target_present == 1 ) %>% pull (target) %>% unique length (geneset_niche1 %>% setdiff ( rownames (ligand_target_matrix))) # [1] 67 length (geneset_niche2 %>% setdiff ( rownames (ligand_target_matrix))) # [1] 427 # 基于新的阈值，得到的每个niche的target的基因数目 length (geneset_niche1) # [1] 1047 length (geneset_niche2) # [1] 4018 3.5 计算配体活性并推断配受体对 # 定在后续分析中选取的前 n 个靶基因的数量 top_n_target <- 250 # 这部分代码创建了一个列表niche_geneset_list，其中包含两个子列表，分别对应两个细胞生态位：pEMT_High_niche 和 pEMT_Low_niche # 每个子列表包含三个元素：# "receiver"：指定该生态位下的接收细胞类型，分别从 niches 列表中获取。

# "geneset"：指定该生态位下的靶基因集，分别为 geneset_niche1 和 geneset_niche2。# "background"：指定背景基因集，这里两个生态位使用相同的背景基因集 background。

niche_geneset_list <- list ( "pEMT_High_niche" <- list ( "receiver" = niches[[ 1 ]] $ receiver, "geneset" = geneset_niche1, "background" = background), "pEMT_Low_niche" <- list ( "receiver" = niches[[ 2 ]] $ receiver, "geneset" = geneset_niche2 , "background" = background) ) # get_ligand_activities_targets 函数来计算配体的活性 # 时间较长，大约十分钟 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) # [1] "Calculate Ligand activities for: Malignant_High" # [1] "Calculate Ligand activities for: Malignant_Low" head (ligand_activities_targets) # # A tibble: 6 × 8 # ligand activity target ligand_target_weight receiver activity_normalized # <chr> <dbl> <chr> <dbl> <chr> <dbl> # 1 A2M 0.103 ANGPTL4 0.00747 Malignant_Hi… 0.275 # 2 A2M 0.103 APP 0.0112 Malignant_Hi… 0.275 # 3 A2M 0.103 ASS1 0.00892 Malignant_Hi… 0.275 # 4 A2M 0.103 ATF3 0.0107 Malignant_Hi… 0.275 # 5 A2M 0.103 BIRC3 0.00988 Malignant_Hi… 0.275 # 6 A2M 0.103 C3 0.00717 Malignant_Hi… 0.275 # # ℹ 2 more variables: scaled_activity_normalized <dbl>, scaled_activity <dbl> 3.6 计算配体、受体、target的表达量 NA。

features_oi <- union (lr_network $ ligand, lr_network $ receptor) %>% union (ligand_activities_targets $ target) %>% setdiff ( NA ) # 对 niches 中所有涉及的细胞类型进行 DotPlot绘图 dotplot <- suppressWarnings ( Seurat :: DotPlot (seurat_obj %>% subset ( idents = niches %>% unlist %>% unique ), features = features_oi, assay = assay_oi)) # 获取每个基因在每种 celltype 的平均表达（avg.exp）、scaled 表达（avg.exp.scaled）和表达细胞百分比（pct.exp） 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)) head (exprs_tbl) # # A tibble: 6 × 5 # celltype gene expression_scaled fraction # <fct> <chr> <dbl> <dbl> <dbl> # 1 myofibroblast_High A2M 6208. 1.04 0.966 # 2 myofibroblast_Low A2M 4346. 0.916 0.885 # 3 CAF_High A2M 882. 0.374 0.588 # 4 Malignant_Low A2M 9.10 -1.15 0.0419 # 5 Malignant_High A2M 21.4 -0.875 0.0677 # 6 CAF_Low A2M 349. 0.0600 0.510 # 将大表 exprs_tbl 拆分为三个小表（ligand、receptor 和 target） # exprs_tbl_ligand：每个 sender 细胞中的 ligand 表达 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：每个 receiver 细胞中的 receptor 表达 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：每个 receiver 细胞中的 target 基因表达 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) # scale_quantile_adapted：自定义分位数标准化函数，用于让不同基因/细胞类型之间的表达值可以比较 # mutate_cond(...)：如果表达比例 ≥ expression_pct（比如 0.10），就截断为 expression_pct，避免高比例值主导评分 # 最终保留标准化后的表达值和表达比例，作为后续网络过滤、配体选择等依据 exprs_tbl_ligand <- exprs_tbl_ligand %>% mutate ( scaled_ligand_expression_scaled = scale_quantile_adapted (ligand_expression_scaled)) %>% ligand_fraction_adapted，初始值等于 ligand_fraction mutate ( ligand_fraction_adapted = ligand_fraction) %>% ligand_fraction 大于等于 expression_pct，则将 ligand_fraction_adapted 的值设为 expression_pct mutate_cond (ligand_fraction >= expression_pct, ligand_fraction_adapted = expression_pct) %>% ligand_fraction_adapted 列进行分位数标准化，得到新的列 scaled_ligand_fraction_adapted mutate ( scaled_ligand_fraction_adapted = scale_quantile_adapted (ligand_fraction_adapted)) # receptor 的处理是同样逻辑 # scale_quantile_adapted：自定义分位数标准化函数，用于让不同基因/细胞类型之间的表达值可以比较 # mutate_cond(...)：如果表达比例 ≥ expression_pct（比如 0.10），就截断为 expression_pct，避免高比例值主导评分 # 最终保留标准化后的表达值和表达比例，作为后续网络过滤、配体选择等依据 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)) 3.7 通过受体表达量评估配受体对互作强度 # 一是将不同的数据框进行连接操作，整合配体、受体、细胞间通讯的差异表达等信息 # 二是对整合后的数据进行分组计算和排序，得出配体 - 受体 - 接收细胞组合下的综合得分 # 使用 inner_join 函数将lr_network与exprs_tbl_ligand按ligand列进行内连接。

内连接意味着只保留两个数据框中ligand列值相同的行，这样就把配体相关的表达信息添加到 lr_network 中 exprs_sender_receiver <- lr_network %>% inner_join (exprs_tbl_ligand, by = c ( "ligand" )) %>% # 接着将上一步结果与exprs_tbl_receptor按receptor列进行内连接，把受体相关的表达信息也整合进来 inner_join (exprs_tbl_receptor, by = c ( "receptor" )) %>% # 最后将前面的结果与DE_sender_receiver按niche、sender和receiver列进行内连接，并且只保留 DE_sender_receiver 中这三列的唯一行 inner_join (DE_sender_receiver %>% distinct (niche, sender, receiver)) # 按照ligand和receiver列对exprs_sender_receiver数据框进行分组，后续的计算将在每个分组内独立进行 ligand_scaled_receptor_expression_fraction_df <- exprs_sender_receiver %>% group_by (ligand, receiver) %>% # 使用dense_rank函数分别对每个分组内的receptor_expression和receptor_fraction列进行排名，得到新的列rank_receptor_expression和rank_receptor_fraction。

dense_rank函数会为相同的值赋予相同的排名，且排名是连续的 mutate ( rank_receptor_expression = dense_rank (receptor_expression), rank_receptor_fraction = dense_rank (receptor_fraction)) %>% # 计算每个分组内的综合得分ligand_scaled_receptor_expression_fraction。

具体做法是先将rank_receptor_fraction和rank_receptor_expression分别除以各自分组内的最大值，将排名归一化到0到1之间，然后取两者的平均值乘以0.5 mutate ( ligand_scaled_receptor_expression_fraction = 0.5 * ( (rank_receptor_fraction / max (rank_receptor_fraction)) + ((rank_receptor_expression / max (rank_receptor_expression))) )) %>% # 去除 ligand、receptor、receiver 和 ligand_scaled_receptor_expression_fraction 列组合中的重复 distinct (ligand, receptor, receiver, ligand_scaled_receptor_expression_fraction) %>% distinct %>% # 取消分组操作，使数据框恢复到未分组的状态 ungroup 3.8 计算配受体对权重 参与权重计算的参数有：Ligand DE score 、 Scaled ligand expression 、 Ligand expression fraction 、 Ligand spatial DE score 、 Receptor DE score 、 Scaled receptor expression 、 Receptor expression fraction 、 Receptor expression strength 、 Receptor spatial DE score 、 Absolute ligand activity 、 Normalized ligand activity。

prioritizing_weights <- c ( "scaled_ligand_score" = 5 , "scaled_ligand_expression_scaled" = 1 , "ligand_fraction" = 1 , = 2, "scaled_receptor_score" = 0.5 , "scaled_receptor_expression_scaled" = 0.5 , "receptor_fraction" = 1 , "ligand_scaled_receptor_expression_fraction" = 1 , = 0, "scaled_activity" = 0 , "scaled_activity_normalized" = 1 ) # 将之前计算得到的多个数据对象（如差异表达结果、配体 - 受体表达分数数据框等）整合到一个名为 output 的列表中，方便后续在自定义函数中调用这些数据 # DE_sender_receiver 发送细胞和接收细胞的差异表达分析整合结果 # ligand_scaled_receptor_expression_fraction_df 配体 - 受体组合的综合表达情况 # ligand_activities_targets 配体对靶基因的活性信息 # DE_receiver_processed_targets 接收细胞靶基因的差异表达结果 # exprs_tbl_ligand 配体在不同发送细胞中的表达水平、标准化后的表达水平以及表达比例等数据 # exprs_tbl_receptor 不同接收细胞中的表达水平、标准化后的表达水平以及表达比例等 # exprs_tbl_target 靶基因在不同接收细胞中的表达水平、标准化后的表达水平以及表达比例等 output <- list ( DE_sender_receiver = DE_sender_receiver, ligand_scaled_receptor_expression_fraction_df = ligand_scaled_receptor_expression_fraction_df, 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) # 修改为：my_get_prioritization_tables <- function (output_nichenet_analysis, prioritizing_weights) { # 数据连接 combined_information = output_nichenet_analysis $ DE_sender_receiver %>% inner_join (output_nichenet_analysis $ ligand_scaled_receptor_expression_fraction_df, by = c ( "receiver" , "ligand" , "receptor" )) %>% inner_join (output_nichenet_analysis $ ligand_activities_targets, by = c ( "receiver" , "ligand" )) %>% left_join (output_nichenet_analysis $ DE_receiver_processed_targets, by = c ( "niche" , "receiver" , "target" )) %>% inner_join (output_nichenet_analysis $ exprs_tbl_ligand, by = c ( "sender" , "ligand" )) %>% inner_join (output_nichenet_analysis $ exprs_tbl_receptor, by = c ( "receiver" , "receptor" )) %>% left_join (output_nichenet_analysis $ exprs_tbl_target, by = c ( "receiver" , "target" )) # 生成新列 combined_information = combined_information %>% mutate ( ligand_receptor = paste (ligand, receptor, sep = "--" )) # # 选择需要的列并去重 combined_information = combined_information %>% select (niche, receiver, sender, ligand_receptor, ligand, receptor, target, ligand_score, ligand_significant, ligand_present, ligand_expression, ligand_expression_scaled, ligand_fraction, receptor_score, receptor_significant, receptor_present, receptor_expression, receptor_expression_scaled, receptor_fraction, ligand_scaled_receptor_expression_fraction, avg_score_ligand_receptor, target_score, target_significant, target_present, target_expression, target_expression_scaled, target_fraction, ligand_target_weight, activity, activity_normalized, scaled_ligand_score, scaled_ligand_expression_scaled, scaled_receptor_score, scaled_receptor_expression_scaled, scaled_avg_score_ligand_receptor, scaled_ligand_fraction_adapted, scaled_receptor_fraction_adapted, scaled_activity, scaled_activity_normalized) %>% distinct # 计算优先级分数 combined_information_prioritized = combined_information %>% dplyr :: mutate ( prioritization_score = ((prioritizing_weights[ "scaled_ligand_score" ] * scaled_ligand_score) + (prioritizing_weights[ "scaled_ligand_expression_scaled" ] * scaled_ligand_expression_scaled) + (prioritizing_weights[ "scaled_receptor_score" ] * scaled_receptor_score) + (prioritizing_weights[ "scaled_receptor_expression_scaled" ] * scaled_receptor_expression_scaled) + (prioritizing_weights[ "ligand_scaled_receptor_expression_fraction" ] * ligand_scaled_receptor_expression_fraction) + (prioritizing_weights[ "scaled_activity" ] * scaled_activity) + (prioritizing_weights[ "scaled_activity_normalized" ] * scaled_activity_normalized) + (prioritizing_weights[ "ligand_fraction" ] * scaled_ligand_fraction_adapted) + (prioritizing_weights[ "receptor_fraction" ] * scaled_receptor_fraction_adapted) ) * ( 1 / length (prioritizing_weights))) %>% dplyr :: arrange ( - prioritization_score) # 生成配体 - 受体优先级表格 prioritization_tbl_ligand_receptor = combined_information_prioritized %>% select (niche, receiver, sender, ligand_receptor, ligand, receptor, ligand_score, ligand_significant, ligand_present, ligand_expression, ligand_expression_scaled, ligand_fraction, receptor_score, receptor_significant, receptor_present, receptor_expression, receptor_expression_scaled, receptor_fraction, ligand_scaled_receptor_expression_fraction, avg_score_ligand_receptor, activity, activity_normalized, scaled_ligand_score, scaled_ligand_expression_scaled, scaled_receptor_score, scaled_receptor_expression_scaled, scaled_avg_score_ligand_receptor, scaled_ligand_fraction_adapted, scaled_receptor_fraction_adapted, scaled_activity, scaled_activity_normalized, prioritization_score) %>% distinct # 生成配体 - 靶基因优先级表格 prioritization_tbl_ligand_target = combined_information_prioritized %>% select (niche, receiver, sender, ligand_receptor, ligand, receptor, target, target_score, target_significant, target_present, target_expression, target_expression_scaled, target_fraction, ligand_target_weight, activity, activity_normalized, scaled_activity, scaled_activity_normalized, prioritization_score) %>% distinct return ( list ( prioritization_tbl_ligand_receptor = prioritization_tbl_ligand_receptor, prioritization_tbl_ligand_target = prioritization_tbl_ligand_target)) } # 调用my_get_prioritization_tables函数，传入output和prioritizing_weights，得到包含两个优先级表格的 prioritization_tables列表 prioritization_tables <- my_get_prioritization_tables (output, prioritizing_weights) # 分别从 prioritization_tables 的两个表格中筛选出 receiver 为 niches 列表中第一个和第二个元素的 receiver 值的记录，并取前 10 行，用于查看特定接收细胞的高优先级配体 - 受体对和配体 - 靶基因对的信息 prioritization_tables $ prioritization_tbl_ligand_receptor %>% # 配体和受体优先级表格 filter (receiver == niches[[ 1 ]] $ receiver) %>% head ( 10 ) # # A tibble: 10 × 32 # niche receiver sender ligand_receptor ligand receptor ligand_score # <chr> <chr> <chr> <chr> <chr> <chr> <dbl> # 1 pEMT_High_niche Malignan… Myelo… C1QB--LRP1 C1QB LRP1 17.1 # 2 pEMT_High_niche Malignan… T.cel… GZMB--CHRM3 GZMB CHRM3 19.1 # 3 pEMT_High_niche Malignan… T.cel… TNF--TNFRSF21 14.7 # 4 pEMT_High_niche Malignan… T.cel… GZMB--MCL1 GZMB MCL1 19.1 # 5 pEMT_High_niche Malignan… Myelo… IL1RN--IL1R2 IL1RN IL1R2 16.1 # 6 pEMT_High_niche Malignan… Myelo… TNF--TNFRSF21 13.5 # 7 pEMT_High_niche Malignan… Myelo… S100B--ALCAM S100B ALCAM 14.1 # 8 pEMT_High_niche Malignan… Myelo… OSM--OSMR 14.5 # 9 pEMT_High_niche Malignan… T.cel… LTB--LTBR 16.1 # 10 pEMT_High_niche Malignan… Myelo… CCL22--DPP4 CCL22 DPP4 12.5 # # ℹ 25 more variables: ligand_significant <dbl>, ligand_present <dbl>, # # ligand_expression <dbl>, ligand_expression_scaled <dbl>, # # ligand_fraction <dbl>, receptor_score <dbl>, receptor_significant <dbl>, # # receptor_present <dbl>, receptor_expression <dbl>, # # receptor_expression_scaled <dbl>, receptor_fraction <dbl>, # # ligand_scaled_receptor_expression_fraction <dbl>, # # avg_score_ligand_receptor <dbl>, activity <dbl>, … prioritization_tables $ prioritization_tbl_ligand_target %>% # 配体和靶基因优先级表格 filter (receiver == niches[[ 1 ]] $ receiver) %>% head ( 10 ) # # A tibble: 10 × 19 # niche receiver sender ligand_receptor ligand receptor target_score # <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> # 1 pEMT_Hig… Maligna… Myelo… C1QB--LRP1 C1QB LRP1 APP 2.16 # 2 pEMT_Hig… Maligna… Myelo… C1QB--LRP1 C1QB LRP1 AREG 3.58 # 3 pEMT_Hig… Maligna… Myelo… C1QB--LRP1 C1QB LRP1 ASS1 1.10 # 4 pEMT_Hig… Maligna… Myelo… C1QB--LRP1 C1QB LRP1 ATF3 1.41 # 5 pEMT_Hig… Maligna… Myelo… C1QB--LRP1 C1QB LRP1 BIRC3 1.02 # 6 pEMT_Hig… Maligna… Myelo… C1QB--LRP1 C1QB LRP1 C3 1.44 # 7 pEMT_Hig… Maligna… Myelo… C1QB--LRP1 C1QB LRP1 CA2 9.94 # 8 pEMT_Hig… Maligna… Myelo… C1QB--LRP1 C1QB LRP1 CAV1 3.06 # 9 pEMT_Hig… Maligna… Myelo… C1QB--LRP1 C1QB LRP1 CCL20 2.54 # 10 pEMT_Hig… Maligna… Myelo… C1QB--LRP1 C1QB LRP1 CD14 3.36 # # ℹ 11 more variables: target_significant <dbl>, target_present <dbl>, # # target_expression <dbl>, target_expression_scaled <dbl>, # # target_fraction <dbl>, ligand_target_weight <dbl>, activity <dbl>, # # activity_normalized <dbl>, scaled_activity <dbl>, # # scaled_activity_normalized <dbl>, prioritization_score <dbl> prioritization_tables $ prioritization_tbl_ligand_receptor %>% # 配体和受体优先级表格 filter (receiver == niches[[ 2 ]] $ receiver) %>% head ( 10 ) # # A tibble: 10 × 32 # niche receiver sender ligand_receptor ligand receptor ligand_score # <chr> <chr> <chr> <chr> <chr> <chr> <dbl> # 1 pEMT_Low_niche Malignant… CAF_L… WIF1--RYK WIF1 RYK 13.2 # 2 pEMT_Low_niche Malignant… CAF_L… IGSF11--CLEC2L IGSF11 CLEC2L 8.09 # 3 pEMT_Low_niche Malignant… CAF_L… LRRN1--LRRC4 LRRN1 LRRC4 3.25 # 4 pEMT_Low_niche Malignant… myofi… IGFBPL1--DCC IGFBP… DCC 2.79 # 5 pEMT_Low_niche Malignant… CAF_L… LRFN5--PTPRS LRFN5 PTPRS 2.47 # 6 pEMT_Low_niche Malignant… CAF_L… IGSF10--IGSF11 IGSF10 IGSF11 3.63 # 7 pEMT_Low_niche Malignant… CAF_L… MMP13--F12 MMP13 F12 7.90 # 8 pEMT_Low_niche Malignant… CAF_L… SLIT2--SDC1 SLIT2 SDC1 2.35 # 9 pEMT_Low_niche Malignant… CAF_L… FGF14--SCN9A FGF14 SCN9A 1.51 # 10 pEMT_Low_niche Malignant… CAF_L… CHL1--TMEM132A CHL1 TMEM132A 2.71 # # ℹ 25 more variables: ligand_significant <dbl>, ligand_present <dbl>, # # ligand_expression <dbl>, ligand_expression_scaled <dbl>, # # ligand_fraction <dbl>, receptor_score <dbl>, receptor_significant <dbl>, # # receptor_present <dbl>, receptor_expression <dbl>, # # receptor_expression_scaled <dbl>, receptor_fraction <dbl>, # # ligand_scaled_receptor_expression_fraction <dbl>, # # avg_score_ligand_receptor <dbl>, activity <dbl>, … prioritization_tables $ prioritization_tbl_ligand_target %>% # 配体和靶基因优先级表格 filter (receiver == niches[[ 2 ]] $ receiver) %>% head ( 10 ) # # A tibble: 10 × 19 # niche receiver sender ligand_receptor ligand receptor target_score # <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> # 1 pEMT_Low… Maligna… CAF_L… WIF1--RYK WIF1 RYK AKAP1 1.53 # 2 pEMT_Low… Maligna… CAF_L… WIF1--RYK WIF1 RYK AXIN2 4.26 # 3 pEMT_Low… Maligna… CAF_L… WIF1--RYK WIF1 RYK BCL11A 2.04 # 4 pEMT_Low… Maligna… CAF_L… WIF1--RYK WIF1 RYK BCL6 1.20 # 5 pEMT_Low… Maligna… CAF_L… WIF1--RYK WIF1 RYK BDNF 4.98 # 6 pEMT_Low… Maligna… CAF_L… WIF1--RYK WIF1 RYK BIRC5 4.70 # 7 pEMT_Low… Maligna… CAF_L… WIF1--RYK WIF1 RYK BMP4 0.821 # 8 pEMT_Low… Maligna… CAF_L… WIF1--RYK WIF1 RYK BMP7 0.866 # 9 pEMT_Low… Maligna… CAF_L… WIF1--RYK WIF1 RYK BRCA1 1.62 # 10 pEMT_Low… Maligna… CAF_L… WIF1--RYK WIF1 RYK CBX5 2.25 # # ℹ 11 more variables: target_significant <dbl>, target_present <dbl>, # # target_expression <dbl>, target_expression_scaled <dbl>, # # target_fraction <dbl>, ligand_target_weight <dbl>, activity <dbl>, # # activity_normalized <dbl>, scaled_activity <dbl>, # # scaled_activity_normalized <dbl>, prioritization_score <dbl> 3.9 差异可视化 3.9.1 获取top配体与top受体 获取每个 niche 的top50配体 # 从 prioritization_tables 中的 prioritization_tbl_ligand_receptor 数据框里筛选出高优先级的配体、受体以及对应的细胞生态位信息。

下面分别对两段代码的功能进行详细解释：# prioritization_tables$prioritization_tbl_ligand_receptor 数据框中选取 niche（细胞生态位）、sender（发送细胞）、receiver（接收细胞）、ligand（配体）、receptor（受体）和 prioritization_score（优先级得分）这些列，生成一个新的数据子集 top_ligand_niche_df <- prioritization_tables $ prioritization_tbl_ligand_receptor %>% select (niche, sender, receiver, ligand, receptor, prioritization_score) %>% # 按照 ligand 列对数据进行分组，后续的操作会在每个配体组内独立执行 group_by (ligand) %>% # 每个配体组内，根据 prioritization_score 列的值选取优先级得分最高的那一行数据 top_n ( 1 , prioritization_score) %>% # 取消之前的分组操作，使数据恢复到未分组状态 ungroup %>% # 从筛选后的数据中选取 ligand、receptor 和 niche 这三列 select (ligand, receptor, niche) %>% # 将 niche 列重命名为 top_niche，最终得到 top_ligand_niche_df 数据框，该数据框记录了每个配体优先级最高的那一行对应的受体和细胞生态位信息 rename ( top_niche = niche) # 和第一段代码一样，从 prioritization_tables$prioritization_tbl_ligand_receptor 数据框中选取指定的列，生成新的数据子集 top_ligand_receptor_niche_df = prioritization_tables $ prioritization_tbl_ligand_receptor %>% select (niche, sender, receiver, ligand, receptor, prioritization_score) %>% # 按照 ligand 和 receptor 两列对数据进行分组，后续操作会在每个配体 - 受体组合组内独立执行 group_by (ligand, receptor) %>% # 每个配体 - 受体组合组内，根据 prioritization_score 列的值选取优先级得分最高的那一行数据 top_n ( 1 , prioritization_score) %>% # 取消分组操作，让数据回到未分组状态 ungroup %>% # 从筛选后的数据中选取 ligand、receptor 和 niche 这三列 select (ligand, receptor, niche) %>% # 将 niche 列重命名为 top_niche，最终得到 top_ligand_receptor_niche_df 数据框，该数据框记录了每个配体 - 受体组合优先级最高的那一行对应的细胞生态位信息 rename ( top_niche = niche) ligand_prioritized_tbl_oi <- prioritization_tables $ prioritization_tbl_ligand_receptor %>% select (niche, sender, receiver, ligand, prioritization_score) %>% # 按照 ligand（配体）和 niche（细胞生态位）对数据进行分组，后续的操作将在每个配体 - 细胞生态位组合的组内进行 group_by (ligand, niche) %>% # 在每个组内，根据 prioritization_score（优先级得分）选取得分最高的那一行数据，即找出每个配体在每个细胞生态位中的最高优先级记 top_n ( 1 , prioritization_score) %>% ungroup %>% distinct %>% # inner_join 函数将前面筛选得到的数据与top_ligand_niche_df数据框进行内连接。

连接的依据是两个数据框中共同的列，这样可以将 top_ligand_niche_df 中的信息整合到当前数据集中 inner_join (top_ligand_niche_df) %>% # 根据 niche 列和 top_niche 列的值进行筛选，只保留 niche 列的值与 top_niche 列的值相等的行。

这一步进一步筛选出符合特定条件的数据 filter (niche == top_niche) %>% # 按照 niche（细胞生态位）对数据进行分组 group_by (niche) %>% # 在每个细胞生态位组内，根据 prioritization_score（优先级得分）选取得分最高的前 50 行数据 top_n ( 50 , prioritization_score) %>% ungroup # 这段代码通过一系列的数据处理和筛选操作，对配体 - 受体的优先级数据进行了整理和分析，帮助快速获取每个细胞生态位中最重要的配体信息 每个配体获取top2的受体：# 对配体 - 受体的优先级数据进行进一步筛选和处理，以获取特定接收细胞（Malignant_High）相关的高优先级配体 - 受体对信息 receiver_oi <- "Malignant_High" # 从 ligand_prioritized_tbl_oi 数据框中筛选出 receiver 列的值等于 receiver_oi（即 "Malignant_High"）的行 # 从筛选后的数据中提取 ligand 列的数据，得到一个包含配体名称的向量 # 去除向量中的重复元素，得到一个包含唯一配体名称的向量 filtered_ligands filtered_ligands <- ligand_prioritized_tbl_oi %>% filter (receiver == receiver_oi) %>% pull (ligand) %>% unique prioritized_tbl_oi <- prioritization_tables $ prioritization_tbl_ligand_receptor %>% # 从 prioritization_tables 列表中的 prioritization_tbl_ligand_receptor 数据框中筛选出 ligand 列的值在 filtered_ligands 向量中的行 filter (ligand %in% filtered_ligands) %>% # 选取筛选后数据中的 niche（细胞生态位）、sender（发送细胞）、receiver（接收细胞）、ligand（配体）、receptor（受体）、ligand_receptor（配体 - 受体组合）和prioritization_score（优先级得分）这些列 select (niche, sender, receiver, ligand, receptor, ligand_receptor, prioritization_score) %>% # 去除重复的行，确保数据的唯一 distinct %>% # 将前面处理的数据与 top_ligand_receptor_niche_df 数据框进行内连接 inner_join (top_ligand_receptor_niche_df) %>% # 按照 ligand（配体）对连接后的数据进行分组，并在每个组内筛选出 receiver 列的值等于 receiver_oi（即 "Malignant_High"）的行 group_by (ligand) %>% filter (receiver == receiver_oi) %>% # 每个配体组内，根据 prioritization_score（优先级得分）选取得分最高的前 2 行数据 top_n ( 2 , prioritization_score) %>% ungroup # 上面这一步你也可以从"pEMT-low"的角度进行可视化：if (F){ receiver_oi = "Malignant_Low" filtered_ligands = ligand_prioritized_tbl_oi %>% filter (receiver == receiver_oi) %>% top_n ( 50 , 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 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 ) lfc_plot } 3.9.2 可视化 查看配体在 Sender 的组间差异，查看受体在 Receiver 中的组间差异：# 绘制配体-受体对的 log fold change（LFC）差异表达热图 # LFC 越高，表示该基因在该 niche 中差异表达越显著 # receiver_oi：receiver（接收细胞类型）"Malignant_High" # prioritized_tbl_oi根据优先级筛选出来的配体-受体表 # prioritization_tables$prioritization_tbl_ligand_receptor：完整的配体-受体优先级评分表 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 ) lfc_plot 查看配体表达量、配体活性以及target genes与配体的关联度 # 可视化在指定接收细胞类型（receiver_oi）中，配体-靶点信号活性和表达信息的整体概况 lfc_cutoff = 3 # 1.配体表达图, # 展示选定的发送细胞中，配体的表达强，来源于 output$exprs_tbl_ligand # 横轴为发送细胞类型，纵轴为配体，颜色表示表达水平 # 2.信号活性图（Ligand Activity Plot）：# 展示每个配体对于目标基因的调控活性（来源于 ligand_target_matrix 计算的活性分值，通常为 Pearson 分数或调控得分）。

# 通常表现为热图或者点图，纵轴为配体，横轴为 target 基因，颜色表示活性分值 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, # 对数倍数变化（Log2 Fold Change）的阈值，用于筛选具有显著表达差异 lfc_cutoff, ligand_target_matrix, plot_legend = FALSE , heights = NULL , widths = NULL ) exprs_activity_target_plot $ combined_plot 展示top20 # 从 ligand_prioritized_tbl_oi 数据框中筛选出receiver列等于receiver_oi（"Malignant_High"）的行 # 对筛选后的行，按照 prioritization_score（优先级得分）从高到低排序，选取前20行 # 提取这些行中的ligand列，并去除重复的配体名称，得到一个包含20个高优先级配体的向量filtered_ligands filtered_ligands <- ligand_prioritized_tbl_oi %>% filter (receiver == receiver_oi) %>% top_n ( 20 , prioritization_score) %>% pull (ligand) %>% unique # 从prioritization_tables$prioritization_tbl_ligand_receptor数据框中筛选出ligand列的值在filtered_ligands 向量中的行 # 选取筛选后数据中的 niche（细胞生态位）、sender（发送细胞）、receiver（接收细胞）、ligand（配体）、receptor（受体）、ligand_receptor（配体 - 受体组合）和 prioritization_score（优先级得分）这些列，并去除重复行 # 将处理后的数据与top_ligand_receptor_niche_df 数据框进行内连接，整合更多相关信息 # 按照 ligand 对连接后的数据进行分组，并在每个组内筛选出receiver列等于 receiver_oi 的行 # 对每个配体组内的行，按照 prioritization_score 从高到低排序，选取前2行，最后取消分组操作，得到一个新的数据框 prioritized_tbl_oi，包含了特定接收细胞相关的高优先级配体 - 受体对信息 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 ) # 整合和可视化细胞间通讯中与特定接收细胞相关的高优先级发送细胞中的配体表达、配体到下游基因的信号调控活性和接收细胞中靶点基因的表达 exprs_activity_target_plot $ combined_plot 配受体对的 circos 图片 # 选择优先级排名top15的配体 filtered_ligands <- ligand_prioritized_tbl_oi %>% filter (receiver == receiver_oi) %>% top_n ( 15 , prioritization_score) %>% pull (ligand) %>% unique # 保留每个配体在当前 receiver（receiver_oi）下得分最高的两个配体-受体对 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 # 设置颜色 # brewer.pal(n, name = "Spectral")：从 "Spectral" 调色板中选取 n 种颜色 # magrittr::set_names：将这些颜色与sender列的不同值一一对应，得到一个命名向量colors_sender，用于在绘图时为不同的发送细胞类型分配颜色 colors_sender <- brewer.pal ( n = prioritized_tbl_oi $ sender %>% unique %>% sort %>% length , name = "Spectral" ) %>% magrittr :: set_names (prioritized_tbl_oi $ sender %>% unique %>% sort ) # receiver统一用lavender色(淡紫色)，这里的receiver为Malignant_High colors_receiver <- c ( "lavender" ) %>% magrittr :: set_names (prioritized_tbl_oi $ receiver %>% unique %>% sort ) # 这个函数把上面筛选出的配体-受体对绘制成一个 circos plot，显示它们的连接关系，并使用 sender/receiver 的颜色编码 circos_output <- make_circos_lr (prioritized_tbl_oi, colors_sender, colors_receiver) Sys.time # [1] "2025-05-15 11:29:40 CST" save.image ( "result_nichenet/nichenet.rdata" )

六、cellcall CellCall是一个通过整合细胞内和细胞间信号来推断细胞间通讯网络和内部调节信号的工具包。(1) CellCall基于KEGG 通路收集配体-受体-转录因子(L-R-TF)数据集。(2) 根据L-R-TF相互作用的先验知识，CellCall通过结合配体/受体的表达和某些L-R对的下游TF活动来推断细胞间通讯。(3) CellCall嵌入了一种通路活性分析方法来识别某些细胞类型之间的通信所涉及的关键通路。

(4) CellCall提供了一套丰富的可视化选项（Circos plot、Sankey plot、bubble plot、ridge plot等）来直观地呈现分析结果。

1、软件安装和加载 # devtools::install_github("ShellyCoder/cellcall") library (cellcall) library (stringr) library (Seurat) library (pheatmap) library (ggplot2) 2、数据准备 2.1 软件自带测试数据 Sys.time # [1] "2025-05-15 11:34:22 CST" rm ( list = ls ) gc # used (Mb) gc trigger (Mb) max used (Mb) # Ncells 9024234 482.0 13029127 695.9 13029127 695.9 # Vcells 100483519 766.7 677141600 5166.2 845333542 6449.4 setwd ( "/data/07_scRNA-seq_cell_communication/" ) # 加载测试数据集 f.tmp <- system.file ( "extdata" , "example_Data.Rdata" , package= "cellcall" ) load (f.tmp) # cellcall需要的分析数据是一个行名为基因，列名为细胞类型的表达矩阵 # 查看数据维度 dim (in.content) # [1] 35135 366 # 查看数据前四行和前四列 in.content[ 1 : 4 , 1 : 4 ] # 1_ST 2_ST 3_ST 4_ST # TSPAN6 2.278 2.031 0.000 0.000 # TNMD 0.000 0.000 0.000 0.000 # DPM1 9.111 0.000 24.846 0.553 # SCYL3 0.000 0.000 0.000 0.000 # 查看细胞种类 table ( str_split ( colnames (in.content), "_" , simplify = T)[, 2 ]) # # SPGed SPGing SSC ST # 42 122 159 43 2.2 从seurat对象开始 seuratObj <- readRDS ( "data/seuratObj.rds" ) # 如遇报错"no slot of name "images" for this object of class "Seurat""，可以使用UpdateSeuratObject更细seurat对象 seuratObj <- UpdateSeuratObject (seuratObj) seuratObj # An object of class Seurat # 13541 features across 5027 samples within 1 assay # Active assay: RNA (13541 features, 1575 variable features) # 3 layers present: counts, data, scale.data # 4 dimensional reductions calculated: cca, cca.aligned, tsne, pca # 抽取部分细胞用于后续分析 seuratObj <- subset (seuratObj, downsample= 300 ) # 查看数据基本信息 # 注释信息 head (seuratObj @ meta.data) # nGene nUMI orig.ident aggregate res.0.6 celltype nCount_RNA # W380404 612 946 LN_SS 3 Treg 938 # W380420 734 1367 LN_SS 2 B 1356 # W380441 965 1912 LN_SS 1 CD8 T 1907 # W380455 541 987 LN_SS 0 CD4 T 984 # W380460 1156 2887 LN_SS 0 CD4 T 2881 # W380467 935 1994 LN_SS 1 CD8 T 1988 # nFeature_RNA # W380404 606 # W380420 723 # W380441 961 # W380455 539 # W380460 1151 # W380467 930 # cellcall需要的分析数据是一个行名为基因，列名为细胞类型的表达矩阵 # 从Seurat对象提取表达矩阵 in.content.new <- as.data.frame ( GetAssayData (seuratObj, slot= "count" )) # 修改列名,组成为细胞barcode_celltype colnames (in.content.new) <- paste0 ( colnames (in.content.new), "_" ,seuratObj @ meta.data $ celltype) in.content.new[ 1 : 4 , 1 : 4 ] # W380404_Treg W380420_B W380441_CD8 T W380455_CD4 T # 0610005C13Rik 0 0 0 0 # 0610007C21Rik 0 0 1 0 # 0610007L01Rik 0 0 0 0 # 0610007P08Rik 0 0 0 0 3、创建cellcall对象 3.1 软件自带测试数据 # 输入的基因表达数据，通常是一个数据框或矩阵形式，每一行表示一个基因，每一列表示一个细胞 mt <- CreateNichConObject ( data= in.content, # 一个基因至少需要出现在 3 个细胞中表达 min.feature = 3 , # 列名由索引和单元格类型组成。

CreateNichConObject中设置names.field=2和拆分字符分隔符为"_" names.field = 2 , names.delim = "_" , # 指定数据的来源是TPM，即每百万转录本的标准化基因表达量 source = "TPM" , # 设置标准化的缩放因子。

TPM通过总转录本数进行标准化，10^6是一个常用的标准化系数，用来将原始读数标准化到每百万单位 scale.factor = 10 ^ 6 , # 指定物种为人类 Org = "Homo sapiens" , # 设定当前项目的名称 project = "Microenvironment" ) 3.2 从seurat对象开始 mt.new <- CreateNichConObject ( data= in.content.new, # 一个基因至少需要出现在 3 个细胞中表达 min.feature = 3 , # 列名由索引和单元格类型组成。

CreateNichConObject中设置names.field=2和拆分字符分隔符为"_" names.field = 2 , names.delim = "_" , # 指定数据的来源是Seurat对象的UMI source = "UMI" , # 设置标准化的缩放因子。

TPM通过总转录本数进行标准化，10^6是一个常用的标准化系数，用来将原始读数标准化到每百万单位 scale.factor = 10 ^ 6 , # 指定物种为小鼠 Org = "Mus musculus" , # 设定当前项目的名称 project = "Seurat" ) 4、推断细胞 - 细胞通信得分 4.1 软件自带测试数据 # 约35min mt <- TransCommuProfile ( # 输入的对象 object = mt, # 使用spearson筛选TF的靶基因，默认值为0.05，意味着只有p值小于0.05的相关性才会被认为是显著的 pValueCor = 0.05 , # 设定相关性值的阈值。

只有大于该阈值的基因之间的相关性才会被考虑。这可以用来去除较小的、可能不太重要的相关性 CorValue = 0.1 , # 指定每个配体 - 受体对保留的目标基因数量。若topTargetCor = 1，则每个配体 - 受体对只保留相关性最强的 1 个目标基因 topTargetCor= 1 , # 这是经过多重检验校正后的p值阈值。在进行大量的相关性检验时，会出现多重比较问题，因此需要对 p 值进行校正。

p.adjust就是用来筛选校正后显著的结果。例如p.adjust=0.05，表示只有校正后的p值小于0.05 的结果才会被保留 p.adjust = 0.05 , # 指定在计算过程中使用的统计量类型。

use.type = "median"表示使用中位数来进行计算 use.type= "median" , # 一种细胞类型中代表该细胞类型的基因表达百分比 probs = 0.9 , # "weighted"、"max"、"mean", 其中 "weighted" 是默认值。

选择合适的方法来评估给定配体 - 受体关系的所有规则的下游配体 - 受体激活 method= "weighted" , # 这是一个布尔型参数，IS_core = TRUE表示将高置信度的核心参考LR数据纳入分析 IS_core = TRUE , # 只当分析的物种为人类 Org = 'Homo sapiens' ) # Ligand_ID Receptor_ID TF_ID Pathway Ligand_Symbol # 1 100506658 2626 hsa04530_3,hsa04530_5 OCLN # 2 100506658 8531 hsa04530_2 OCLN # 3 10344 12309 hsa04062_5 CCL26 # 4 10344 1230 3551 hsa04062_5 CCL26 # Receptor_Symbol TF_Symbol # 1 OCLN GATA4 # 2 OCLN YBX3 # 3 CCR1 FOXO3 # 4 CCR1 IKBKB # [1] "step1: compute means of gene" # [1] "ST" "SSC" "SPGing" "SPGed" # [1] "ST" # [1] "SSC" # [1] "SPGing" # [1] "SPGed" # [1] "step2: filter tf-gene with correlation, then scoregulons" # [1] "ST" # [1] 26 # [1] "SSC" # [1] 19 # [1] "SPGing" # [1] 1 # [1] "SPGed" # [1] 1 # [1] "ST" # [1] 26 # [1] "SSC" # [1] 19 # [1] "SPGing" # [1] 1 # [1] "SPGed" # [1] 1 # [1] "step3: get distance between receptor and tf in pathway" # [1] "step4: score downstream activation of ligand-receptor all regulons of given ligand-receptor relation (weighted, max, or mean) # [1] "step5: softmax for ligand" # [1] "step6: score ligand-receptor relation (weighted, max, or mean) # 创建保存结果的文件夹 dir.create ( "result_cellcall/" ) # 保存cellcall对象 saveRDS (mt, "result_cellcall/cellcall.rds" ) # 下次可以直接读取保存的rds文件继续后续流程的分析 # mt <- readRDS("result_cellcall/cellcall.rds") 4.2 从seurat对象开始 mt.new <- TransCommuProfile ( # 输入的对象 object = mt.new, # 使用spearson筛选TF的靶基因，默认值为0.05，意味着只有p值小于0.05的相关性才会被认为是显著的 pValueCor = 0.05 , # 设定相关性值的阈值。

只有大于该阈值的基因之间的相关性才会被考虑。这可以用来去除较小的、可能不太重要的相关性 CorValue = 0.1 , # 指定每个配体 - 受体对保留的目标基因数量。若topTargetCor = 1，则每个配体 - 受体对只保留相关性最强的 1 个目标基因 topTargetCor= 1 , # 这是经过多重检验校正后的p值阈值。在进行大量的相关性检验时，会出现多重比较问题，因此需要对 p 值进行校正。

p.adjust就是用来筛选校正后显著的结果。例如p.adjust=0.05，表示只有校正后的p值小于0.05 的结果才会被保留 p.adjust = 0.05 , # 指定在计算过程中使用的统计量类型。

use.type = "median"表示使用中位数来进行计算 use.type= "median" , # 一种细胞类型中代表该细胞类型的基因表达百分比 probs = 0.9 , # "weighted"、"max"、"mean", 其中 "weighted" 是默认值。

选择合适的方法来评估给定配体 - 受体关系的所有规则的下游配体 - 受体激活 method= "weighted" , # 这是一个布尔型参数，IS_core = TRUE表示将高置信度的核心参考LR数据纳入分析 IS_core = TRUE , # 只当分析的物种为小鼠 Org = 'Mus musculus' ) # 创建保存结果的文件夹 dir.create ( "result_cellcall/" ) # 保存cellcall对象 saveRDS (mt.new, "result_cellcall/cellcall.new.rds" ) # 下次可以直接读取保存的rds文件继续后续流程的分析 # mt.new <- readRDS("result_cellcall/cellcall.new.rds") 5、通路活动分析 从第五步开始后续的分析步骤，可以根据第四步生成的结果，进行计算和展示，这里以第四步生成的mt对象为例，也可以使用生成的mt.new对象分析，大家根据自己定义的变量名称灵活调整就可以了，注意物种也要做对应的调整 # CellCall 嵌入了一种通路活性分析方法，帮助探索特定细胞之间通信所涉及的主要通路 # 其中mt为Cellcall S4 对象，函数 CreateNichConObject 和 TransCommuProfile 的结果 # expr_l_r_log2_scale通信评分的数据框，其中行名称为配体受体，列名称为 cellA-cellB n <- mt @ data $ expr_l_r_log2_scale # getHyperPathway获取每个配体和受体互作细胞对富集的通路 pathway.hyper.list <- lapply ( colnames (n), function (i){ print (i) tmp <- tryCatch ({ getHyperPathway ( # 具有行 LR 和列 cellA-cellB 的通信评分数据框 data = n, # 修改对象名称 object = mt, # 探索发送单元A和接收单元B之间的LR, 例如 “A-B” cella_cellb = i, # 修改物种名称 Org = "Homo sapiens" ) }, error = function (e) { # 如果没有找到满足条件的 ligand-receptor 对就跳过 message ( sprintf ( "skip" )) return ( NULL ) }) return (tmp) }) # [1] "ST-ST" # [1] "ST-SSC" # [1] "ST-SPGing" # [1] "ST-SPGed" # [1] "SSC-ST" # [1] "SSC-SSC" # [1] "SSC-SPGing" # [1] "SSC-SPGed" # [1] "SPGing-ST" # [1] "SPGing-SSC" # [1] "SPGing-SPGing" # [1] "SPGing-SPGed" # [1] "SPGed-ST" # [1] "SPGed-SSC" # [1] "SPGed-SPGing" # [1] "SPGed-SPGed" # 对于路径活动分析，采用气泡图来呈现分析结果。

# 使用 getForBubble 函数合并数据 myPub.df <- getForBubble (pathway.hyper.list, cella_cellb= colnames (n)) head (myPub.df) # pathway cc pvalue NES # 1 Notch signaling pathway ST-ST 1.399286e-14 1.140034 # 2 Breast cancer ST-ST 4.477741e-11 1.067615 # 3 Th1 and Th2 cell differentiation ST-ST 7.733967e-06 1.674364 # 4 PI3K-Akt signaling pathway ST-ST 4.558940e-03 0.000000 # 5 Jak-STAT signaling pathway ST-ST 6.383580e-01 0.000000 # 6 Proteoglycans in cancer ST-ST 7.287498e-01 0.000000 # 使用 plotBubble 绘制气泡图 p <- plotBubble (myPub.df) + theme ( axis.text.y = element_text ( color = "black" )) + coord_flip # 横纵坐标反转 p 6、可视化 6.1 圈图 # 创建一个数据框 cell_color，用于存储不同细胞类型对应的颜色信息，这里提供的颜色数目和细胞类型数目要一致 cell_color <- data.frame ( color= c ( ' , ' , ' , ' ), # 确保数据框中的 color 列被存储为字符类型，而不是因子类型 stringsAsFactors = FALSE ) # 设置数据框行名称 rownames (cell_color) <- c ( "SSC" , "SPGing" , "SPGed" , "ST" ) cell_color # color # SSC # SPGing # SPGed # ST ViewInterCircos ( # 绘图对象 object = mt, # 字体大小 font = 2 , # 细胞类型颜色 cellColor = cell_color, # 配体和受体颜色 lrColor = c ( " , " ), # 设置连线箭头的类型为 "big.arrow"，即大箭头样式 arr.type = "big.arrow" , # 设置箭头的长度 arr.length = 0.04 , # 外侧轨道高度 trackhight1 = 0.05 , # 用特定槽的数据绘制图表 slot= "expr_l_r_log2_scale" , # 表示连线的颜色根据发送方细胞类型来确定 linkcolor.from.sender = TRUE , # 如果 linkcolor.from.sender 为 TRUE，则此参数通常设置为 NULL linkcolor = NULL , # 设置环形图中不同细胞类型之间的间隔角度 gap.degree = 2 , # 指定细胞类型在环形图中的排列顺序 order.vector= c ( 'ST' , "SSC" , "SPGing" , "SPGed" ), # 设置第二个轨道高度 trackhight2 = 0.032 , # 设置第二个轨道的边距，是一个包含两个值的向量，分别表示上下边距 track.margin2 = c ( 0.01 , 0.12 ), # 表示不使用自定义模式，使用函数的默认设置进行绘图 DIY = FALSE ) # [1] 0 # [1] 0 # [1] 0 # [1] 0 # [1] 0.02067826 # [1] 0.02234305 # [1] 0.06404499 # [1] 0.09506311 # [1] 0.2732887 # [1] 0.3378797 # [1] 0.498962 # [1] 0.5228232 # [1] 0.5467831 # [1] 0.6220295 # [1] 0.8290871 # [1] 1 # 第二次调用 ViewInterCircos 函数 ViewInterCircos ( # 直接传入 mt 对象中 data 槽位下的 expr_l_r_log2_scale 矩阵作为输入数据 object = mt @ data $ expr_l_r_log2_scale, # 中间部分参数和前面第一次调用设置意思一致 font = 2 , cellColor = cell_color, lrColor = c ( " , " ), arr.type = "big.arrow" , arr.length = 0.04 , trackhight1 = 0.05 , slot= "expr_l_r_log2_scale" , linkcolor.from.sender = TRUE , linkcolor = NULL , gap.degree = 2 , order.vector= c ( 'ST' , "SSC" , "SPGing" , "SPGed" ), trackhight2 = 0.032 , track.margin2 = c ( 0.01 , 0.12 ), # 启用自定义模式，允许用户根据自己的需求对绘图进行更多的自定义设置 # 逻辑值，如果为 TRUE, 参数对象应为数据框，并设置 slot="expr_l_r_log2_scale"。

否则，对象应为 cellcall对象 DIY = T) # [1] 0 # [1] 0 # [1] 0 # [1] 0 # [1] 0.02067826 # [1] 0.02234305 # [1] 0.06404499 # [1] 0.09506311 # [1] 0.2732887 # [1] 0.3378797 # [1] 0.498962 # [1] 0.5228232 # [1] 0.5467831 # [1] 0.6220295 # [1] 0.8290871 # [1] 1 6.2 热图 # 矩阵的每一行代表一个配体 - 受体对 # 矩阵的每一列代表不同细胞类型之间的相互作用组合。

列名的格式为 源细胞类型 - 目标细胞类型 head (mt @ data $ expr_l_r_log2_scale) # ST-ST-SSC ST-SPGing ST-SPGed SSC-ST SSC-SSC # DLL3-NOTCH1 0.1388065 0.38629302 0 0 0.3341891 0.4253563 # DLL3-NOTCH2 0.6705814 0.02475957 0 0 0.6929619 0.1910709 # DLL3-NOTCH3 0.5927026 0.00000000 0 0 0.6248139 0.0000000 # DLL3-NOTCH4 0.3381416 0.09253086 0 0 0.4340570 0.2236393 # NRG3-ERBB4 0.0000000 0.00000000 0 0 0.0000000 0.7625596 # ADCYAP1-ADCYAP1R1 0.0000000 0.58994436 0 0 0.0000000 0.0000000 # SSC-SPGing SSC-SPGed SPGing-ST SPGing-SSC SPGing-SPGing # DLL3-NOTCH1 0 0 0.3813847 0.4399023 0 # DLL3-NOTCH2 0 0 0.7017636 0.2342321 0 # DLL3-NOTCH3 0 0 0.6370433 0.0000000 0 # DLL3-NOTCH4 0 0 0.4641267 0.2609279 0 # NRG3-ERBB4 0 0 0.0000000 0.0000000 0 # ADCYAP1-ADCYAP1R1 0 0 0.0000000 0.0000000 0 # SPGing-SPGed-ST SPGed-SSC SPGed-SPGing SPGed-SPGed # DLL3-NOTCH1 0 0.3363955 0.4259846 0 0 # DLL3-NOTCH2 0 0.6933370 0.1930689 0 0 # DLL3-NOTCH3 0 0.6253393 0.0000000 0 0 # DLL3-NOTCH4 0 0.4354055 0.2253410 0 0 # NRG3-ERBB4 0 0.0000000 0.4667228 0 0 # ADCYAP1-ADCYAP1R1 0 0.0000000 0.0000000 0 0 dim (mt @ data $ expr_l_r_log2_scale) # [1] 264 16 # 由于互作对非常多，所以第一步我们先不展示行名 viewPheatmap ( # 输入绘图数据 object = mt, slot= "expr_l_r_log2_scale" , # 不显示行名 show_rownames = F, # 显示列名 show_colnames = T, # 设置热图行方向聚类树的高度。

将其设为0意味着在行方向上不显示聚类树 treeheight_row= 0 , # 设置热图列方向聚类树的高度。

这里设置为0意味着在列方向上不显示聚类树 treeheight_col= 0 , # 行聚类 cluster_rows = T, # 列不聚类 cluster_cols = F, # 字体大小 fontsize = 12 , # 字体角度 angle_col = "90" , # 图片标题 main= "score" ) # 筛选前十个互作对单独展示热图 selected_matrix <- mt @ data $ expr_l_r_log2_scale[ 1 : 10 , ] pheatmap (selected_matrix, # 行名称和列名称均展示出来 show_rownames = T, show_colnames = T, # 聚类树的高度设置为0 treeheight_row= 0 , treeheight_col= 0 , # 行列均不聚类 cluster_rows = F, cluster_cols = F, # 设置字体大小、展示角度和图片主标题 fontsize = 12 , angle_col = "90" , main= "score" , border_color = "white" ) 6.3 桑基图 # 进行配体 - 受体（LR）到转录因子（TF）的分析，探索细胞间通讯中配体 - 受体相互作用如何影响转录因子的调控 mt <- LR2TF ( object = mt, # 指定发送信号的细胞类型为 "ST" sender_cell= "ST" , # 指定接收信号的细胞类型为 "SSC" recevier_cell= "SSC" , slot= "expr_l_r_log2_scale" , # 分析物种 org= "Homo sapiens" ) # 查看文件前几行 head (mt @ reductions $ sankey) # Ligand Receptor TF weight1 weight2 # 2 BTC ERBB4 STAT5B 0.9076341 1.756423 # 3 BTC EGFR STAT5B 0.9020591 1.756423 # 4 INHBA ACVR1C SMAD3 0.8815241 1.539244 # 5 INHBA ACVR2B SMAD3 0.8576777 1.539244 # 6 INHBA ACVR1B SMAD3 0.8506733 1.539244 # 7 INHBA ACVR2A SMAD3 0.8448782 1.539244 # 加载R包 if ( ! require (networkD3)){ BiocManager :: install ( "networkD3" ) } # 用于创建桑基图对象，展示细胞间通讯的配体 - 受体 - 转录因子关系 sank <- LRT.Dimplot (mt, # 设置桑基图中文字的字体大小为 8 fontSize = 8 , # 设置桑基图中节点的宽度为 30 nodeWidth = 30 , # 设置桑基图的宽度为 1200 像素 width = 1200 , # 控制桑基图的流向，设置为 FALSE 时，桑基图的流向可能不是从左到右 sinksRight= FALSE , # 不使用自定义颜色，即使用函数默认的颜色方案 DIY.color = FALSE ) # 保存绘图结果 networkD3 :: saveNetwork (sank, "result_cellcall/ST-SSC_full.html" ) # 加载 magrittr 和 dplyr 包，dplyr 包用于数据处理和筛选，magrittr 包提供了管道操作符 %>% library (magrittr) library (dplyr) # 将 mt 对象中 reductions 槽位下的 sankey 数据赋值给 tmp <- mt @ reductions $ sankey # 使用 dplyr 包的 filter 函数筛选出 weight1 大于 0 的数据行，存储在 tmp1 中 tmp1 <- dplyr :: filter (tmp, weight1 > 0 ) # 调用 trans2tripleScore 函数将筛选后的数据 tmp1 转换为适合绘制桑基图的三元组数据格式，存储在 tmp.df 中 tmp.df <- trans2tripleScore (tmp1) head (tmp.df) # Ligand Receptor TF value # 1 sender:BTC receiver:ERBB4 TF:STAT5B 0.9076341 # 2 sender:BTC receiver:EGFR TF:STAT5B 0.9020591 # 3 sender:INHBA receiver:ACVR1C TF:SMAD3 0.8815241 # 4 sender:INHBA receiver:ACVR2B TF:SMAD3 0.8576777 # 5 sender:INHBA receiver:ACVR1B TF:SMAD3 0.8506733 # 6 sender:INHBA receiver:ACVR2A TF:SMAD3 0.8448782 # 设置桑基图的颜色 mycol.vector = c ( ' , ' , ' , ' , ' , ' , ' , ' , ' ) # 计算 tmp.df 中所有唯一元素的数量，即需要设置颜色的元素个数 elments.num <- tmp.df %>% unlist %>% unique %>% length # 将 mycol.vector 重复足够的次数，以确保有足够的颜色来为所有元素设置颜色 # ceiling(elements.num / length(mycol.vector))：计算最少需要重复多少次 mycol.vector 才能覆盖 elements.num 个元素（向上取整） mycol.vector.list <- rep (mycol.vector, times= ceiling (elments.num / length (mycol.vector))) # 由于整个数据框比较大，所以我们这里选择数据前30行绘制桑基图 sankey_graph ( df = tmp.df[ 1 : 30 ,], # 如果 sankey 的三次实现，设1:3, 否则1:2。

默认值为 1:3 axes= 1 : 3 , # 为桑基图的节点和边设置颜色，使用之前生成的颜色列表的前 elments.num 个颜色 mycol = mycol.vector.list[ 1 : elments.num], # 节点在 x 轴方向的偏移量，设置为 NULL 时使用默认偏移量 nudge_x = NULL , # 设置桑基图中文字的字体大小为 4 font.size = 4 , # 设置桑基图节点的边框颜色为白色 boder.col= "white" , # 如果为FALSE,第三个轴继承前两个轴，并且只考虑关系而不是分数。

否则，每个轴只继承第一个轴ggtitle, 并且只考虑分数 isGrandSon = F) 6.4 TF富集分析图 # 提取与细胞类型 “SSC” 相关的基因集（gene sets）名称，并查看前 3 个基因集的名称 ssc.tf <- names (mt @ data $ gsea.list $ SSC @ geneSets) ssc.tf[ 1 : 3 ] # [1] "FOXO1" "FOXO3" "GLI2" getGSEAplot ( # 问 mt 对象中 data 槽位下的 gsea.list 元素，该元素是一个列表，包含了不同细胞类型的 GSEA 结果 gsea.list= mt @ data $ gsea.list, # 指定要绘制可视化图的基因集 ID geneSetID= ssc.tf[ 1 : 3 ], # 指定要分析的细胞类型为 “SSC” myCelltype= "SSC" , # 传入包含基因表达倍数变化（fold change）信息的列表，用于在GSEA图中标记基因的表达变化情况 fc.list= mt @ data $ fc.list, # 指定要在图中突出显示的基因 ID，这里选择了 “SSC” 细胞类型的第一个基因集中的前 10 个基因 selectedGeneID = mt @ data $ gsea.list $ SSC @ geneSets[[ 1 ]][ 1 : 10 ], # 使用默认的配色方案 mycol = NULL ) 6.5 山脊图 # SSC细胞类型gsea富集分析结果 egmt <- mt @ data $ gsea.list $ SSC # 将egmt转换为数据框格式 egmt.df <- data.frame (egmt) head (egmt.df[, 1 : 6 ]) # ID Description setSize enrichmentScore NES pvalue # PPARG 481 0.4751570 1.635721 4.166312e-09 # STAT5B 257 0.5320769 1.778999 1.312772e-08 # TCF7 127 0.5990502 1.882576 1.732135e-07 # FOXO3 309 0.4618091 1.562171 1.452494e-05 # SMAD3 220 0.4658016 1.545993 1.248319e-04 # SMAD9 12 0.8304408 1.725800 4.905222e-04 # 保留p.adjust < 0.05的通路,并获取通路所在的位置索引 flag.index <- which (egmt.df $ p.adjust < 0.05 ) # 绘制山脊图 ridgeplot.DIY ( x= egmt, # 按照p.adjust填充 fill= "p.adjust" , # 展示满足条件的通路 showCategory= flag.index, # 使用核心富集的基因，即对富集结果贡献最大的基因 core_enrichment = T, # 按照NES值排序 orderBy = "NES" , # 按照升序排序 decreasing = FALSE ) 7、结果保存 Sys.time # [1] "2025-05-15 12:13:59 CST" save.image ( "result_cellcall/cellcall.rdata" )

七、icellnet 1、软件安装和加载 # install.packages("jetset") # install_github("soumelis-lab/ICELLNET",ref="master", subdir="icellnet") library (devtools) library ( "hgu133plus2.db" ) library (jetset) library (icellnet) library (Seurat) 2、数据准备 # 读取数据库 # 这一步骤可能在线读取不成功，因此我已经下载至本地，可以直接从本地读取数据库 # db=as.data.frame(read.csv(curl::curl(url="https://raw.githubusercontent.com/soumelis-lab/ICELLNET/master/data/ICELLNETdb.tsv"), sep="\t",header = T, check.names=FALSE, stringsAsFactors = FALSE, na.strings = "")) Sys.time # [1] "2025-05-15 12:14:27 CST" rm ( list = ls ) gc # used (Mb) gc trigger (Mb) max used (Mb) # Ncells 9408180 502.5 16491765 880.8 16491765 880.8 # Vcells 40256363 307.2 433370624 3306.4 845333542 6449.4 setwd ( "/data/07_scRNA-seq_cell_communication/" ) db <- read.table ( "data/ICELLNETdb.tsv" , sep= " \t " , # 文件分隔符为制表符 header = T, # 第一行作为列名称 check.names= FALSE , # 直接使用文件中的原始列名，不会对其进行修改 stringsAsFactors = FALSE , # 字符型数据会以字符向量的形式存储在数据框中，而不会被转换为因子 na.strings = "" ) # 将空字符串""定义为缺失值 # name.lr.couple函数从配体 / 受体相互作用的数据库中检索相应的基因符号。

它返回一个表，其中包含配对 L/R 的名称以及相互作用所属分子的家族 / 亚家族。

db.name.couple <- name.lr.couple (db, type= "Family" ) head (db.name.couple) # Pair Family # [1,] "TIMP1 + MMP9 / LRP1" NA # [2,] "LTB + LTA / LTBR" "Cytokine" # [3,] "ITGAV + ITGB5 / ADGRB1" "Cell adhesion" # [4,] "ITGAV + ITGB3 / ADGRA2" "Cell adhesion" # [5,] "ITGAV + ITGB3 / ADGRE5" "Cell adhesion" # [6,] "ITGA5 + ITGB1 / ADGRE5" "Cell adhesion" # 读取单细胞数据 # 数据https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE222703 # 人类透明细胞肾细胞癌中血管原依赖性肿瘤 # 数据同样已经保存至本地，可以直接从本地读取rds文件 seurat <- readRDS ( "data/GSE222703_SeuratObj_ALL_integrated.rds" ) seurat # An object of class Seurat # 28140 features across 27963 samples within 1 assay # Active assay: RNA (28140 features, 2000 variable features) # 3 layers present: counts, data, scale.data # 4 dimensional reductions calculated: pca, umap, harmony, harmony_umap seurat <- NormalizeData (seurat) seurat <- ScaleData (seurat) seurat # An object of class Seurat # 28140 features across 27963 samples within 1 assay # Active assay: RNA (28140 features, 2000 variable features) # 3 layers present: counts, data, scale.data # 4 dimensional reductions calculated: pca, umap, harmony, harmony_umap # 展示降维图 DimPlot (seurat, reduction = 'umap' , group.by = 'Celltype_Harmony2' , split.by= "Tissue" , label = T) 3、数据预处理 # 选择来源于肿瘤样本的数据 Idents (seurat) <- seurat $ Celltype_Harmony2 seurat.tum <- subset (seurat, cells= which (seurat $ Tissue == "Tumor" )) seurat.tum # An object of class Seurat # 28140 features across 15113 samples within 1 assay # Active assay: RNA (28140 features, 2000 variable features) # 3 layers present: counts, data, scale.data # 4 dimensional reductions calculated: pca, umap, harmony, harmony_umap DimPlot (seurat.tum, reduction = 'umap' , label = T) # 过滤 scRNAseq 数据表达比例小于0.1的基因 # 计算每个基因在每个细胞类型的表达比例和表达均值 # 删除 ICellNet 数据库中不存在的所有基因 # 约30min dir.create ( "result_icellnet" ) data <- sc.data.cleaning ( object = seurat.tum, # 用于分析的配体-受体数据库 db = db, # 基因必须在至少 10% 的细胞中被检测到，才会保留 filter.perc = 10 , # 保存结果文件 save_file = T, # 强制覆盖已有的文件 force.file = T, path= "/data/07_scRNA-seq_cell_communication/result_icellnet" ) # [1] "Filling intermediate table: percentage of expressing cell per cluster per gene, and mean of expression" # [1] "Intermediate table were saved as scRNAseq_statsInfo_for_ICELLNET.csv." # [1] "Filtering done" # 对输入的基因表达数据 data 进行缩放（scaling），将表达值转换为统一尺度（0~1 之间） files <- list.files ( "/data/07_scRNA-seq_cell_communication" ) # 筛选出文件名包含statsInfo的文件 icellnet_files <- files[ grepl ( "statsInfo" , files)] # 移动文件 file.rename ( from = file.path ( "/data/07_scRNA-seq_cell_communication" , icellnet_files), to = file.path ( "/data/07_scRNA-seq_cell_communication/result_icellnet" , icellnet_files)) # [1] TRUE data.icell <- as.data.frame ( gene.scaling (data, n= 1 , db= db)) 4、在感兴趣的集群上应用icellnet管道 研究来自ccRCC1和ccRCC2簇的肿瘤细胞与TME中其他细胞的通讯。

这意味着我们要考虑肿瘤细胞表达的配体和TME中其他细胞表达的受体。

# CC：配体细胞，这里指定两种 CC1 <- "ccRCC1" CC2 <- "ccRCC2" # PC：受体细胞类型 PC <- c ( "B_cell" , "PlasmaC" , "Treg" , "CD4T" , "Prolif" , "CD8T" , "NK_T" , "NK" , "NK_CD160" , "Mono_nc" , "Mono_c" , "Mac_C1QA" , "cDC2" , "cDC1" , "pDC" , "Neutrop" , "Mast" , "PT_GPX3" , "PT_MT1G" , "ccRCC1" , "ccRCC2" , "Epith" , "Endoth" , "Fibro" ) # 信号得分方向 # "out" 表示从发送细胞（如 CC1、CC2）向 PC 中细胞释放配体 direction <- "out" # 数据筛选 # 从标准化表达数据 data.icell 中提取两列数据，CC1的表达值列和Symbol列 CC1.data <- as.data.frame (data.icell[, c (CC1, "Symbol" )], row.names = rownames (data.icell)) # 从标准化表达数据 data.icell 中提取两列数据，一列是 CC2的表达值列和Symbol列 CC2.data <- as.data.frame (data.icell[, c (CC2, "Symbol" )], row.names = rownames (data.icell)) # 从标准化表达数据 data.icell 中提取两列数据，所有PC中包含的细胞类型的表达值列和Symbol列 PC.data <- as.data.frame (data.icell[, c (PC, "Symbol" )], row.names = rownames (data.icell)) 5、计算细胞间通讯分数 # 计算第一组细胞通讯数据的得分 score.computation .1 <- icellnet.score ( # 分析 ccRCC1 向外释放信号（即它是配体发送方） direction= direction, # 所有接收细胞的表达数据 PC.data= PC.data, # ccRCC1 的表达数据 CC.data= CC1.data, # 接收细胞名称向量 PC= PC, # 配体-受体数据库 db = db) # 返回结果是一个列表，包含两个部分 # [[1]]：细胞通信得分矩阵 # [[2]]：对应的配体-受体对信息 score1 <- as.data.frame (score.computation .1 [[ 1 ]]) lr1 <- score.computation .1 [[ 2 ]] # 计算第二组细胞通讯数据的得分 score.computation .2 <- icellnet.score ( # 分析 ccRCC2 向外释放信号（即它是配体发送方） direction= direction, # 所有接收细胞的表达数据 PC.data= PC.data, # ccRCC2 的表达数据 CC.data= CC2.data, # 接收细胞名称向量 PC= PC, # 配体-受体数据库 db = db) # 返回结果是一个列表，包含两个部分 # [[1]]：细胞通信得分矩阵 # [[2]]：对应的配体-受体对信息 score2 <- as.data.frame (score.computation .2 [[ 1 ]]) lr2 = score.computation .2 [[ 2 ]] # 将两个得分矩阵（score1 和 score2）按列合并 Scores <- cbind (score1, score2) # 结果 Scores 是一个数据框，行是受体细胞类型，列是每个发送细胞的打分 # 可以直观看出哪一个发送细胞对不同受体细胞具有更高的通信能力 colnames (Scores) <- c ( paste0 (CC1, "_" , direction), paste0 (CC2, "_" , direction)) 6、结果可视化 6.1 热图 # 调整y轴，计算最大得分 +1，用于图形美观或坐标设定 ymax = max (Scores) + 1 # 对配体-受体对按“家族/类型”归类统计打分并求和 # 分析哪些配体-受体“家族”（例Cytokine等）在细胞通信中占据主导 LR.family.score ( lr= lr1, db.couple= db.name.couple, plot= NULL ) # B_cell PlasmaC Treg CD4T Prolif CD8T # Cytokine 165.905854 103.05001 117.177639 112.889115 144.676093 153.03583 # Cell adhesion 14.085539 26.77468 36.509652 25.915522 35.376775 27.37436 # Innate immune 0.000000 0.00000 0.000000 0.000000 0.000000 0.00000 # ECM interaction 0.000000 0.00000 0.000000 0.000000 0.000000 0.00000 # Growth factor 0.000000 0.00000 0.000000 0.000000 6.834163 0.00000 # Wnt pathway 0.000000 0.00000 0.000000 0.000000 0.000000 0.00000 # HLA recognition 30.253765 27.74285 6.464834 9.308346 6.217381 0.00000 # Chemokine 55.700981 85.07079.631575 72.847661 58.804015 67.74178 # Checkpoint 6.996102 73.07775 36.365148 7.771715 210.412861 261.97994 # Notch pathway 0.000000 0.00000 0.000000 0.000000 0.000000 0.00000 # other 125.449652 70.07518 69.448787 69.824719 151.554209 97.71089 # NK_T NK_CD160 Mono_nc Mono_c Mac_C1QA # Cytokine 126.09611 99.37361 118.72436 148.455641 170.436521 164.38255 # Cell adhesion 22.82098 25.88223 22.86799 21.470079 30.796470 35.79270 # Innate immune 0.00000 0.00000 0.00000 100.854638 119.235387 61.70336 # ECM interaction 0.00000 0.00000 0.00000 0.000000 0.000000 0.00000 # Growth factor 0.00000 0.00000 0.00000 0.000000 0.000000 63.05492 # Wnt pathway 0.00000 0.00000 0.00000 0.000000 0.000000 0.00000 # HLA recognition 19.40808 69.38913 64.69043 184.959258 107.068337 72.60399 # Chemokine 56.62170 38.73531 48.09783 38.677382 99.356130 42.19831 # Checkpoint 97.49204 300.66331 47.43266 3.856249 6.590529 11.34261 # Notch pathway 0.00000 0.00000 0.00000 0.000000 0.000000 0.00000 # other 79.38896 85.89410 72.53888 268.238086 243.004426 151.14624 # cDC2 cDC1 pDC Neutrop Mast PT_GPX3 # Cytokine 188.10688 203.36807 183.04583 262.67747 40.296735 10.017060 # Cell adhesion 21.76423 21.37588 10.60357 0.00000 41.202943 18.215311 # Innate immune 22.94667 11.79295 0.00000 193.78318 0.000000 0.000000 # ECM interaction 0.00000 0.00000 0.00000 0.00000 0.000000 0.000000 # Growth factor 0.00000 0.00000 20.16138 0.00000 16.218563 0.000000 # Wnt pathway 0.00000 0.00000 0.00000 0.00000 0.000000 0.000000 # HLA recognition 48.50512 15.53833 44.54312 47.31583 7.061686 0.000000 # Chemokine 39.69463 31.90878 50.08308 58.57526 10.754984 0.000000 # Checkpoint 13.31953 20.10259 6.26097 0.00000 12.212570 0.000000 # Notch pathway 0.00000 0.00000 0.00000 0.00000 0.000000 0.000000 # other 150.49162 135.33151 107.83579 207.75686 127.417397 8.160692 # PT_MT1G ccRCC1 ccRCC2 Epith Endoth Fibro # Cytokine 34.61355 33.469314 132.248973 15.897350 133.66357 55.688991 # Cell adhesion 0.00000 0.000000 108.267066 71.447661 192.75620 164.829587 # Innate immune 0.00000 0.000000 0.000000 0.000000 14.10676 0.000000 # ECM interaction 0.00000 0.000000 0.000000 0.000000 0.00000 0.000000 # Growth factor 0.00000 34.808569 158.012798 0.000000 697.47607 179.171817 # Wnt pathway 0.00000 0.000000 0.000000 0.000000 0.00000 0.000000 # HLA recognition 0.00000 0.000000 0.000000 0.000000 46.43686 5.864988 # Chemokine 0.00000 5.129455 7.512258 4.282749 16.13692 0.000000 # Checkpoint 0.00000 0.000000 6.044988 0.000000 0.00000 0.000000 # Notch pathway 0.00000 0.000000 0.000000 0.000000 0.00000 0.000000 # other 0.00000 5.000082 125.313711 121.596316 226.48782 116.976253 # 生成配体-受体家族热图 # y轴：配体-受体家族; x轴：不同细胞类型 # 热图颜色强度代表该家族的通信活跃度，可以直观看到哪些家族主导了ccRCC1 → PC的通信 # 柱形图表示每个细胞类型的得分总和 LR.family.score ( lr= lr1, db.couple= db.name.couple, plot= "heatmap" ) 6.2 柱形图 # 使用icellnet中的LR.family.score函数绘制配体-受体家族的柱状图（barplot），并为每个家族指定不同的颜色 # 使用 RColorBrewer 的调色板 "Set2"（共8种柔和颜色） # colorRampPalette(...) 用来扩展这些颜色，使颜色数量足够多 # 输出的 colors 是一个颜色向量，长度与家族数目一致 colors = grDevices :: colorRampPalette (RColorBrewer :: brewer.pal ( 8 , "Set2" ))( length ( unique (db $ Family))) LR.family.score ( lr= lr1, db.couple= db.name.couple, plot= "barplot" , # 绘制柱形图 family.col= colors) # 使用指定颜色 7、在TME内识别ccRCC2特定的通信信道 7.1 计算ccRCC2的特定相互作用 # 构建 global_LR 矩阵（配体-受体对在不同细胞中的强度矩阵） global_LR = matrix ( nrow= length ( rownames (db)), ncol= length (PC)) # 行：所有的配体-受体对（来自 db.name.couple[,1]） rownames (global_LR) = db.name.couple[, 1 ] # 列：所有的接收细胞类型（PC） colnames (global_LR) = colnames ( # LR.viz(..., int = i)：对每一个配体-受体对 i，调用 LR.viz 计算它在各个细胞间的通信矩阵（如 sender → receiver），得到的矩阵行为发送细胞（sender cell），列为接收细胞（receiver cell） LR.viz ( # 从 data.icell 中筛选出所有的目标细胞（即 PC 向量中列出的细胞类型）和 Symbol（即基因名）列，构建出一个只包含感兴趣细胞和基因的表达矩阵 data = data.icell[, colnames (data.icell) %in% c (PC, "Symbol" )], # 选择一个特定的配体-受体对（这里是第一对 L-R 对），并将其转换为字符型数据格式 int= as.character (db.name.couple[ 1 , 1 ]), db= db, plot= F)) direction = "out" # 对每一个配体-受体对 i，调用 LR.viz 计算它在各个细胞间的通信矩阵（如 sender → receiver） for (i in db.name.couple[, 1 ]){ if (direction == "out" ){ # 使用 rowSums（方向 out）或 colSums（方向 in）来提取每个细胞在该 L-R 对中发出或接收信号的总量 # 将结果存入 global_LR[i, ]，即这个配体-受体对在所有目标细胞中的强度 global_LR[i,] = rowSums ( LR.viz ( data = data.icell[, colnames (data.icell) %in% c (PC, "Symbol" )], int= i, db= db, plot= F)) } if (direction == "in" ){ global_LR[i,] = colSums ( LR.viz ( data = data.icell[, colnames (data.icell) %in% c (PC, "Symbol" )], int= i, db= db, plot= F)) } } # 评估相互作用的细胞特异性 int.spe <- function (LR.mat, cell, thresh){ # 找到每一对 L-R 的最大响应细胞，通过计算在哪个细胞上通信得分最高 colMax = apply (LR.mat, 1 , which.max) # 只保留在目标细胞 cell 上作用最强的那些 L-R 对 a = which ( colnames (LR.mat) == cell) LR.mat2 = LR.mat[ which (colMax == a),] “cell” 的最大值 value = data.frame ( "vmax" = LR.mat2[,a]) rownames (value) = rownames (LR.mat2) # 计算它们与“次强细胞”的比值 # 次强细胞得分 +0.1（防除0） value $ vmax2 = apply (LR.mat2[, - a], 1 , function (x) ( max (x) + 0.1 )) value $ vmax2_nb = apply (LR.mat2[, - a], 1 , which.max) value $ vmax2_cell = colnames (LR.mat2[, - a])[value $ vmax2_nb] value = value %>% dplyr :: select ( - c ( "vmax2_nb" )) value $ vmean_pos = apply ( LR.mat2[, - a], 1 , function (x) mean (x[x != 0 ])) # 得分比值（越大表示越特异） value $ ratio = value $ vmax / value $ vmax2 # 过滤掉不够特异的配体-受体对 value = value %>% filter (ratio > thresh) # 最终返回的是对目标细胞最特异性的 L-R 对及其相关信息 return (value) } thresh = 1.5 # 给出一个相互作用的列表，其中 ccRCC2 的贡献是其他细胞类型的 1.5 倍 specific_ccRCC2 <- int.spe (global_LR, cell= "ccRCC2" , thresh= thresh) 7.2 可视化 7.2.1 热图展示所有特异互作对 热图展示ccRCC2专门用于与TME其他细胞类型相互作用的所有相互作用 LR.heatmap ( lr = lr2[ rownames ( na.omit (specific_ccRCC2)),], thresh = 0 , # 控制热图上显示的配体-受体对的数量 topn= 20 , # "var" 表示按照配体-受体对的通信得分的方差进行排序 sort.by= "var" , # 图片标题 title= "ccRCC2 specific interactions" ) 7.2.2 热图展示指定互作对 spe1 <- LR.viz ( data= data.icell, db = db, int= "HHLA2 / TMIGD2" , plot= T) spe1 spe2 <- LR.viz ( data= data.icell, db = db, int= "CD70 / CD27" , plot= T) spe2 spe3 <- LR.viz ( data= data.icell, db = db, int= "VEGFA / NRP1" , plot= T) spe3 8、结果保存 Sys.time # [1] "2025-05-15 12:45:10 CST" save.image ( "result_icellnet/icellnet.rdata" ) 环境版本信息 sessionInfo # R version 4.4.2 (2024-10-31) # Platform: x86_64-pc-linux-gnu # Running under: Ubuntu 20.04.6 LTS # # Matrix products: default # BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 # LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0 # # locale: # [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C # [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 # [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 # [7] LC_PAPER=en_US.UTF-8 LC_NAME=C # [9] LC_ADDRESS=C LC_TELEPHONE=C # [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C # # time zone: Asia/Shanghai # tzcode source: system (glibc) # # attached base packages: # [1] stats4 grid stats graphics grDevices utils datasets # [8] methods base # # other attached packages: # [1] icellnet_2.2.1 jetset_3.4.0 hgu133plus2.db_3.13.0 # [4] org.Hs.eg.db_3.16.0 AnnotationDbi_1.68.0 IRanges_2.32.0 # [7] S4Vectors_0.36.0 devtools_2.4.5 usethis_2.1.6 # [10] magrittr_2.0.3 networkD3_0.4 pheatmap_1.0.12 # [13] cellcall_1.0.7 RColorBrewer_1.1-3 forcats_1.0.0 # [16] stringr_1.5.1 purrr_1.0.2 readr_2.1.5 # [19] tidyr_1.3.1 tibble_3.2.1 tidyverse_1.3.2 # [22] Seurat_5.0.3 SeuratObject_5.0.1 sp_2.1-4 # [25] nichenetr_2.2.0 ComplexHeatmap_2.14.0 cowplot_1.1.3 # [28] ggalluvial_0.12.5 patchwork_1.3.0 CellChat_1.6.1 # [31] bigmemory_4.6.4 Biobase_2.66.0 BiocGenerics_0.52.0 # [34] ggplot2_3.5.1 igraph_2.0.1.1 dplyr_1.1.4 # # loaded via a namespace (and not attached): # [1] ica_1.0-3 plotly_4.10.1 Formula_1.2-5 # [4] zlibbioc_1.52.0 tidyselect_1.2.1 rvest_1.0.3 # [7] bit_4.5.0.1 doParallel_1.0.17 clue_0.3-66 # [10] lattice_0.22-6 rjson_0.2.23 urlchecker_1.0.1 # [13] blob_1.2.4 rngtools_1.5.2 parallel_4.4.2 # [16] caret_6.0-93 png_0.1-8 cli_3.6.3 # [19] ggplotify_0.1.0 usdata_0.3.1 registry_0.5-1 # [22] goftest_1.2-3 gargle_1.2.1 textshaping_0.3.6 # [25] BiocNeighbors_2.0.1 ggnetwork_0.5.13 uwot_0.1.14 # [28] shadowtext_0.1.4 mime_0.12 evaluate_1.0.3 # [31] tidytree_0.4.1 leiden_0.4.3 openintro_2.4.0 # [34] stringi_1.8.4 pROC_1.18.0 backports_1.5.0 # [37] lubridate_1.9.0 httpuv_1.6.6 clusterProfiler_4.6.0 # [40] splines_4.4.2 prodlim_2019.11.13 jpeg_0.1-10 # [43] ggraph_2.1.0 DT_0.26 sctransform_0.4.1 # [46] ggbeeswarm_0.7.2 sessioninfo_1.2.2 DBI_1.2.3 # [49] jquerylib_0.1.4 withr_3.0.2 class_7.3-20 # [52] systemfonts_1.1.0 enrichplot_1.18.1 lmtest_0.9-40 # [55] ggnewscale_0.5.0 tidygraph_1.2.2 BiocManager_1.30.25 # [58] htmlwidgets_1.5.4 fs_1.6.5 ggrepel_0.9.6 # [61] cherryblossom_0.1.0 statnet.common_4.11.0 labeling_0.4.3 # [64] cellranger_1.1.0 reticulate_1.26 zoo_1.8-12 # [67] XVector_0.38.0 knitr_1.49 network_1.18.0 # [70] airports_0.1.0 timechange_0.3.0 foreach_1.5.2 # [73] caTools_1.18.3 visNetwork_2.1.2 data.table_1.16.4 # [76] timeDate_4041.110 ggtree_3.6.2 psych_2.4.12 # [79] RSpectra_0.16-2 irlba_2.3.5.1 DiagrammeR_1.0.11 # [82] ggrastr_1.0.1 fastDummies_1.7.4 gridGraphics_0.5-1 # [85] ellipsis_0.3.2 lazyeval_0.2.2 yaml_2.3.10 # [88] survival_3.4-0 scattermore_1.2 crayon_1.5.3 # [91] RcppAnnoy_0.0.22 progressr_0.15.1 tweenr_2.0.3 # [94] later_1.4.1 profvis_0.3.7 ggridges_0.5.6 # [97] codetools_0.2-20 base64enc_0.1-3 GlobalOptions_0.1.2 # [100] KEGGREST_1.38.0 Rtsne_0.17 shape_1.4.6.1 # [103] limma_3.62.2 foreign_0.8-88 pkgconfig_2.0.3 # [106] xml2_1.3.6 ggpubr_0.5.0 aplot_0.1.8 # [109] spatstat.sparse_3.1-0 ape_5.8-1 viridisLite_0.4.2 # [112] gridBase_0.4-7 xtable_1.8-4 interp_1.1-6 # [115] car_3.1-1 plyr_1.8.9 httr_1.4.4 # [118] tools_4.4.2 globals_0.16.3 hardhat_1.4.0 # [121] pkgbuild_1.3.1 beeswarm_0.4.0 htmlTable_2.4.1 # [124] broom_1.0.7 checkmate_2.3.2 nlme_3.1-160 # [127] HDO.db_0.99.1 dbplyr_2.5.0 crosstalk_1.2.1 # [130] digest_0.6.37 Matrix_1.7-1 farver_2.1.2 # [133] tzdb_0.4.0 reshape2_1.4.4 ModelMetrics_1.2.2.2 # [136] yulab.utils_0.1.9 viridis_0.6.5 rpart_4.1.24 # [139] glue_1.8.0 cachem_1.1.0 polyclip_1.10-7 # [142] Hmisc_4.7-2 generics_0.1.3 Biostrings_2.66.0 # [145] googledrive_2.0.0 presto_1.0.0 parallelly_1.41.0 # [148] pkgload_1.4.0 mnormt_2.1.1 statmod_1.5.0 # [151] RcppHNSW_0.6.0 ragg_1.4.0 carData_3.0-5 # [154] pbapply_1.7-2 spam_2.11-0 gson_0.1.0 # [157] utf8_1.2.4 gower_1.0.2 graphlayouts_0.8.3 # [160] readxl_1.4.3 ggsignif_0.6.4 gridExtra_2.3 # [163] shiny_1.10.0 lava_1.7.0 GenomeInfoDbData_1.2.13 # [166] RCurl_1.98-1.16 memoise_2.0.1 rmarkdown_2.29 # [169] downloader_0.4 scales_1.3.0 googlesheets4_1.0.1 # [172] future_1.34.0 svglite_2.1.0 RANN_2.6.2 # [175] Cairo_1.6-2 bigmemory.sri_0.1.8 spatstat.data_3.1-4 # [178] rstudioapi_0.17.1 cluster_2.1.8 spatstat.utils_3.1-2 # [181] hms_1.1.3 fitdistrplus_1.2-2 munsell_0.5.1 # [184] fdrtool_1.2.18 colorspace_2.1-1 FNN_1.1.4.1 # [187] rlang_1.1.4 GenomeInfoDb_1.34.3 ipred_0.9-13 # [190] dotCall64_1.2 ggforce_0.4.2 circlize_0.4.16 # [193] xfun_0.50 coda_0.19-4.1 e1071_1.7-16 # [196] sna_2.7-2 modelr_0.1.10 remotes_2.5.0 # [199] recipes_1.0.3 iterators_1.0.14 matrixStats_1.5.0 # [202] abind_1.4-8 randomForest_4.7-1.2 GOSemSim_2.24.0 # [205] treeio_1.22.0 bitops_1.0-9 ps_1.8.1 # [208] promises_1.3.2 scatterpie_0.1.8 RSQLite_2.3.9 # [211] qvalue_2.38.0 reprex_2.0.2 fgsea_1.24.0 # [214] proxy_0.4-27 GO.db_3.16.0 compiler_4.4.2 # [217] prettyunits_1.2.0 listenv_0.9.1 Rcpp_1.0.14 # [220] tensor_1.5 MASS_7.3-64 uuid_1.2-1 # [223] BiocParallel_1.40.0 spatstat.random_3.0-1 R6_2.5.1 # [226] fastmap_1.2.0 fastmatch_1.1-6 rstatix_0.7.1 # [229] vipor_0.4.7 ROCR_1.0-11 nnet_7.3-20 # [232] gtable_0.3.6 KernSmooth_2.23-26 latticeExtra_0.6-30 # [235] miniUI_0.1.1.1 deldir_2.0-4 htmltools_0.5.8.1 # [238] bit64_4.5.2 spatstat.explore_3.0-5 lifecycle_1.0.4 # [241] processx_3.8.5 callr_3.7.6 sass_0.4.9 # [244] vctrs_0.6.5 spatstat.geom_3.0-3 DOSE_3.24.1 # [247] haven_2.5.4 NMF_0.28 ggfun_0.0.8 # [250] future.apply_1.11.3 bslib_0.9.0 pillar_1.10.1 # [253] magick_2.8.5 jsonlite_1.8.9 GetoptLong_1.0.5 参考：https://www.nature.com/articles/s41392-024-01888-z https://www.nature.com/articles/s41467-021-21246-9.pdf https://github.com/ShellyCoder/cellcall https://github.com/soumelis-lab/ICELLNET?tab=readme-ov-file