CellChat通过综合信号配体、受体及其辅因子基因的表达只与它们之间互作的先验知识对细胞通讯概率建模。在推断出细胞间通讯网络后，CellChat提供了进一步数据探索、分析和可视化功能。

CellChat工作流程图：

一、单样本分析(可以是同一组的多个生物学重复一起分析)

1.数据准备
library(Seurat)
library(SeuratData)
library(tidyverse)
library(CellChat)
library(NMF)
library(ggalluvial)
library(patchwork)
library(ggplot2)
library(svglite)
options(stringsAsFactors = F)

pbmc3k.final <- readRDS("pbmc.rds")
library(mindr)
(out <- capture.output(str(pbmc3k.final)))
out2 <- paste(out, collapse = "\n")
mm(gsub("\\.\\.@", "# ", gsub("\\.\\. ", "#", out2)), type = "text", root = "Seurat")
2.创建一个CellChat对象

从Seurat对象直接创建：构建CellChat对象时，输入的是log后的数据

cellchat <- createCellChat(pbmc3k.final@assays$RNA@data, meta = pbmc3k.final@meta.data, group.by = "cell_type")
summary(cellchat)
levels(cellchat@idents)
groupSize <- as.numeric(table(cellchat@idents))
groupSize
3.导入配体受体数据库
CellChatDB <- CellChatDB.human
str(CellChatDB)
# 包含interaction、complex、cofactor和geneInfo这四个dataframe
colnames(CellChatDB$interaction)
CellChatDB$interaction[1:4,1:4]
showDatabaseCategory(CellChatDB)

在CellChat中，我们还可以先选择特定的信息描述细胞间的相互作用，可以理解为从特定的侧面来刻画细胞间相互作用，比用一个大的配体库又精细了许多。

unique(CellChatDB$interaction$annotation)
CellChatDB.use <- subsetData(CellChatDB, search = "Secreted Signaling")
cellchat@DB <- CellChatDB.use
4.预处理库

对表达数据进行预处理，用于细胞间的通讯分析。首先在一个细胞组中识别过表达的配体或受体，然后将基因表达数据投射到蛋白-蛋白相互作用网络(PPI)上。如果配体或受体过表达，则识别过表达配体和受体之间的相互作用。

# 在矩阵的所有基因中提取signaling gene，结果保存在data.signaling(13714个基因，过滤完只有270个)
cellchat <- subsetData(cellchat)
# 相当于Seurat的FindMarkers，找每个细胞群中高表达的配体受体
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
# identify over-expressed ligand-receptor interactions (pairs) within the used CellChatDB
cellchat <- identifyOverExpressedInteractions(cellchat)
# 上一步运行的结果存储在cellchat@LR$LRsig
cellchat <- projectData(cellchat, PPI.human)
# 找到配体受体关系后，projectData将配体受体对的表达质投射到PPI上，来对@data.signaling中的表达值进行校正。
# 结果保存在@data.project
5.推断细胞通讯网络

通过为每个相互作用分配一个概率值并进行置换检验来推断生物意义上的细胞-细胞通讯。

通过计算与每个信号通路相关的所有配体-受体相互作用的通信概率来推断信号通路水平上的通信概率。

# 根据表达值推测细胞互作的概率(cellphonedb是用平均表达值代表互作强度)
# 如果不想用上一步PPI矫正的结果，raw.use = TRUE即可
cellchat <- computeCommunProb(cellchat, raw.use = F, population.size = T)
# Filter out the cell-cell communication if there are only few number of cells in certain cell groups
cellchat <- filterCommunication(cellchat, min.cells = 10)
df.net <- subsetCommunication(cellchat)
openxlsx::write.xlsx(df.net, "net_lr.xlsx")

- 推断信号通路水平细胞通讯网络 (结果储存存在net下面，有一个概率值和对应的pval)

我们可以通过计算链路的数量或汇总通信概率来计算细胞间的聚合通信网络。

cellchat <- computeCommunProbPathway(cellchat)
df.netp <- subsetCommunication(cellchat, slot.name = "netP")
openxlsx::write.xlsx(df.netp, file = "net_pathway.xlsx")

至此，统计分析部分已经完成。

