细胞通讯-5 实操

作者: oceanandshore | 来源:发表于2022-08-25 22:13 被阅读0次

细胞通讯-5 实操
细胞通讯-5 实操完整代码
【单细胞高级绘图】08.细胞通讯_两组比较_气泡图
Supercell代码实操（二）：COVID-19数据的整合分析
【王兴国营养特训班】第5期第11课7组虞文施
109-小days6第三课#裂变增长实验室#
宫氏脑针29期公告
monocle实操+细胞功能富集分析
单细胞分析：质控实操（五）
Scissor实操（二）：逻辑回归分析流程解读

这篇重点写下实操过程中一些注意事项。

视频教程里面说用的是 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

#以上的计算依赖的都是细胞通讯的可能性(强度)，接下来我们将通过配受体对基因表达的角度进一步研究
#以通路为单位
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做

细胞通讯-5 实操
这篇重点写下实操过程中一些注意事项。视频教程里面说用的是 human skin ，但是没说用哪个。从Seura...
细胞通讯-5 实操完整代码
所在位置 D:\细胞通讯-0808 复现视频1.R 1、单样本 2、多组别所在位置 D:\细胞通讯-08...
【单细胞高级绘图】08.细胞通讯_两组比较_气泡图
细胞通讯分析，我很早之前就写过，当时详细介绍了CellPhoneDB的原理、实操、可视化。单细胞分析实录(18):...
Supercell代码实操（二）：COVID-19数据的整合分析
前言在上一期推文：Supercell代码实操（一）：5种肿瘤细胞系的整合分析[https://mp.weixin...
【王兴国营养特训班】第5期第11课7组虞文施
细胞通讯和分化、干细胞、细胞衰老和死亡、癌细胞细胞通讯细胞通讯的方式：直接通讯：间隙连接、膜表面分子接触通讯、...
109-小days6第三课#裂变增长实验室#
以下是实操截图： 2、被动引流5人
宫氏脑针29期公告
培训时间:第一阶段2018年5月29日报到，5月30日-6月3日上课5天，理论+实操+解剖参加过实操班的学员，...
monocle实操+细胞功能富集分析
monocle实操 1 软件安装 1.1 conda安装——Linux https://anaconda.org/...
单细胞分析：质控实操（五）
1. 学习目标构建质量控制指标并评估数据质量适当的应用过滤器去除低质量的细胞 2. 过滤目标过滤数据以仅包含...
Scissor实操（二）：逻辑回归分析流程解读
前言在上一期推文：Scissor实操（一）：利用Cox回归模型识别特征细胞[https://mp.weixin....