1. 数据的准备

#加载需要的包
library(CellChat)
library(patchwork)

#创建图片保存的目录
data.dir <- './Result_CellChat'
dir.create(data.dir)
setwd(data.dir)

#加载对应数据集的cellchat对象
cellchat.NL <- readRDS(url("https://ndownloader.figshare.com/files/25954199"))
cellchat.LS <- readRDS(url("https://ndownloader.figshare.com/files/25956518"))

#合并数据集
object.list <- list(NL = cellchat.NL, LS = cellchat.LS)

cellchat <- mergeCellChat(object.list, add.names = names(object.list))
#需要在每个数据集上单独运行 CellChat，然后合并对象。如果使用较早版本（< 0.5.0）获得的 CellChat 对象，先updateCellChat

2.预测细胞通信

在比较多种生物条件下的细胞-细胞通信时，可以得知细胞通信是否增强；细胞类型显著变化之间的相互作用；不同条件下，source和target的变化

2.1 比较交互总数和交互强度

gg1 <- compareInteractions(cellchat, show.legend = F, group = c(1,2))
gg2 <- compareInteractions(cellchat, show.legend = F, group = c(1,2), measure = "weight")
gg1 + gg2

2.2 比较不同细胞亚群之间的相互作用数量和强度

边的宽度表示交互的相对数量或交互强度。与第一个数据集相比，第二个数据集中的红色（或蓝色）边缘表示增加（或减少）的信号。

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

可以使用热图在更大的细节中显示交互的差异数或交互强度。顶部彩色条形图表示热图（传入信号）中显示的列值的总和。右边的彩色条形图表示一行值（传出信号）的总和。在色条中红色或蓝色表示第二个数据集中与第一个数据集相比增加或[减少]信号。

gg1 <- netVisual_heatmap(cellchat, measure = "count")

gg2 <- netVisual_heatmap(cellchat, measure = "weight")

gg1 + gg2
#measure   "count" or "weight". "count": comparing the number of interactions; "weight": comparing the total interaction weights (strength)

差异网络分析仅适用于配对数据集。如果有更多的数据集进行比较，直接显示每个数据集中任意两个细胞群之间的交互次数或交互强度。

为了更好地控制不同数据集中推断网络的节点大小和边缘权重，计算每个细胞组的最大细胞数量以及所有数据集中交互（或交互权重）的最大数量。

weight.max <- getMaxWeight(object.list, attribute = c("idents","count"))
par(mfrow = c(1,2), xpd=TRUE)
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]))
}

2.3 比较 2D 空间中的主要source和target

比较二D空间中的传出和传入交互强度，可以识别不同数据集之间显著变化的发送或接收信号的细胞群。

num.link <- sapply(object.list, function(x) {rowSums(x@net$count) + colSums(x@net$count)-diag(x@net$count)})

weight.MinMax <- c(min(num.link), max(num.link)) 
# 控制不同数据集中的气泡大小
gg <- list()
for (i in 1:length(object.list)) {
  gg[[i]] <- netAnalysis_signalingRole_scatter(object.list[[i]], title = names(object.list)[i], weight.MinMax = weight.MinMax)
}
#从所有信号通路分析聚合细胞-细胞通信网络的信号作用
patchwork::wrap_plots(plots = gg)

从散点图中，可以看到与 NL 相比，Inflam.DC 和 cDC1 成为 LS 中的主要source和target之一。纤维细胞群也成为LS的主要source。

3.识别保守和环境特异的信号通路

CellChat 可以根据对象在多种生物条件下的细胞通信网络，识别差异较大（或更少）的信号网络、信号组以及基于其细胞通信网络的保守和环境特异的信号通路。

3.1 根据信号/结构的相似性识别差异较大（或更少）的信号网络以及信号组

CellChat 根据推断的通信网络的功能和拓扑相似性，对其进行联合多重学习和分类。NB：此类分析适用于两个以上的数据集。

功能相似性：功能相似度高表示主要source和target相似，可解释为两个信号通路或两个配体受体对具有相似的作用。功能相似性分析不适用于具有不同细胞类型成分的多个数据集。

结构相似性：结构相似性用于比较其信号网络结构，而不考虑source和target的相似性。结构相似性分析适用于具有相同细胞类型组成或截然不同的细胞类型组成多个数据集。

在这里，我们可以根据功能相似性运行多重和分类学习分析，因为两个数据集具有相同的单元类型组成。

根据信号组的功能相似性识别信号组

cellchat <- computeNetSimilarityPairwise(cellchat, type = "functional")
#计算数据集信号网络相似性
cellchat <- netEmbedding(cellchat, type = "functional")
#数据集网络的多重学习
cellchat <- netClustering(cellchat, type = "functional")
#数据集网络的分类
# 2D 可视化
netVisual_embeddingPairwise(cellchat, type = "functional", label.size = 3.5)
# netVisual_embeddingZoomIn(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)

计算和可视化通路距离

