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