6.细胞互作关系展示
6.1所有细胞群总体观：细胞互作数量与强度统计分析
# 统计细胞和细胞之间通信的数量(有多少个配体-受体对)和强度(概率)
cellchat <- aggregateNet(cellchat)
# 计算每种细胞各有多少个
groupSize <- as.numeric(table(cellchat@idents))
par(mfrow = c(1,2), 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")

左图：外周各种颜色圆圈的大小表示细胞的数量，圈越大，细胞数量越多。发出箭头的细胞表达配体，箭头指向的细胞表达受体。配体-受体对越多，线越粗。右图：互作的概率/强度值(强度就是概率值相加)

检查每种细胞发出的信号

mat <- cellchat@net$count
par(mfrow = c(3,3), 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, arrow.width = 0.2,
                   arrow.size = 0.1, edge.weight.max = max(mat), title.name = rownames(mat)[i])
}

mat <- cellchat@net$weight
par(mfrow = c(3,3), xpd = T)
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, arrow.width = 0.2,
                   arrow.size = 0.1, edge.weight.max = max(mat), title.name = rownames(mat)[i])
}

mat <- cellchat@net$weight
par(mfrow = c(2,3), xpd = T)
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, arrow.width = 0.2,
                   arrow.size = 0.1, edge.weight.max = max(mat), title.name = rownames(mat)[i])
}
6.2单个信号通路或配体-受体介导的细胞互作可视化(层次图、网络图、和弦图、热图)

比如在前面的功能富集分析或case control的比较中找到了一些信号通路差异，就可以进一步在细胞水平上验证。

cellchat@netP$pathways
pathways.show <- c("TGFb")

层次图(Hierarchy plot)

levels(cellchat@idents)
# define a numeric vector (淋巴系细胞) giving the index of the celltype as targets
vertex.receiver <- c(1,2,4,6)
netVisual_aggregate(cellchat, signaling = pathways.show, vertex.receiver = vertex.receiver)

在层次图中，实体圆和空心圆分别表示源和目标。圆的大小与每个细胞组的细胞数成比例。线越粗，互作信号越强。左图中间的target是我们选定的靶细胞。右图是选中的靶细胞之外的另外一组放在中间看互作。

网络图(circle plot)

