细胞通讯-5 实操
2022-08-25 本文已影响0人
oceanandshore
这篇重点写下实操过程中一些注意事项。
视频教程里面说用的是 human skin ，但是没说用哪个。
从Seurat对象直接创建：
⚠️：构建Cell Chat对象时，输入的是log后的数据。
1、单样本
# Chord diagram
group.cellType <- c(rep("FIB", 4), rep("DC", 4), rep("TC", 4)) # grouping cell clusters into fibroblast, DC and TC cells
names(group.cellType) <- levels(cellchat@idents)
par(mfrow = c(1,2), xpd=TRUE)
netVisual_chord_cell(cellchat, signaling = pathways.show[1],
group = group.cellType, #分组展示的顺序
title.name = paste0(pathways.show[1], " signaling network"))
## Plot the aggregated cell-cell communication network at the signaling pathway level
netVisual_chord_cell(cellchat, signaling = pathways.show[1],
title.name = paste0(pathways.show[1], " signaling network"))
## Plot the aggregated cell-cell communication network at the signaling pathway level
# p1加了分组顺序，p2没加，可以看出，细胞分类的话，类与类之间有个小间隔，更容易分辨
配受体展示
#####配受体展示
p1 <- netAnalysis_contribution(cellchat, signaling = pathways.show[1],
title = pathways.show[1])#展现对特定通路的贡献程度
p2 <- netAnalysis_contribution(cellchat, signaling = pathways.show)
cowplot::plot_grid(p1, p2, align = "h",ncol=2)
指定通路画连线图
#指定通路画连线图
pairLR.CXCL <- extractEnrichedLR(cellchat, signaling = pathways.show[1],
geneLR.return = FALSE) #提取显著作用的配受体对
LR.show <- pairLR.CXCL[1,] # show one ligand-receptor pair
vertex.receiver = seq(1,4) # a numeric vector 取1到4的时候，source就只有4种细胞，其实也就是取多少个levels作为source
p1<-netVisual_individual(cellchat, signaling = pathways.show,
pairLR.use = LR.show, vertex.receiver = vertex.receiver,
layout = 'hierarchy')
vertex.receiver = seq(1,7) # a numeric vector 取1到7的时候，source就只有7种细胞
p2<-netVisual_individual(cellchat, signaling = pathways.show,
pairLR.use = LR.show, vertex.receiver = vertex.receiver,
layout = 'hierarchy')
vertex.receiver = seq(1,4) 即source就只有4种细胞
vertex.receiver = seq(1,7) 即source就只有7种细胞
netVisual_individual(cellchat, signaling = pathways.show[1],
pairLR.use = LR.show, layout = "circle")
# 取的是LR.show，也就是最显著的一个配受体，所以这个函数可以指定配受体
#LR.show
#[1] "FGF7_FGFR1"
image.png
多个细胞通讯通路气泡图可视化
#####多个细胞通讯通路气泡图可视化
#指定发送者与接收者，展示所有配受体对
levels(cellchat@idents)
## [1] "APOE+ FIB" "FBN1+ FIB" "COL11A1+ FIB" "Inflam. FIB" "cDC1"
## [6] "cDC2" "LC" "Inflam. DC" "TC" "Inflam. TC"
## [11] "CD40LG+ TC" "NKT"
par(mfrow = c(1,3), xpd=TRUE)
p1 <- netVisual_bubble(cellchat, sources.use = 4,
targets.use = c(5:11), remove.isolate = FALSE) # targets.use 只展示想要的目标
## Comparing communications on a single object
p2 <- netVisual_bubble(cellchat, sources.use = 4,
targets.use = c("cDC1","cDC2","LC","Inflam. DC","TC","Inflam. TC",
"CD40LG+ TC"), remove.isolate = FALSE)
## Comparing communications on a single object
p3 <-netVisual_bubble(cellchat, sources.use = 4,
targets.use = c("cDC1","cDC2"), remove.isolate = FALSE)
## Comparing communications on a single object
#指定信号通路
p4 <-netVisual_bubble(cellchat, sources.use = 4,
targets.use = c(5:11), signaling = c("CCL","CXCL"),
remove.isolate = FALSE)
## Comparing communications on a single object
cowplot::plot_grid(p1, p2,p3,p4, align = "h",ncol=4)
#纵坐标展示的是配受体的名称；p1和p2两张图是一模一样的，一个是数字，一个是用标签标识，结果都一样 是 Rplot13-气泡图
Rplot13-气泡图
指定信号通路提取配受体
#指定信号通路提取配受体
par(mfrow = c(1,3), xpd=TRUE)
pairLR.use <- extractEnrichedLR(cellchat, signaling = c("CCL","CXCL","FGF"))
pairLR.use
## interaction_name
## 1 CCL19_CCR7
## 2 CXCL12_CXCR4
## 3 CXCL12_ACKR3
## 4 FGF7_FGFR1
#也即是说"CCL","CXCL","FGF"这3个通路当中，有4对显著的配受体
p1 <-netVisual_bubble(cellchat, sources.use = c(3,4), targets.use = c(5:8),
pairLR.use = pairLR.use, remove.isolate = TRUE) #remove.isolate = TRUE 去除空的通路
## Comparing communications on a single object
p2 <-netVisual_bubble(cellchat, sources.use = 4, targets.use = c(5:11),
signaling = c("CCL","CXCL"), remove.isolate = FALSE)
## Comparing communications on a single object
p3 <-netVisual_bubble(cellchat, sources.use = 4, targets.use = c(5:11),
signaling = c("CCL","CXCL"), remove.isolate = T)
## Comparing communications on a single object
cowplot::plot_grid(p1, p2,p3, align = "h",ncol=3)
#p1是指定了配受体对（pairLR.use），并且指定了2种细胞sources.use = c(3,4)；
#p2指定的是通路的名称，指定了一种细胞sources.use = 4；
#p3里面用了targets.use = c(5:11)
#> levels(meta$labels)[5:11]
#[1] "cDC1" "cDC2" "LC" "Inflam. DC" "TC" "Inflam. TC" "CD40LG+ TC"
Rplot14
2、多样本
这部分的内容主要分析目的为：
1、细胞通讯在不同组别中是否发生变化
2、细胞通讯在不同细胞类型中是否发生变化
3、细胞通讯的“发起者”与接收者是否随组别发生变化
首先还是做点准备工作
library(CellChat)
library(patchwork)
library(cowplot)
dir.create('./comparison')
setwd('./comparison')
#加载两个数据集
#这个数据集相当于我们上次课程中已经进行细胞通讯计算的cellchat对象，没有这部分基础的同学请回顾：https://www.bilibili.com/video/BV1Ab4y1W7qx?p=3
#cellchat.NL <- readRDS(url("https://ndownloader.figshare.com/files/25954199"))
#cellchat.LS <- readRDS(url("https://ndownloader.figshare.com/files/25956518"))
#dir.create('data')
#saveRDS(cellchat.LS,'./data/cellchat.LS.rds')
#saveRDS(cellchat.NL,'./data/cellchat.NL.rds')
cellchat.LS <- readRDS('./data/cellchat.LS.rds')
cellchat.NL <- readRDS('./data/cellchat.NL.rds')
object.list <- list(NL = cellchat.NL, LS = cellchat.LS) #构建整合list
cellchat <- mergeCellChat(object.list, add.names = names(object.list)) #进行合并，相当于seurat里面的merge函数
## 合并的数据槽：
##Merge the following slots: 'data.signaling','net', 'netP','meta', 'idents', 'var.features' , 'DB', and 'LR'.
cellchat
## An object of class CellChat created from a merged object with multiple datasets
## 555 signaling genes.
## 7563 cells. # 这里的7563个细胞就是合并前加起来的和
#最简单的展示，查看细胞互作的数量在不同条件下是否有差异
gg1 <- compareInteractions(cellchat, show.legend = F, group = c(1,2))
gg2 <- compareInteractions(cellchat, show.legend = F, group = c(1,2), measure = "weight")
gg1 + gg2
左图是细胞通讯数量，右边是细胞通讯的权重。
image.png
#查看细胞通路在两组间的富集程度
gg1 <- rankNet(cellchat, mode = "comparison", stacked = T, do.stat = TRUE)
gg2 <- rankNet(cellchat, mode = "comparison", stacked = F, do.stat = TRUE)
gg1 + gg2
image.png
#以circle plot的形式展示第二个组别中相较于第一个组别细胞通讯发生的变化，红色为上调蓝色为下调
par(mfrow = c(1,2), xpd=TRUE)
netVisual_diffInteraction(cellchat, weight.scale = T)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "weight")
#不是以分组拆开的方式去展示，而是直接以两组数据相减（第二组减去第一组）得到的结果展示
#红色代表在第二个组别当中上调，蓝色代表在第二个组别当中这些通路下调
#如何看哪个是第一个组别，哪个是第二个组别，用names ( object.list)看 NL是第一组
#> names ( object.list)
#[1] "NL" "LS"
image.png
#以上是直接展示二组“相减”的结果，当然你也可以直接将两组分开展示:
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]))
}
image.png
#同理，用heatmap也可以进行展示
gg1 <- netVisual_heatmap(cellchat)
## Do heatmap based on a merged object
gg2 <- netVisual_heatmap(cellchat, measure = "weight")
## Do heatmap based on a merged object
gg1 + gg2
#同样是第二个组减去第一个组，左侧是细胞通讯数量，右侧是细胞通讯强度，和上面的圈图一样的。
image.png
#展示特定通路， 老三样
#cicle
pathways.show <- c("CXCL")
weight.max <- getMaxWeight(object.list, slot.name = c("netP"), attribute = pathways.show) # control the edge weights across different datasets
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]))
}
# 展示CXCL通路在组间的差异，可以看到在LS组有更多的细胞参与这个通路
image.png
#heatmap
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]))
}
## Do heatmap based on a single object
##
## Do heatmap based on a single object
ComplexHeatmap::draw(ht[[1]] + ht[[2]], ht_gap = unit(0.5, "cm"))
#Rplot7
image.png
#弦图
pathways.show <- c("CXCL")
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
netVisual_aggregate(object.list[[i]], signaling = pathways.show,
layout = "chord",
signaling.name = paste(pathways.show, names(object.list)[i]))
}
image.png
#基因表达情况
cellchat@meta$datasets = factor(cellchat@meta$datasets,
levels = c("NL", "LS")#设定组别的level
)
plotGeneExpression(cellchat, signaling = "CXCL", split.by = "datasets", colors.ggplot = T)
## Registered S3 method overwritten by 'spatstat.geom':
## method from
## plot.imlist imager
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##
## This message will be shown once per session.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
image.png
#研究细胞类型之间的通讯差异
#按照组别拆分整合
levels(object.list[[1]]@idents)
## [1] "APOE+ FIB" "FBN1+ FIB" "COL11A1+ FIB" "Inflam. FIB" "cDC1"
## [6] "cDC2" "LC" "Inflam. DC" "TC" "Inflam. TC"
## [11] "CD40LG+ TC" "NKT"
levels(object.list[[2]]@idents)
## [1] "APOE+ FIB" "FBN1+ FIB" "COL11A1+ FIB" "Inflam. FIB" "cDC1"
## [6] "cDC2" "LC" "Inflam. DC" "TC" "Inflam. TC"
## [11] "CD40LG+ TC" "NKT"
group.cellType <- c(rep("FIB", 4), rep("DC", 4), rep("TC", 4))
group.cellType <- factor(group.cellType, levels = c("FIB", "DC", "TC"))
object.list <- lapply(object.list, function(x) {mergeInteractions(x, group.cellType)})#合并细胞类型
cellchat <- mergeCellChat(object.list, add.names = names(object.list))#合并两组数据
## Merge the following slots: 'data.signaling','net', 'netP','meta', 'idents', 'var.features' , 'DB', and 'LR'.
#直接展示相减的结果
par(mfrow = c(1,2), xpd=TRUE)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "count.merged", label.edge = T)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "weight.merged", label.edge = T)
# 这里还是两组数据相减的结果，把12种细胞类型汇聚成3大类。DC相较其他两类，变化最大
image.png
#两组分别展示
weight.max <- getMaxWeight(object.list, slot.name = c("idents", "net", "net"), attribute = c("idents","count", "count.merged"))
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
netVisual_circle(object.list[[i]]@net$count.merged, weight.scale = T, label.edge= T, edge.weight.max = weight.max[3], edge.width.max = 12, title.name = paste0("Number of interactions - ", names(object.list)[i]))
}
#Rplot11 DC与TC之间的通讯，显然LS组更强一些
image.png
#展示两组间各类细胞incoming与outcoming通讯的强度
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)) # control the dot size in the different datasets
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)
}
## Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways
## Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways
patchwork::wrap_plots(plots = gg)
#Rplot12
image.png
展示两组间的差异
报错，视频说需要环境相同，视频在17：30秒处
#展示两组间的差异
gg1 <- netAnalysis_signalingChanges_scatter(cellchat,
idents.use = "Inflam. DC",
signaling.exclude = "MIF")
## Visualizing differential outgoing and incoming signaling changes from NL to LS
gg2 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "cDC1",
signaling.exclude = c("MIF"))
## Visualizing differential outgoing and incoming signaling changes from NL to LS
patchwork::wrap_plots(plots = list(gg1,gg2))
#报错，视频说需要环境相同，视频在17：30秒处
#继续探索outgoing与incoming的组间模式差异
suppressMessages(library(ComplexHeatmap))
i = 1
#outgoing
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"))
# cDC1 在两组中差异较大
image.png
#incoming
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"))
#Rplot14
image.png
#overall
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"))
#Rplot15
image.png
#从配受体的角度画一画气泡图
netVisual_bubble(cellchat, sources.use = 4,
targets.use = c(5:11), comparison = c(1, 2), angle.x = 45)
## Comparing communications on a merged object
#Rplot16
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)
## Comparing communications on a merged object
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)
## Comparing communications on a merged object
gg1 + gg2
#Rplot17
16
17
#以上的计算依赖的都是细胞通讯的可能性(强度)，接下来我们将通过配受体对基因表达的角度进一步研究
#以通路为单位
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)
## Use the joint cell labels from the merged CellChat object
#提取细胞通讯预测数据框
net <- netMappingDEG(cellchat, features.name = features.name)
#提取LS组上调的配体以及NL上调的受体，
net.up <- subsetCommunication(cellchat, net = net, datasets = "LS",
ligand.logFC = 0.2, receptor.logFC = NULL)
#反之亦然
net.down <- subsetCommunication(cellchat, net = net, datasets = "LS",
ligand.logFC = -0.2, receptor.logFC = -0.1)
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]))
## Comparing communications on a merged object
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]))
## Comparing communications on a merged object
gg1 + gg2
#Rplot18
image.png
#弦图
par(mfrow = c(1,2), xpd=TRUE)
netVisual_chord_gene(object.list[[1]], sources.use = 4, targets.use = c(5:11),
slot.name = 'net', net = net.down, lab.cex = 0.8, small.gap = 3.5,
title.name = paste0("Down-regulated signaling in ",
names(object.list)[1]))
## You may try the function `netVisual_chord_cell` for visualizing individual signaling pathway
netVisual_chord_gene(object.list[[2]], sources.use = 4, targets.use = c(5:11),
slot.name = 'net', net = net.up, lab.cex = 0.8, small.gap = 3.5,
title.name = paste0("Up-regulated signaling in ",
names(object.list)[2]))
#Rplot19
image.png
#基因表达情况
cellchat@meta$datasets = factor(cellchat@meta$datasets,
levels = c("NL", "LS")#设定组别的level
)
plotGeneExpression(cellchat, signaling = "CXCL", split.by = "datasets", colors.ggplot = T)
#这部分结果展示建议用Seurat做