我们可以根据信号网络在共享双维空间中的欧几里德距离来识别差异较大（或更少）的信号网络。更大的距离意味着两个数据集之间的通信网络在功能或结构相似性方面存在更大的差异。只计算两个数据集之间重叠信号通路的距离。此处未考虑仅在一个数据集中标识的信号通路。如果有三个以上的数据集，可以通过在函数rankSimilarity中定义comparison进行对比。

#计算欧几里德距离
rankSimilarity(cellchat, type = "functional")

3.2 识别并可视化保守和环境特异的信号通路

通过比较每个信号通路的信息流/交互，可以识别信号通路的状态：关闭，减少，打开，增加。
比较每个信号通路的整体信息流

可以通过简单地比较每个信号通路的信息流来识别保守和环境特异的信号通路，该信息流由推断网络中所有一对细胞群之间的通信概率之和（即网络中的总权重）定义。

此条形图可在堆叠模式下绘制。根据 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

比较与每个细胞群相关的传出（或传入）信号

上述分析将传出和传入信号的信息汇总在一起。还可以比较两个数据集之间的传出（或传入）信号模式，从而识别显示不同信号模式的信号通路/受配体。

可以将来自不同数据集的所有已识别的信号通路进行组合，从而并排比较它们，包括传出信号、传入信号和整体信号，方法是将传出和传入信号聚合在一起。rankNet还显示了整体信号的比较，但它没有显示特定细胞群中的信号强度。

library(ComplexHeatmap)

i = 1
# 组合来自不同数据集的所有已识别的信号通路 
pathway.union <- union(object.list[[i]]@netP$pathways, object.list[[i+1]]@netP$pathways)
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"))

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

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

4. 识别上调和下调的信号配体对

可以比较由某些细胞群到其他细胞亚群的配体受体对调节的通信概率。可以通过设置comparison在函数netVisual_bubble中来完成。

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

气泡图中显示的配体受体对可以通过signaling.LSIncreased = gg1$data 访问。

通过比较每个 L-R 对和每对细胞组的两个数据集之间的通信概率，可以采用上述方法来识别上调和下调的信号。另外，可以根据微分基因表达分析来识别上调和下调的信号配体对。具体来说，对每个细胞组执行两种生物条件（即NL和LS）之间的微分表达分析，然后根据发送者细胞中配体和接收器细胞中受体的折叠变化获得上调和下调的信号。此类分析可如下所示。

# 定义一个阳性数据集
pos.dataset = "LS"
# 定义一个字符名称，用于存储差分表达式分析的结果
features.name = pos.dataset
# 进行差异表达分析
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)

#将差异表达分析的结果映射到推断的细胞-细胞通信上
net <- netMappingDEG(cellchat, features.name = features.name)
# 提取LS中上调的受配体对
net.up <- subsetCommunication(cellchat, net = net, datasets = "LS",ligand.logFC = 0.2, receptor.logFC = NULL)
# 提取NS中上调的受配体对
net.down <- subsetCommunication(cellchat, net = net, datasets = "NL",ligand.logFC = -0.1, receptor.logFC = -0.1)

由于信号基因在多细胞中可能很复杂，我们可以使用net.upnet.down进一步的来获得单个信号基因。

gene.up <- extractGeneSubsetFromPair(net.up, cellchat)
gene.down <- extractGeneSubsetFromPair(net.down, cellchat)

然后，使用气泡图或和弦图可视化上调和向下调的信号配体对。

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 = 90, remove.isolate = T,title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))

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 = 90, remove.isolate = T,title.name = paste0("Down-regulated signaling in ", names(object.list)[2]))

gg1 + gg2

5. 使用层次结构图、圆图或和弦图可视比较细胞-细胞通信

同单个数据集可视化类似，可以使用层次结构图、圆图或和弦图可视化细胞通信网络。

边缘颜色/重量、节点颜色/大小/形状：在所有可视化图中，边缘颜色与发送者源一致，边缘权重与交互强度成正比。较厚的边缘线表示信号更强。

在层次结构图和圆图中，圆的大小与每个细胞组中的细胞数量成正比。在层次图中，实心和开放的圆分别代表源和目标。

在和弦图中，内条颜色表示从相应的外条接收信号的目标。内条大小与目标接收的信号强度成正比。这种内条有助于解释复杂的和弦图。

展示某个通路的细胞通信

pathways.show <- c("CXCL") 
weight.max <- getMaxWeight(object.list, slot.name = c("netP"), attribute = pathways.show) 
# 控制不同数据集的边的权重
par(mfrow = c(1,2), xpd=TRUE)
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()
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]))
}

ComplexHeatmap::draw(ht[[1]] + ht[[2]], ht_gap = unit(0.5, "cm"))

6. 比较不同数据集之间的信号基因表达分布

同样使用seurat包装的函数plotGeneExpression绘制与信号通路相关的信号基因的基因表达分布图。

cellchat@meta$datasets = factor(cellchat@meta$datasets, levels = c("NL", "LS")) # set factor level
plotGeneExpression(cellchat, signaling = "CXCL", split.by = "datasets", colors.ggplot = T)

保存合并的cellchat对象

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

参考来源

CellChat三部曲2：使用CellChat 对多个数据集细胞通讯进行比较分析

致谢

I thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

THE END