细胞通讯-5 实操完整代码
2022-08-25 本文已影响0人
oceanandshore
所在位置 D:\细胞通讯-0808 复现视频1.R
1、单样本
suppressMessages(if(!require(CellChat))devtools::install_github("sqjin/CellChat"))
#需要安装rtools才能运行
#suppressMessages(reticulate::py_install(packages = 'umap-learn'))
###通过minicnda安装帮助后续评估通路间得相似性
suppressMessages(if(!require(ggplot2))install.packages('ggplot2'))
suppressMessages(if(!require(patchwork))install.packages('patchwork') )
suppressMessages(if(!require(ggalluvial))install.packages('ggalluvial'))
suppressMessages(if(!require(igraph))install.packages('igraph'))
suppressMessages(if(!require(dplyr))install.packages('dplyr'))
suppressMessages(options(stringsAsFactors = FALSE))
# library packages
library(Seurat)
library(dplyr)
library(SeuratData)
library(patchwork) #最强大的拼图包
library(ggplot2)
library(CellChat)
library(ggalluvial)
library(svglite)
rm(list=ls()) #清空所有变量
options(stringsAsFactors = F) #输入数据不自动转换成因子（防止数据格式错误）
getwd()
#load(url("https://ndownloader.figshare.com/files/25950872"))#读取示例数据集
#View(data_humanSkin)
#saveRDS(data_humanSkin,'data_humanSkin.rds')
data_humanSkin <- readRDS('data_humanSkin.rds')
class(data_humanSkin)
## [1] "list"
#示例数据：来源于人类皮肤
#作者发表文献 #https://www.nature.com/articles/s41467-021-21246-9.pdf
#把表达矩阵取出来，用于输入
data.input <- data_humanSkin$data # 标准化过的矩阵
#单细胞对象的话，是取这个pbmc[["RNA"]]@data
# 把临床信息，也就是注释信息取出来
meta <- data_humanSkin$meta # a dataframe with rownames containing cell mata data
unique(meta$condition)
## [1] "LS" "NL"
#原文的作者只取了 LS 进行下游分析
cell.use <- rownames(meta)[meta$condition == "LS"] # 按指定的变量提取细胞
data.input <- data.input[, cell.use]#取出对应细胞
meta = meta[cell.use, ]#取出对应细胞的meta信息
# meta = data.frame(labels = meta$labels[cell.use], row.names = colnames(data.input)) # manually create a dataframe consisting of the cell labels
unique(meta$labels)
## [1] Inflam. FIB FBN1+ FIB APOE+ FIB COL11A1+ FIB cDC2
## [6] LC Inflam. DC cDC1 CD40LG+ TC Inflam. TC
## [11] TC NKT
## 12 Levels: APOE+ FIB FBN1+ FIB COL11A1+ FIB Inflam. FIB cDC1 cDC2 ... NKT
#创建CellChat对象
cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels")
## Create a CellChat object from a data matrix
## Set cell identities for the new CellChat object
## The cell groups used for CellChat analysis are APOE+ FIB FBN1+ FIB COL11A1+ FIB Inflam. FIB cDC1 cDC2 LC Inflam. DC TC Inflam. TC CD40LG+ TC NKT
#创建celllchat对象，group.by指定通讯间的对象，用meta中的注释作为分组依据
######如果上述没有指定meta=meta，可以如下手动加入
#cellchat <- createCellChat(object = data.input, group.by = "labels")
#cellchat <- addMeta(cellchat, meta = meta)
cellchat <- setIdent(cellchat, ident.use = "labels") # set "labels" as default cell identity
#Idents(pbmc) <- 'group'
levels(cellchat@idents) # show factor levels of the cell labels
## [1] "APOE+ FIB" "FBN1+ FIB" "COL11A1+ FIB" "Inflam. FIB" "cDC1"
## [6] "cDC2" "LC" "Inflam. DC" "TC" "Inflam. TC"
## [11] "CD40LG+ TC" "NKT"
groupSize <- as.numeric(table(cellchat@idents)) # number of cells in each cell group
groupSize
## [1] 1228 813 181 484 121 294 67 81 765 266 630 81
###载入数据库并开始计算
CellChatDB <- CellChatDB.human # use CellChatDB.mouse if running on mouse data
showDatabaseCategory(CellChatDB)
#展示细胞组成的比例:
dplyr::glimpse(CellChatDB$interaction)#展示互作的记录
## Rows: 1,939
## Columns: 11
## $ interaction_name <chr> "TGFB1_TGFBR1_TGFBR2", "TGFB2_TGFBR1_TGFBR2", "TGFB…
## $ pathway_name <chr> "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TG…
## $ ligand <chr> "TGFB1", "TGFB2", "TGFB3", "TGFB1", "TGFB1", "TGFB2…
## $ receptor <chr> "TGFbR1_R2", "TGFbR1_R2", "TGFbR1_R2", "ACVR1B_TGFb…
## $ agonist <chr> "TGFb agonist", "TGFb agonist", "TGFb agonist", "TG…
## $ antagonist <chr> "TGFb antagonist", "TGFb antagonist", "TGFb antagon…
## $ co_A_receptor <chr> "", "", "", "", "", "", "", "", "", "", "", "", "",…
## $ co_I_receptor <chr> "TGFb inhibition receptor", "TGFb inhibition recept…
## $ evidence <chr> "KEGG: hsa04350", "KEGG: hsa04350", "KEGG: hsa04350…
## $ annotation <chr> "Secreted Signaling", "Secreted Signaling", "Secret…
## $ interaction_name_2 <chr> "TGFB1 - (TGFBR1+TGFBR2)", "TGFB2 - (TGFBR1+TGFBR2)…
#CellChatDB.use <- CellChatDB
#这一步类似取子集的操作，单独去除“分泌”这一条通路，类似取子集 pbmc < - subset(pbmc, ident ='T cell')
CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling")#取出相应分类用作分析数据库
#CellChatDB.use <- CellChatDB # simply use the default CellChatDB
cellchat@DB <- CellChatDB.use#将数据库内容载入cellchat对象中
#表达量预处理
#这一步是取出表达数据，如果你有感兴趣的基因可以填在features里，没有的话，就写NULL
cellchat <- subsetData(cellchat,features = NULL)#取出表达数据
cellchat <- identifyOverExpressedGenes(cellchat)#寻找高表达的基因#
#类似寻找高变基因FindVariableFeatures()
cellchat <- identifyOverExpressedInteractions(cellchat)#寻找高表达的通路
cellchat <- projectData(cellchat, PPI.human)#投影到PPI
cellchat <- computeCommunProb(cellchat, raw.use = T)#默认计算方式为type = "truncatedMean",
#默认cutoff的值为20%，即表达比例在25%以下的基因会被认为是0， trim = 0.1可以调整比例阈值
cellchat <- filterCommunication(cellchat, min.cells = 10)
#去掉通讯数量很少的细胞
df.net <- subsetCommunication(cellchat)#将细胞通讯预测结果以数据框的形式取出
View(df.net)
DT::datatable(df.net)
write.csv(df.net,'01.df.net.csv')
#这种方式只取通路，数据结构更简单
df.netp <- subsetCommunication(cellchat,slot.name = "netP")
#这一步类似取子集，source可以以细胞类型的名称定义，也可以按照细胞名称中的顺序以数值向量直接取
df.net <- subsetCommunication(cellchat, sources.use = c(1,2), targets.use = c(4,5))
# levels(meta$labels)
#[1] "APOE+ FIB" "FBN1+ FIB" "COL11A1+ FIB" "Inflam. FIB" "cDC1" "cDC2"
#[7] "LC" "Inflam. DC" "TC" "Inflam. TC" "CD40LG+ TC" "NKT"
# c(1,2) 就是"APOE+ FIB" "FBN1+ FIB"这两种细胞，这里也可以换成c("APOE+ FIB","FBN1+ FIB")
#指定输入与输出的细胞集群
#df.net <- subsetCommunication(cellchat, signaling = c("WNT", "TGFb"))#指定通路提取
cellchat <- computeCommunProbPathway(cellchat)
#每对配受体的预测结果存在net中，每条通路的预测结果存在netp中
cellchat <- aggregateNet(cellchat)
#计算联路数与通讯概率，可用sources.use and targets.use指定来源与去向
##################可视化######################
groupSize
## [1] 1228 813 181 484 121 294 67 81 765 266 630 81
library(patchwork)
par(mfrow = c(1,3), xpd=TRUE) #c(1,3)的意思是，画一个 一行三列的图片
netVisual_circle(cellchat@net$count, vertex.weight = groupSize, weight.scale = T,
label.edge= F, title.name = "Number of interactions") # 展示互作数量
netVisual_circle(cellchat@net$weight, vertex.weight = groupSize, weight.scale = T,
label.edge= F, title.name = "Interaction weights/strength") #展示互作权重
netVisual_circle(cellchat@net$weight, vertex.weight = groupSize, weight.scale = T,
label.edge= F, title.name = "Interaction weights/strength",
targets.use = 'cDC2') #只展示与CD2相关的互作，箭头指向DC2
#互作数量与重要性图
mat <- cellchat@net$weight
par(mfrow = c(3,4), xpd=TRUE)
for (i in 1:nrow(mat)) {#
mat2 <- matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))
mat2[i, ] <- mat[i, ]
netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T,
edge.weight.max = max(mat), title.name = rownames(mat)[i])
} #循环绘出每种与其他细胞之间互作的关系
#for (i in 1:nrow(mat)) {
mat2 <- matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))
mat2[i, ] <- mat[i, ]
netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T,
edge.weight.max = max(mat), title.name = rownames(mat)[i])
#####进阶可视化
pathways.show <- df.net$pathway_name#计算到的所有通路
# Hierarchy plot
par(mfrow = c(1,2), xpd=TRUE)
vertex.receiver = seq(1,4)
netVisual_aggregate(cellchat, signaling =pathways.show[1],
vertex.receiver = vertex.receiver,layout = 'hierarchy')
#左侧显示只当细胞的的通讯，右侧显示剩余细胞的通讯情况
# Circle plot
netVisual_aggregate(cellchat, signaling = pathways.show[1], layout = "circle")
#以circle的方式展示通讯
pdf('Rplot05.pdf',width = 7, height = 7)
netVisual_aggregate(cellchat, signaling = pathways.show[1], layout = "chord")
dev.off()
#把所有的通路都画出来
dir.create('chord/')
for(i in 1:length(pathways.show)){
pdf(paste0('chord/',pathways.show[i],'.pdf'),width = 7,height = 7)
netVisual_aggregate(cellchat, signaling = pathways.show[i], layout = "chord")
dev.off()
}
#> length(pathways.show)
#[1] 4
#视频里 length 长度271条通路，为什么我只有4条？？？
# 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')
cowplot::plot_grid(p1, p2 ,align = "h",ncol=2)
netVisual_individual(cellchat, signaling = pathways.show[1],
pairLR.use = LR.show, layout = "circle")
# 取的是LR.show，也就是最显著的一个配受体，所以这个函数可以指定配受体
#LR.show
#[1] "FGF7_FGFR1"
## [[1]]
# Circle plot
netVisual_individual(cellchat, signaling = pathways.show,
pairLR.use = LR.show, layout = "chord")
## [[1]]
# Chord diagram
#####通路展示
pathways.show.all <- cellchat@netP$pathways
# check the order of cell identity to set suitable vertex.receiver
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"
vertex.receiver = seq(1,4)
dir.create('04.pathwat.show')
for (i in 1:length(pathways.show.all)) {
# Visualize communication network associated with both signaling pathway and individual L-R pairs
# netVisual(cellchat, signaling = pathways.show.all[i],
# vertex.receiver = vertex.receiver, layout = "hierarchy")#直接出pdf与svg
# Compute and visualize the contribution of each ligand-receptor pair to the overall signaling pathway
gg <- netAnalysis_contribution(cellchat, signaling = pathways.show.all[i])
ggsave(filename=paste0('04.pathwat.show/',pathways.show.all[i], "_L-R_contribution.pdf"),
plot=gg, width = 3, height = 2, units = 'in', dpi = 300)
}
#循环存出连线图
gg
#####多个细胞通讯通路气泡图可视化
#指定发送者与接收者，展示所有配受体对
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两张图是一模一样的，一个是数字，一个是用标签标识，结果都一样
#指定信号通路提取配受体
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"
#指定通路，发出、接收的点图
2、多组别
所在位置 D:\细胞通讯-0808 多组别的细胞通讯.R
#这部分的内容主要分析目的为：
#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
#查看细胞通路在两组间的富集程度
gg1 <- rankNet(cellchat, mode = "comparison", stacked = T, do.stat = TRUE)
gg2 <- rankNet(cellchat, mode = "comparison", stacked = F, do.stat = TRUE)
gg1 + gg2
#以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"
#以上是直接展示二组“相减”的结果，当然你也可以直接将两组分开展示:
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]))
}
#同理，用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
#同样是第二个组减去第一个组，左侧是细胞通讯数量，右侧是细胞通讯强度，和上面的圈图一样的。
#展示特定通路， 老三样
#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组有更多的细胞参与这个通路
#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
#弦图
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]))
}
#基因表达情况
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.
#研究细胞类型之间的通讯差异
#按照组别拆分整合
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相较其他两类，变化最大
#两组分别展示
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组更强一些
#展示两组间各类细胞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
#展示两组间的差异
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 在两组中差异较大
#Rplot13
#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
#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
#从配受体的角度画一画气泡图
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
#弦图
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
#基因表达情况
cellchat@meta$datasets = factor(cellchat@meta$datasets,
levels = c("NL", "LS")#设定组别的level
)
plotGeneExpression(cellchat, signaling = "CXCL", split.by = "datasets", colors.ggplot = T)
#这部分结果展示建议用Seurat做
## 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.
#版本信息
#sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ComplexHeatmap_2.10.0 cowplot_1.1.1 patchwork_1.1.1
## [4] CellChat_1.1.3 Biobase_2.54.0 BiocGenerics_0.40.0
## [7] ggplot2_3.3.5 igraph_1.2.8 dplyr_1.0.7
## [10] imager_0.42.13 magrittr_2.0.1
##
## loaded via a namespace (and not attached):
## [1] circlize_0.4.13 systemfonts_1.0.3 NMF_0.23.0
## [4] plyr_1.8.6 lazyeval_0.2.2 splines_4.1.2
## [7] listenv_0.8.0 scattermore_0.7 gridBase_0.4-7
## [10] digest_0.6.29 foreach_1.5.1 htmltools_0.5.2
## [13] magick_2.7.3 tiff_0.1-11 ggalluvial_0.12.3
## [16] fansi_1.0.2 tensor_1.5 cluster_2.1.2
## [19] doParallel_1.0.16 ROCR_1.0-11 sna_2.6
## [22] globals_0.14.0 matrixStats_0.61.0 svglite_2.0.0
## [25] spatstat.sparse_2.0-0 jpeg_0.1-9 colorspace_2.0-2
## [28] ggrepel_0.9.1 xfun_0.29 crayon_1.4.2
## [31] jsonlite_1.7.3 spatstat.data_2.1-0 survival_3.2-13
## [34] zoo_1.8-9 iterators_1.0.13 glue_1.6.0
## [37] polyclip_1.10-0 registry_0.5-1 gtable_0.3.0
## [40] leiden_0.3.9 GetoptLong_1.0.5 future.apply_1.8.1
## [43] shape_1.4.6 abind_1.4-5 scales_1.1.1
## [46] DBI_1.1.1 rngtools_1.5.2 miniUI_0.1.1.1
## [49] Rcpp_1.0.8 viridisLite_0.4.0 xtable_1.8-4
## [52] clue_0.3-60 spatstat.core_2.3-1 reticulate_1.22
## [55] stats4_4.1.2 htmlwidgets_1.5.4 httr_1.4.2
## [58] FNN_1.1.3 RColorBrewer_1.1-2 ellipsis_0.3.2
## [61] Seurat_4.0.5 ica_1.0-2 pkgconfig_2.0.3
## [64] farver_2.1.0 uwot_0.1.10 deldir_1.0-6
## [67] sass_0.4.0 utf8_1.2.2 tidyselect_1.1.1
## [70] labeling_0.4.2 rlang_0.4.12 reshape2_1.4.4
## [73] later_1.3.0 munsell_0.5.0 tools_4.1.2
## [76] generics_0.1.1 statnet.common_4.5.0 ggridges_0.5.3
## [79] evaluate_0.14 stringr_1.4.0 fastmap_1.1.0
## [82] goftest_1.2-3 yaml_2.2.1 knitr_1.37
## [85] fitdistrplus_1.1-6 purrr_0.3.4 RANN_2.6.1
## [88] readbitmap_0.1.5 nlme_3.1-155 pbapply_1.5-0
## [91] future_1.23.0 mime_0.12 compiler_4.1.2
## [94] plotly_4.10.0 png_0.1-7 spatstat.utils_2.2-0
## [97] tibble_3.1.6 bslib_0.3.1 stringi_1.7.6
## [100] highr_0.9 RSpectra_0.16-0 lattice_0.20-45
## [103] Matrix_1.4-0 vctrs_0.3.8 pillar_1.6.4
## [106] lifecycle_1.0.1 spatstat.geom_2.3-0 lmtest_0.9-39
## [109] jquerylib_0.1.4 GlobalOptions_0.1.2 RcppAnnoy_0.0.19
## [112] data.table_1.14.2 irlba_2.3.3 httpuv_1.6.3
## [115] R6_2.5.1 promises_1.2.0.1 network_1.17.1
## [118] gridExtra_2.3 KernSmooth_2.23-20 bmp_0.3
## [121] IRanges_2.28.0 parallelly_1.30.0 codetools_0.2-18
## [124] MASS_7.3-54 assertthat_0.2.1 pkgmaker_0.32.2
## [127] rjson_0.2.21 withr_2.4.3 SeuratObject_4.0.3
## [130] sctransform_0.3.2 S4Vectors_0.32.2 mgcv_1.8-38
## [133] parallel_4.1.2 rpart_4.1-15 tidyr_1.1.4
## [136] coda_0.19-4 rmarkdown_2.11 Rtsne_0.15
## [139] shiny_1.7.1