par(mfrow = c(1,1))
netVisual_aggregate(cellchat, signaling = pathways.show, layout = 'circle)

和弦图(chord diagram)

par(mfrow = c(1,1))
netVisual_aggregate(cellchat, signaling = pathways.show, layout = 'chord)

热图(Heatmap)

par(mfrow = c(1,1))
netVisual_heatmap(cellchat, signaling = pathways.show, color.heatmap = 'Reds)

纵轴是发出信号的细胞，横轴是接收信号的细胞，热图颜色深浅代表信号强度。上侧和右侧的柱子是纵轴和横轴强度的累计

6.3配体-受体层级的可视化(计算各个ligand-receptor pair对信号通路的贡献)
# 计算配体受体对选定信号通路的贡献值(在这里就是查看哪条信号通路对TGFβ贡献最大)
netAnalysis_contribution(cellchat, signaling = pathways.show)
# 提取对TGFb有贡献的所有配体受体
pairLR.TGFb <- extractEnrichedLR(cellchat, signaling = pathways.show, geneLR.return = F)

层次图(Hierarchy plot)

# 提取对这个通路贡献最大的配体受体对来展示(也可以选择其他的配体受体对)
LR.show <- pairLR.TGFb[1,]
vertex.receiver <- c(1,2,4,6)
netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, vertex.receiver = vertex.receiver)

网络图(Circle plot)

netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = 'circle)

和弦图(chord diagram)

netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = 'chord)
6.4自动(批量)保存每个信号通路的互作结果
# Access all the signaling pathways showing signigficant communication
pathways.show.all <- cellchat@netP$pathways
# check the order of cell identity to set suitable vertex.receiver
levels(cellchat@idents)
vertex.receiver <- c(1,2,4,6)
dir.create("all_pathways_com_circle")
setwd("all_pathways_com_circle")
for(i in 1:length(pathways.show.all)){
  # Visualize communication network associated with both signaling pathway and individual L-R pairs
  netVisual(cellchat, signaling = pathways.show.all[i], out.format = c('pdf'),
            vertex.receiver = vertex.receiver, layout = 'circle')
  # compute and visualize the contribution of each ligand-receptor pair to the overall signaling pathway
  gg <- netAnalysis_contribution(cellchat, signaling = pathways.show.all[i])
  ggsave(filename = paste0(pathways.show.all[i],'_L-R_contribution.pdf'), plot = gg, width = 5, height = 2.5, units = 'in', dpi = 300)
}
setwd("../")
6.5多个配体-受体介导的细胞互作关系可视化

气泡图(全部配体受体)

levels(cellchat@idents)
# show all the significant interactions (LR pairs)
p <- netVisual_bubble(cellchat, sources.use = c(3,5,7,8,9), targets.use = c(1,2,4,6), remove.isolate = F)
ggsave("Mye_Lymph_bubble.pdf",p,width = 8, height = 12)  # 髓系对淋巴的调节

气泡图(指定信号通路或配体-受体)

# 比如制定CCL和CXCL这两个信号通路
netVisual_bubble(cellchat, sources.use = c(3,5,7,8,9),targets.use = c(1,2,4,6),signaling = c("CCL","CXCL"),remove.isolate = F)

气泡图(指定信号通路或配体-受体并制定细胞)

# show all the significant interactions (LR pairs) based on user's input
pairLR.use <- extractEnrichedLR(cellchat, signaling = c("CCL","CXCL","TGFb"))
netVisual_bubble(cellchat, sources.use = c(3,6,8),targets.use = c(1,4,5),pairLR.use = pairLR.use, remove.isolate = T)

参与某条信号通路(如TGFb)的所有基因在细胞群中的表达情况展示(小提琴图和气泡图)

# plot the signaling gene expression distribution
p <- plotGeneExpression(cellchat, signaling = "TGFb")
ggsave("TGFb_GeneExpression_vln.pdf",p,width = 8,height = 8)
p <- plotGeneExpression(cellchat, signaling = "TGFb", type = "dot")
ggsave("TGFb_GeneExpression_dot.pdf",p,width = 8,height = 6)

7.通讯网络系统分析(可信度有待考量)

通讯网络系统分析使用了三种算法：社会网络分析、NMF分析和流行学习与分析 注意：不同的算法算出来的结果可能互相矛盾，需要结合生物学知识加以判断。

7.1社会网络分析(通讯网络中的角色识别)

计算网络中心权重

cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")

通过计算每个细胞群的网络中心性指标，识别每类细胞在信号通路中的角色/作用C(发送者、接收者、调节者和影响者)

netAnalysis_signalingRole_network(cellchat, signaling = pathways.show, width = 15, height = 6, font.size = 10)

识别细胞的信号流模式

ht1 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing", font.size = 5)
ht2 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming", font.size = 5)
ht1 + ht2

上图横轴是细胞类型，纵轴是pathway。左图是各个细胞类型中各个通路发出信号的强度，右图是各个细胞类型中各个通路接受信号的强度

7.2非负矩阵分解(NMF)识别细胞的通讯模式(这里是一个比较标准的NMF的应用方式)
  • 信号输出细胞的模式识别
# 计算分解成几个因子(pattern)比较合适(这一步运行比较慢，在使用NMF对细胞进行亚群细分时，如果不测试的话，最好选择比细胞类型多一点的值)
selectK(cellchat, pattern = "outgoing")

热图

nPatterns = 2 # 挑选曲线中第一个出现下降的点(从2就开始下降了)
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "outgoing", k = nPatterns, width = 5, height = 9, font.size = 6)

按selectK算出来的pattern值分为了2个大的pattern，左图纵轴是细胞类型，算法自动将DC、Memory CD4T、FCGR3A+Mono和Platelet分为了一个pattern(和Pattern2对于)，剩下的五种细胞分为了一个pattern(和pattern1对应)。右图纵轴是信号通路，也是被分为了两个大的pattern，上面一个部分是在Pattern1中比较活跃的通路，下面一个部分是在Pattern2中比较活跃的通路

冲击图/河流图(river plot)

netAnalysis_river(cellchat, pattern = "outgoing")

气泡图

netAnalysis_dot(cellchat, pattern = "outgoing")

  • 信号输入细胞的模式识别
selectK(cellchat, pattern = "incoming")

热图

nPatterns = 2
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "incoming", k = nPatterns, 
                                          width = 5, height = 9, font.size = 6)

冲击图

netAnalysis_river(cellchat, pattern = "incoming")

气泡图

netAnalysis_dot(cellchat, pattern = "incoming")

