空转 | CellChat-V2，揭秘空间转录组数据的细胞通讯分析

options(stringsAsFactors = FALSE)library(CellChat)library(Seurat)library(tidyverse)library(viridis)library(RColorBrewer)#载入数据load("Brain_ST_scRNA.sBio.Rdata")#查看数据情况Brain_ST@meta.data$celltype <- paste("C",Brain_ST$seurat_clusters,sep = "")Idents(Brain_ST) <- "celltype"head(Brain_ST)#可定义颜色color.use <- scPalette(nlevels(Brain_ST))names(color.use) <- levels(Brain_ST)SpatialDimPlot(Brain_ST, label = TRUE, label.size = 3, cols = color.use)

#矩阵信息data.input = Seurat::GetAssayData(Brain_ST, slot = "data", assay = "SCT") #meta信息meta = data.frame(labels = Idents(Brain_ST), #名字自定义                  row.names = names(Idents(Brain_ST))) # manually create a dataframe consisting of the cell labelsunique(meta$labels)
# 空间图像信息spatial.locs = Seurat::GetTissueCoordinates(Brain_ST, scale = NULL,                                             cols = c("imagerow", "imagecol")) # Scale factors and spot diameters 信息 scale.factors = jsonlite::fromJSON(txt =                                    file.path("E:/bioinformation/sc_ST/data/Brain_35707680/spatial", 'scalefactors_json.json'))scale.factors = list(spot.diameter = 65, spot = scale.factors$spot_diameter_fullres, # these two information are required                     fiducial = scale.factors$fiducial_diameter_fullres, hires = scale.factors$tissue_hires_scalef, lowres = scale.factors$tissue_lowres_scalef # these three information are not required)

cellchat <- createCellChat(object = data.input,                            meta = meta,                            group.by = "labels", #前面的meta ，定义的名字是labels                           datatype = "spatial", ###                           coordinates = spatial.locs,                            scale.factors = scale.factors)

CellChatDB <- CellChatDB.human # use CellChatDB.mouse if running on mouse data
# use a subset of CellChatDB for cell-cell communication analysis#CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling", key = "annotation") # use Secreted Signaling# use all CellChatDB for cell-cell communication analysisCellChatDB.use <- CellChatDB # simply use the default CellChatDB
# set the used database in the objectcellchat@DB <- CellChatDB.use

# subset the expression data of signaling genes for saving computation costcellchat <- subsetData(cellchat) # This step is necessary even if using the whole databasefuture::plan("multisession", workers = 1) #笔记本可以选1##识别过表达基因cellchat <- identifyOverExpressedGenes(cellchat)#识别过表达配体受体对cellchat <- identifyOverExpressedInteractions(cellchat)

cellchat <- computeCommunProb(cellchat,                               type = "truncatedMean", trim = 0.1,                               distance.use = TRUE,                               scale.distance = 0.01)cellchat <- filterCommunication(cellchat, min.cells = 10)

#计算每个信号通路相关的所有配体-受体相互作用的通信结果cellchat <- computeCommunProbPathway(cellchat)#计算整合的细胞类型之间通信结果cellchat <- aggregateNet(cellchat)

groupSize <- as.numeric(table(cellchat@idents))par(mfrow = c(1,2), xpd=TRUE)netVisual_circle(cellchat@net$count, vertex.weight = rowSums(cellchat@net$count),                  weight.scale = T, label.edge= F, title.name = "Number of interactions")netVisual_circle(cellchat@net$weight, vertex.weight = rowSums(cellchat@net$weight),                  weight.scale = T, label.edge= F, title.name = "Interaction weights/strength")

p1 <- netVisual_heatmap(cellchat, measure = "count", color.heatmap = "Blues")
p2 <- netVisual_heatmap(cellchat, measure = "weight", color.heatmap = "Blues")
p1 + p2

cellchat@netP$pathwayspathways.show <- c("SPP1")levels(cellchat@idents)   #[1] "C9" "C3" "C5" "C4" "C0" "C7" "C6" "C8" "C1" "C2"     vertex.receiver = c(1,2,4,5)  #选择的是levels(cellchat@idents) 中的netVisual_aggregate(cellchat, signaling = pathways.show,                      vertex.receiver = vertex.receiver,layout = "hierarchy")

# Circle plotpar(mfrow=c(1,1))netVisual_aggregate(cellchat, signaling = pathways.show, layout = "circle")
# Spatial plotpar(mfrow=c(1,1))netVisual_aggregate(cellchat, signaling = pathways.show, layout = "spatial",                     edge.width.max = 2, vertex.size.max = 1,                     alpha.image = 0.2, vertex.label.cex = 3.5)

# Compute the network centrality scorescellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP") # the slot 'netP' means the inferred intercellular communication network of signaling pathways# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groupspar(mfrow=c(1,1))netAnalysis_signalingRole_network(cellchat, signaling = pathways.show,                                   width = 8, height = 2.5, font.size = 10)

#指定受体-配体细胞类型netVisual_bubble(cellchat, sources.use = c(3,5),                  targets.use = c(1,2,4,6), remove.isolate = FALSE)

netVisual_bubble(cellchat, sources.use = c(3,5), targets.use = c(1,2,4,6),                                   signaling = c("TGFb","SPP1"), remove.isolate = FALSE)

# Take an input of a ligand-receptor pair and show expression in binaryspatialFeaturePlot(cellchat, pairLR.use = "IGF1_IGF1R", point.size = 1.5,                    do.binary = TRUE, cutoff = 0.05, enriched.only = F,                    color.heatmap = "Reds", direction = 1)