🤩 Cellchat | 空间转录组也可以做细胞通讯啦！~（二）（多个数据集篇）

rm(list = ls())
library(tidyverse)
library(CellChat)
library(patchwork)
library(Seurat)

seu1 <- readRDS("./Cellchat_spatial_multi/Spatial_A1_adult_with_predictions.RDS")

seu2 <- readRDS("./Cellchat_spatial_multi/Spatial_A2_adult_with_predictions.RDS")

# assign label to each spot based on the maximum predicted probabilities
assignLabels <- function(object, prediction = "predictions") {
  pred <- object[[prediction]]@data
  pred <- pred[1:(nrow(pred)-1), ]
  # label each spot based on the maximum prediction probability
  labels = rownames(pred)[apply(pred, 2, which.max)]
  names(labels) <- colnames(pred)
  object$labels <- factor(labels)
  Idents(object) <- "labels"
  return(object)
}
seu1 <- assignLabels(seu1, prediction = "adult.predictions")
seu2 <- assignLabels(seu2, prediction = "adult.predictions")

# show the image and annotated spots
color.use <- scPalette(nlevels(seu1)); names(color.use) <- levels(seu1)
p1 <- Seurat::SpatialDimPlot(seu1, label = F, label.size = 3, cols = color.use)
color.use <- scPalette(nlevels(seu2)); names(color.use) <- levels(seu2)
p2 <- Seurat::SpatialDimPlot(seu2, label = F, label.size = 3, cols = color.use) + NoLegend()
p1 + p2

# Prepare input data for CelChat analysis
data.input1 = Seurat::GetAssayData(seu1, slot = "data", assay = "SCT") # normalized data matrix

data.input2 = Seurat::GetAssayData(seu2, slot = "data", assay = "SCT") 

genes.common <- intersect(rownames(data.input1), rownames(data.input2))
colnames(data.input1) <- paste0("A1_", colnames(data.input1))
colnames(data.input2) <- paste0("A2_", colnames(data.input2))
data.input <- cbind(data.input1[genes.common, ], data.input2[genes.common, ])

meta1 = data.frame(labels = Idents(seu1), samples = "A1") # manually create a dataframe consisting of the cell labels
meta2 = data.frame(labels = Idents(seu2), samples = "A2") 

meta <- rbind(meta1, meta2)
rownames(meta) <- colnames(data.input)

meta$labels <- factor(meta$labels, levels = levels(Idents(seu1)))
meta$samples <- factor(meta$samples, levels = c("A1", "A2"))
unique(meta$labels) # check the cell labels
unique(meta$samples) # check the sample labels

# load spatial transcriptomics information
spatial.locs1 = Seurat::GetTissueCoordinates(seu1, scale = NULL, cols = c("imagerow", "imagecol")) 
spatial.locs2 = Seurat::GetTissueCoordinates(seu2, scale = NULL, cols = c("imagerow", "imagecol")) 
spatial.locs <- rbind(spatial.locs1, spatial.locs2)
rownames(spatial.locs) <- colnames(data.input)

# Scale factors of spatial coordinates
scalefactors1 = jsonlite::fromJSON(txt = file.path("./Cellchat_spatial_multi/", 'A1_scalefactors_json.json'))

spot.size = 65 # the theoretical spot size (um) in 10X Visium
conversion.factor1 = spot.size/scalefactors1$spot_diameter_fullres
spatial.factors1 = data.frame(ratio = conversion.factor1, tol = spot.size/2)

scalefactors2 = jsonlite::fromJSON(txt = file.path("Cellchat_spatial_multi/", 'A2_scalefactors_json.json'))
conversion.factor2 = spot.size/scalefactors2$spot_diameter_fullres
spatial.factors2 = data.frame(ratio = conversion.factor2, tol = spot.size/2)

spatial.factors <- rbind(spatial.factors1, spatial.factors2)
rownames(spatial.factors) <- c("A1", "A2")

cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels",
                           datatype = "spatial", coordinates = spatial.locs, spatial.factors = spatial.factors)

cellchat

CellChatDB <- CellChatDB.human # use CellChatDB.human if running on human data

# use a subset of CellChatDB for cell-cell communication analysis
CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling", key = "annotation") # use Secreted Signaling

# set the used database in the object
cellchat@DB <- CellChatDB.use

# subset the expression data of signaling genes for saving computation cost

cellchat <- subsetData(cellchat) # This step is necessary even if using the whole database
future::plan("multisession", workers = 4) 
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)

cellchat <- computeCommunProb(cellchat, type = "truncatedMean", trim = 0.1, 
                              distance.use = FALSE, interaction.range = 250, scale.distance = NULL,
                              contact.dependent = TRUE, contact.range = 100)

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")

netVisual_heatmap(cellchat, measure = "count", color.heatmap = "Blues")
netVisual_heatmap(cellchat, measure = "weight", color.heatmap = "Blues")

pathways.show <- c("EGF") 
# Circle plot
par(mfrow=c(1,1), xpd=TRUE)
netVisual_aggregate(cellchat, signaling = pathways.show, layout = "circle")

# Spatial plot
par(mfrow=c(1,1))
netVisual_aggregate(cellchat, signaling = pathways.show, sample.use = "A1", layout = "spatial", edge.width.max = 2, vertex.size.max = 1, alpha.image = 0.2, vertex.label.cex = 0)

# Compute the network centrality scores
cellchat <- 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 groups
par(mfrow=c(1,1))
netAnalysis_signalingRole_network(cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10)

par(mfrow=c(1,1))
netVisual_aggregate(cellchat, signaling = pathways.show, sample.use = "A1", layout = "spatial", edge.width.max = 2, alpha.image = 0.2, vertex.weight = "incoming", vertex.size.max = 6, vertex.label.cex = 0)

netAnalysis_contribution(cellchat, signaling = pathways.show)

# Take an input of a few genes
spatialFeaturePlot(cellchat, features = c("AREG","EGFR"), sample.use = "A1", point.size = 0.8, color.heatmap = "Reds", direction = 1)
spatialFeaturePlot(cellchat, features = c("AREG","EGFR"), sample.use = "A2",point.size = 0.8, color.heatmap = "Reds", direction = 1)

spatialFeaturePlot(cellchat, pairLR.use = "AREG_EGFR", sample.use = "A1", point.size = 0.5, do.binary = FALSE, cutoff = 0.05, enriched.only = F, color.heatmap = "Reds", direction = 1)

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