多组别cellchat

不同分组之间的配对分析

??：配对分析必须保证细胞类型是一样的，才可以进行配对。如果 两个样本的细胞类型不一样又想进行配对分析时，可以用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())
options(stringsAsFactors = FALSE)

## 创建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")
# load("../cco.rda")
# cco.pbmc <- setIdent(cco.pbmc, ident.use = "celltype")
# cco.til <- setIdent(cco.til, ident.use = "celltype")

• 2.2 分析样本cco.pbmc的细胞通讯网络
  ??：cellchat不太稳定，identifyOverExpressedGenes常出错（不出现进度条就是出错了），重启Rstudio后再运算。

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 = TRUE, population.size = TRUE)
#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, "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 = TRUE, population.size = TRUE)
#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.til <- cellchat
saveRDS(cco.til, "cco.til.rds")

• 2.4 合并cellchat对象

cco.list <- list(pbmc=cco.pbmc, til=cco.til)
cellchat <- mergeCellChat(cco.list, add.names = names(cco.list), cell.prefix = TRUE)

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（实际上有数值的，只是过小，显示为0）

数量与强度差异网络图

par(mfrow = c(1,2))
netVisual_diffInteraction(cellchat, weight.scale = T)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "weight")
# save as Diff_number_strength_net.pdf

红色是case相对于control上调的，蓝色是下调的。

数量与强度差异热图

par(mfrow = c(1,1))
h1 <- netVisual_heatmap(cellchat)
h2 <- netVisual_heatmap(cellchat, measure = "weight")
h1+h2
# save as Diff_number_strength_heatmap.pdf

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]))
}
# save as Counts_Compare_net.pdf

左图是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]))
# save as Counts_Compare_select.pdf 10*6.5

3.3 保守和特异性信号通路的识别与可视化

## 通路信号强度对比分析
gg1 <- rankNet(cellchat, mode = "comparison", stacked = T, do.stat = TRUE)
gg2 <- rankNet(cellchat, mode = "comparison", stacked = F, do.stat = TRUE)
p <- gg1 + gg2
ggsave("Compare_pathway_strengh.pdf", p, width = 10, height = 6)

左图最下面5个信号通路是case组独有的

3.4 流行学习识别差异信号通路

这里function的图出不来，只有structural的图可以出来

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, "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"))
# save as Compare_signal_pattern_all.pdf  10*6

左图下面几个信号通路在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"))
# save as Compare_signal_pattern_outgoing.pdf  10*6

输入信号模式对比

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"))
# save as Compare_signal_pattern_incoming.pdf  10*6

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]))
}
# save as Compare_IL16_net.pdf  10*6.5

热图

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"))
# save as Compare_IL16_heatmap.pdf  12*6.5

和弦图

par(mfrow = c(1,2), xpd=TRUE)
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]))
}
# save as Compare_IL16_chord.pdf  10*6.5

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)

case和control配对的图，文章中常见。cellphoneDB找到的结果，也可以用这种方式呈现。

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

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=TRUE)
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]))
}
# save as Compare_LR_chord.pdf  10*6.5

小思考：但组间细胞类型占比变化似乎对CellChat组间分析结果影响较大，尤其是细胞互作强度。
如处理组和对照组相比，巨噬细胞显著减少，单核细胞显著增多的情况下，蛋壳图中巨噬细胞和巨噬细胞的互作强度（注意 不是数量）就会显著下降，单核和单核的互作强度就会显著上升。
由于互作强度是根据表达的受体配体数决定的，那么细胞占比增多，强度就会增加，细胞占比下降，强度就下降？这样似乎不是很科学