空间转录组细胞通讯分析：CellChat v2（Nature Protocols，IF: 16.0/Q1）

rm(list=ls())
options(stringsAsFactors = FALSE)
library(Seurat)
library(CellChat)
packageVersion("CellChat")
library(patchwork)

## 加载数据
# Here we load a Seurat object of 10X Visium mouse cortex data and its associated cell meta data
# load("visium_mouse_cortex_annotated_full.RData")
load("visium_mouse_cortex_annotated.RData")
ls()

# show the image and annotated spots
nlevels(visium.brain)
color.use <- scPalette(nlevels(visium.brain)); 
names(color.use) <- levels(visium.brain)
color.use

# 看一下切片注释结果
Seurat::SpatialDimPlot(visium.brain, label = T, label.size = 3, cols = color.use, pt.size.factor = 1600)

# Prepare input data for CelChat analysis
# normalized data matrix
data.input = Seurat::GetAssayData(visium.brain, slot = "data", assay = "SCT") 
data.input[1:5,1:5]
dim(data.input)

# define the meta data: 
# a column named `samples` should be provided for spatial transcriptomics analysis, 
# which is useful for analyzing cell-cell communication by aggregating multiple samples/replicates. 
# Of note, for comparison analysis across different conditions, 
# users still need to create a CellChat object seperately for each condition.  

# manually create a dataframe consisting of the cell labels
meta = data.frame(labels = Seurat::Idents(visium.brain), samples = "sample1", 
                  row.names = names(Seurat::Idents(visium.brain))) 
meta$samples <- factor(meta$samples)
# check the cell labels
unique(meta$labels) 

# check the sample labels
unique(meta$samples) 

# 加载空间信息
# Spatial locations of spots from full (NOT high/low) resolution images are required. For 10X Visium, 
# this information is in `tissue_positions.csv`. 
spatial.locs = Seurat::GetTissueCoordinates(visium.brain, scale = NULL, cols = c("imagerow", "imagecol")) 
head(spatial.locs)

# Spatial factors of spatial coordinates
# For 10X Visium, the conversion factor of converting spatial coordinates from Pixels to Micrometers can be computed as the ratio of the theoretical spot size (i.e., 65um) over the number of pixels that span the diameter of a theoretical spot size in the full-resolution image 

# (i.e., 'spot_diameter_fullres' in pixels in the 'scalefactors_json.json' file). 

# Of note, the 'spot_diameter_fullres' factor is different from the `spot` in Seurat object and thus users still need to get the value from the original json file. 

scalefactors = jsonlite::fromJSON(txt = file.path('../data/V1_Mouse_Brain_Sagittal_Anterior/spatial/scalefactors_json.json'))
scalefactors
scalefactors$spot_diameter_fullres
spot.size = 65 # the theoretical spot size (um) in 10X Visium
conversion.factor = spot.size/scalefactors$spot_diameter_fullres
spatial.factors = data.frame(ratio = conversion.factor, tol = spot.size/2)
spatial.factors

d.spatial <- computeCellDistance(coordinates = spatial.locs, ratio = spatial.factors$ratio, tol = spatial.factors$tol)
# this value should approximately equal 100um for 10X Visium data
min(d.spatial[d.spatial!=0]) 

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

cellchat

## 配体受体库
# use CellChatDB.human if running on human data
CellChatDB <- CellChatDB.mouse 
showDatabaseCategory(CellChatDB)

# Show the structure of the database
dplyr::glimpse(CellChatDB$interaction)

# 使用其中一部分
# use a subset of CellChatDB for cell-cell communication analysis
# use Secreted Signaling
CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling", key = "annotation")

# Only uses the Secreted Signaling from CellChatDB v1
#  CellChatDB.use <- subsetDB(CellChatDB, search = list(c("Secreted Signaling"), c("CellChatDB v1")), key = c("annotation", "version"))

# use all CellChatDB except for "Non-protein Signaling" for cell-cell communication analysis
# CellChatDB.use <- subsetDB(CellChatDB)

# use all CellChatDB for cell-cell communication analysis
# CellChatDB.use <- CellChatDB # simply use the default CellChatDB. We do not suggest to use it in this way because CellChatDB v2 includes "Non-protein Signaling" (i.e., metabolic and synaptic signaling) that can be only estimated from gene expression data. 

# 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, variable.both = F)
cellchat

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

cellchat <- filterCommunication(cellchat, min.cells = 10)

# summarizing the communication probabilities of all ligands-receptors interactions associated with each signaling pathway
cellchat <- computeCommunProbPathway(cellchat)

# summarizing the communication probability.
cellchat <- aggregateNet(cellchat)

# 提取推断结果
df.net <- subsetCommunication(cellchat)
#df.net <- subsetCommunication(cellchat, sources.use = c(1,2), targets.use = c(4,5))
#df.net <- subsetCommunication(cellchat, signaling = c("WNT", "TGFb"))

# All the signaling pathways showing significant communications can be accessed by 
df <- cellchat@netP$pathways

# Spatial plot
par(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)

# USER can show this information on the spatial transcriptomics when visualizing a signaling network, e.g., bigger circle indicates larger incoming signaling
par(mfrow=c(1,1))
netVisual_aggregate(cellchat, signaling = pathways.show, layout = "spatial", edge.width.max = 2, alpha.image = 0.2, vertex.weight = "incoming", vertex.size.max = 4, vertex.label.cex = 3.5)

# Take an input of a few genes
spatialFeaturePlot(cellchat, features = c("Igf1","Igf1r"), point.size = 0.8, color.heatmap = "Reds", direction = 1)

# Take an input of a ligand-receptor pair
spatialFeaturePlot(cellchat, pairLR.use = "IGF1_IGF1R", 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 = "IGF1_IGF1R", point.size = 1, do.binary = TRUE, cutoff = 0.05, enriched.only = F, color.heatmap = "Reds", direction = 1)

saveRDS(cellchat, file = "cellchat_visium_mouse_cortex.rds")