7.3信号网络的流行学习与分类

把共同作用的信号通路归纳在一起，分为基于功能的归纳和基于拓扑结构的归纳

  • 基于功能相似性的流行学习与分类
cellchat <- computeNetSimilarity(cellchat, type = "functional")
cellchat <- netEmbedding(cellchat, type = "functional")
cellchat <- netClustering(cellchat, type = "functional")
p <- netVisual_embedding(cellchat, type = "functional", label.size = 3.5)
ggsave("Manifold_functional_cluster.pdf", p, width = 8, height = 6)
  • 基于拓扑相似性的流行学习与分类
cellchat <- computeNetSimilarity(cellchat, type = "structural")
cellchat <- netEmbedding(cellchat, type = "structural")
cellchat <- netClustering(cellchat, type = "structural")
p <- netVisual_embedding(cellchat, type = "structural", label.size = 3.5)
ggsave("Manifold_structural_cluster.pdf", p, width = 8, height = 6)
saveRDS(cellchat, file = "cellchat.rds")

二、不同分组之间的配对分析

注意：配对分析必须保证细胞类型是一样的，才可以进行配对。如果两个样本的细胞类型不一样又想进行配对分析时，可以用subset把两个样本的细胞类型取成一致的。

1.数据准备，分别创建CellChat对象
Sys.setenv(RETICULATE_PYTHON='/usr/bin/python3')
library(Seurat)
library(tidyverse)
library(CellChat)
library(NMF)
library(ggalluvial)
library(patchwork)

rm(list = ls());gc()
options(stringsAsFactors = F)

# 创建cellchat对象
# 提取数据子集
scRNA <- readRDS("~/project/Integrate/scRNA.classified.rds")
scRNA$celltype <- scRNA$SingleR
table(scRNA$celltype)
Idents(scRNA) <- 'celltype'
scRNA <- subset(scRNA, idents = c("B cells", "CD4+ T cells", "CD8+ T cells", "Dendritic cells", "Monocytes", "NK cells"))
scRNA$celltype <- as.factor(as.character(scRNA$celltype))
table(scRNA$orig.ident)
Idents(scRNA) <- "orig.ident"
sco.til <- subset(scRNA, idents = c("HNC01TIL", "HNC10TIL", "HNC20TIL"))
sco.pbmc <- subset(scRNA, idents = c("HNC01PBMC", "HNC10PBMC", "HNC20PBMC"))
# 创建cellchat对象
cco.til <- createCellChat(sco.til@assays$RNA@data, meta = sco.til@meta.data, group.by = "celltype")
cco.pbmc <- createCellChat(sco.pbmc@assays$RNA@data, meta = sco.pbmc@meta.data, group.by = "celltype")
save(cco.til, cco.pbmc, file = "cco.rda")
2.细胞通讯网络分析
  • 2.1数据准备和路径切换
dir.create("./Compare")
setwd("./Compare")
  • 2.2分析样本cco.pbmc的细胞通讯网络
cellchat <- cco.pbmc
cellchat@DB <- CellChatDB.human
cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
# cellchat <- projectData(cellchat, PPI.human)
cellchat <- computeCommunProb(cellchat, raw.use = T, population.size = T)
# cellchat <- filterCommunication(cellchat, min.cells = 5)
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")
# cellchat <- computeNetSimilarity(cellchat, type = "functional")
# cellchat <- netEmbedding(cellchat, type = "functional")
# cellchat <- netClustering(cellchat, type = "functional")
# cellchat <- computeNetSimilarity(cellchat, type = "structural")
# cellchat <- netEmbedding(cellchat, type = "structural")
# cellchat <- netClustering(cellchat, type = "structural")
cco.pbmc <- cellchat
saveRDS(cco.pbmc, file = "cco.pbmc.rds")
  • 2.3分析样本cco.til的细胞通讯网络
cellchat <- cco.til
cellchat@DB <- CellChatDB.human
cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
# cellchat <- projectData(cellchat, PPI.human)
cellchat <- computeCommunProb(cellchat, raw.use = T, population.size = T)
cellchat <- filterCommunication(cellchat)
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")
# cellchat <- computeNetSimilarity(cellchat, type = "functional")
# cellchat <- netEmbedding(cellchat, type = "functional")
# cellchat <- netClustering(cellchat, type = "functional")
# cellchat <- computeNetSimilarity(cellchat, type = "structural")
# cellchat <- netEmbedding(cellchat, type = "structural")
# cellchat <- netClustering(cellchat, type = "structural")
cco.til <- cellchat
saveRDS(cco.til, file = "cco.til.rds")
  • 2.4合并cellchat对象
