第 5/8 章

六、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" )

← 上一章 下一章 →