不同分组之间的配对分析
⚠️：配对分析必须保证细胞类型是一样的，才可以进行配对。如果 两个样本的细胞类型不一样又想进行配对分析时，可以用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组间分析结果影响较大，尤其是细胞互作强度。
如处理组和对照组相比，巨噬细胞显著减少，单核细胞显著增多的情况下，蛋壳图中巨噬细胞和巨噬细胞的互作强度（注意 不是数量）就会显著下降，单核和单核的互作强度就会显著上升。
由于互作强度是根据表达的受体配体数决定的，那么细胞占比增多，强度就会增加，细胞占比下降，强度就下降？这样似乎不是很科学
![【移动通讯】【MIMO】[P1]【科普篇】](https://img-blog.csdnimg.cn/direct/887b630c28db4687b1dd88d329243248.png)