cco.list <- list(pbmc = cco.pbmc, til = cco.list)
cellchat <- mergeCellChat(cco.list, add.names = names(cco.list), cell.prefix = T)
3.可视化
3.1.所有细胞群总体观：通讯数量与强度对比
gg1 <- compareInteractions(cellchat, show.legend = F, group = c(1,2), measure = "count")
gg2 <- compareInteractions(cellchat, show.legend = F, group = c(1,2), measure = "weight")
p <- gg1 + gg2
ggsave("Overview_number_strength.pdf", p, width = 6, height = 4)

左图展示通讯数量之间的差异，右图展示通讯强度之间的差异。本例中信号通路强度weight值过低，导致显示时均为0

数量与强度差异网络图

par(mfrow = c(1,2))
netVisual_diffInteraction(cellchat, weight.scale = T)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "weight")

数量与强度差异热图

par(mfrow = c(1,1))
h1 <- netVisual_heatmap(cellchat)
h2 <- netVisual_heatmap(cellchat, measure = "weight")
h1 + h2

case和control对比，红色是上调，蓝色是下调

细胞互作数量对比网络图

par(mfrow = c(1,2))
weight.max <- getMaxWeight(cco.list, attribute = c("idents", "count"))
for(i in 1:length(cco.list)){
  netVisual_circle(cco.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(cco.list)[i]))
}

左图是control，右图是case，可以直接对比数量变化

3.2.制定细胞互作数量对比网络图
par(mfrow = c(1,2))
s.cell <- c("CD4+ T cells", "CD8+ T cells", "Monocytes")
count1 <- cco.list[[1]]@net$count[s.cell, s.cell]
count2 <- cco.list[[2]]@net$count[s.cell, s.cell]
weight.max <- max(max(count1), max(count2))
netVisual_circle(count1, weight.scale = T, label.edge = T, edge.weight.max = weight.max, edge.width.max = 12,
                 title.name = paste0("Number of interactions-", names(cco.list)[1]))
netVisual_circle(count2, weight.scale = T, label.edge = T, edge.weight.max = weight.max, edge.width.max = 12,
                 title.name = paste0("Number of interactions-", names(cco.list)[2]))

3.3.保守和特异性信号通路的识别与可视化
# 通路信号强度对比分析
gg1 <- rankNet(cellchat, mode = "comparison", stacked = T, do.stat = T)
gg2 <- rankNet(cellchat, mode = "comparison", stacked = T, do.stat = T)
p <- gg1 + gg2
ggsave("Compare_pathway_strength.pdf",p,width = 10,height = 6)

3.4.流行学习识别差异信号通路
cellchat <- computeNetSimilarityPairwise(cellchat, type = "functional")
cellchat <- netEmbedding(cellchat, type = "functional")
cellchat <- netClustering(cellchat, type = "functional")
# netVisual_embeddingPairwise(cellchat, type = "functional", label.size = 3.5)
# netVisual_embeddingPairwiseZoomIn(cellchat, type = "functional", nCol = 2)
cellchat <- computeNetSimilarityPairwise(cellchat, type = "structural")
cellchat <- netEmbedding(cellchat, type = "structural")
cellchat <- netClustering(cellchat, type = "structural")
netVisual_embeddingPairwise(cellchat, type = "structural", label.size = 3.5)
netVisual_embeddingPairwiseZoomIn(cellchat, type = "structural", nCol = 2)
p <- rankSimilarity(cellchat, type = "structural") + ggtitle("structural similarity of pathway")
ggsave("Pathway_Similarity.pdf", p, width = 8, height = 5)

saveRDS(cellchat, file = "cellchat.rds")

case和control之间信号通路差异相差成度排行，在这张图中，ICAM相差最大，其次是SELPLG。

3.5.细胞信号模式对比
library(ComplexHeatmap)

总体信号模式对比

pathway.union <- union(cco.list[[1]]@netP$pathways, cco.list[[2]]@netP$pathways)
ht1 <- netAnalysis_signalingRole_heatmap(cco.list[[1]], pattern = "all", signaling = pathway.union,
                                         title = names(cco.list)[1], width = 8, height = 10)
ht2 <- netAnalysis_signalingRole_heatmap(cco.list[[2]], pattern = "all", signaling = pathway.union,
                                         title = names(cco.list)[2], width = 8, height = 10)
draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))

左图下面几个信号通路在PBMC中没有，和3.3中的图相符

输出信号模式对比

pathway.union <- union(cco.list[[1]]@netP$pathways, cco.list[[2]]@netP$pathways)
ht1 <- netAnalysis_signalingRole_heatmap(cco.list[[1]], pattern = "outgoing", signaling = pathway.union,
                                         title = names(cco.list)[1], width = 8, height = 10)
ht2 <- netAnalysis_signalingRole_heatmap(cco.list[[2]], pattern = "outgoing", signaling = pathway.union,
                                         title = names(cco.list)[2], width = 8, height = 10)
draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))

输入信号模式对比

pathway.union <- union(cco.list[[1]]@netP$pathways, cco.list[[2]]@netP$pathways)
ht1 <- netAnalysis_signalingRole_heatmap(cco.list[[1]], pattern = "incoming", signaling = pathway.union,
                                         title = names(cco.list)[1], width = 8, height = 10)
ht2 <- netAnalysis_signalingRole_heatmap(cco.list[[2]], pattern = "incoming", signaling = pathway.union,
                                         title = names(cco.list)[2], width = 8, height = 10)
draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))

3.6.特定信号通路的对比

网络图

pathways.show <- c("IL16")
weight.max <- getMaxWeight(cco.list, slot.name = c("netP"), attribute = pathways.show)
par(mfrow = c(1,2), xpd = TRUE)
for(i in 1:length(cco.list)){
  netVisual_aggregate(cco.list[[i]], signaling = pathways.show, layout = "circle", edge.weight.max = weight.max[1],
                      edge.width.max = 10, signaling.name = paste(pathways.show, names(cco.list)[i]))
}

热图

par(mfrow = c(1,2), xpd = TRUE)
ht <- list()
for(i in 1:length(cco.list)){
  ht[[i]] <- netVisual_heatmap(cco.list[[i]], signaling = pathways.show, color.heatmap = "Reds",
                               title.name = paste(pathways.show, "signaling ", names(cco.list)[i]))
}
ComplexHeatmap::draw(ht[[1]] + ht[[2]], ht_gap = unit(0.5, "cm"))

和弦图

par(mfrow = c(1,2), xpd = T)
for(i in 1:length(cco.list)){
  netVisual_aggregate(cco.list[[i]], signaling = pathways.show, layout = "chord", pt.title = 3,
                      title.space = 0.05, vertex.label.cex = 0.6, signaling.name = paste(pathways.show, names(cco.list)[i]))
}
3.7.配体-受体对比分析

气泡图展示所有配体受体对的差异

levels(cellchat@idents$joint)
p <- netVisual_bubble(cellchat, sources.use = c(4,5), targets.use = c(1,2,3,6), comparison = c(1,2), angle.x = 45)
ggsave("Compare_LR_bubble.pdf", p, width = 12, height = 8)

气泡图展示上调或下调的配体受体对

p1 <- netVisual_bubble(cellchat, sources.use = c(4,5), targets.use = c(1,2,3,6), comparison = c(1,2),
                       max.dataset = 2, title.name = "Increased signaling in TIL", angle.x = 45, remove.isolate = T)
p2 <- netVisual_bubble(cellchat, sources.use = c(4,5), targets.use = c(1,2,3,6), comparison = c(1,2),
                       max.dataset = 1, title.name = "Decreased signaling in TIL", angle.x = 45, remove.isolate = T)
pc <- p1 + p2
ggsave("Compare_LR_regulated.pdf", pc, width = 12, height = 5.5)

和弦图

par(mfrow = c(1,2), xpd = T)
for(i in 1:length(cco.list)){
  netVisual_chord_gene(cco.list[[i]], sources.use = c(4,5), targets.use = c(1,2,3,6), signaling = "MHC-I",
                       lab.cex = 0.6, legend.pos.x = 10, legend.pos.y = 20,
                       title.name = paste0("Signaling from Treg - ", names(cco.list)[i]))
}