CellChat 细胞通讯分析教程：从 Seurat 对象到通讯网络、通路与多组比较

# 方案 A：conda 装依赖（本页实测，macOS arm64，约 4 分钟 + 编译 2.5 分钟）
micromamba create -y -n cellchat -c conda-forge -c bioconda \
  r-base=4.5 r-seurat r-remotes r-rcppeigen c-compiler cxx-compiler fortran-compiler make \
  r-future r-future.apply r-pbapply r-irlba r-ggalluvial r-svglite r-ggrepel r-circlize \
  bioconductor-complexheatmap bioconductor-biocneighbors bioconductor-biobase \
  r-rspectra r-reticulate r-sna r-fnn r-shape r-ggpubr r-ggnetwork r-plotly r-shiny r-bslib \
  r-collapse r-pkgmaker r-registry r-rngtools r-gridbase r-doparallel r-cowplot r-uwot
micromamba run -n cellchat Rscript -e '
options(repos = c(CRAN = "https://cloud.r-project.org"))
install.packages("NMF")                                        # conda-forge 没有 r-nmf
remotes::install_github("immunogenomics/presto", upgrade = "never")
remotes::install_github("jinworks/CellChat", upgrade = "never")
packageVersion("CellChat")'

# 方案 B：已有 R 4.4 以上
# install.packages(c("BiocManager", "remotes", "NMF"))
# BiocManager::install(c("ComplexHeatmap", "BiocNeighbors", "Biobase"))
# remotes::install_github(c("immunogenomics/presto", "jinworks/CellChat"), upgrade = "never")

library(Seurat)
library(CellChat)

seu <- readRDS("pbmc3k_annot.rds")          # 已注释的 Seurat v5 对象
seu <- JoinLayers(seu)                       # 多样本对象先合并 layer，否则只取到第一层
seu$celltype <- droplevels(factor(seu$celltype))
stopifnot(!anyNA(seu$celltype))              # 标签不能有 NA

data.input <- LayerData(seu, assay = "RNA", layer = "data")   # log 归一化值，不是 counts
stopifnot(min(data.input) >= 0)              # 不能用 scale.data 或 integrated
meta <- data.frame(
  labels  = seu$celltype,
  samples = factor(seu$orig.ident),          # v2 需要 samples 列
  row.names = colnames(seu)
)
# 数字聚类号要加前缀，CellChat 不接受标签 "0"
# meta$labels <- factor(paste0("C", seu$seurat_clusters))

cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels")
table(cellchat@idents)

CellChatDB <- CellChatDB.human               # 小鼠用 CellChatDB.mouse
cellchat@DB <- subsetDB(CellChatDB)          # 去掉 Non-protein Signaling 的全部蛋白类配受体
# cellchat@DB <- subsetDB(CellChatDB, search = "Secreted Signaling", key = "annotation")

future::plan("sequential")                   # 小数据集顺序执行更快
cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat)          # 默认用 presto 做 Wilcoxon
cellchat <- identifyOverExpressedInteractions(cellchat)

cellchat <- computeCommunProb(cellchat,
                              type = "truncatedMean", trim = 0.1,  # trim 只在 truncatedMean 下生效
                              population.size = FALSE)              # 未分选样本可设 TRUE
cellchat <- filterCommunication(cellchat, min.cells = 10)
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)

df.net  <- subsetCommunication(cellchat)                      # 配体-受体层面
df.netP <- subsetCommunication(cellchat, slot.name = "netP")  # 通路层面
nrow(df.net); cellchat@netP$pathways
saveRDS(cellchat, "cellchat_pbmc3k.rds")

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

pathways.show <- "MIF"                       # 先确认它在 cellchat@netP$pathways 里
netVisual_aggregate(cellchat, signaling = pathways.show, layout = "chord")
netVisual_heatmap(cellchat, signaling = pathways.show, color.heatmap = "Reds")
netVisual_bubble(cellchat, sources.use = c("CD14 Mono", "FCGR3A Mono"),
                 targets.use = c("CD8 T", "NK", "B"), remove.isolate = FALSE)

cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")
netAnalysis_signalingRole_scatter(cellchat)
netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing") +
  netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming")

library(NMF); library(ggalluvial)
selectK(cellchat, pattern = "outgoing", k.range = 2:6, nrun = 10)   # 默认 2:10、nrun = 30 很慢
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "outgoing", k = 3)
netAnalysis_river(cellchat, pattern = "outgoing")

run_cellchat <- function(obj) {
  meta <- data.frame(labels = droplevels(obj$celltype),  # 本组没有的类型必须去掉
                      samples = factor(obj$sample), row.names = colnames(obj))
  cc <- createCellChat(LayerData(obj, assay = "RNA", layer = "data"), meta = meta, group.by = "labels")
  cc@DB <- subsetDB(CellChatDB.human, search = "Secreted Signaling", key = "annotation")
  cc <- subsetData(cc)
  cc <- identifyOverExpressedGenes(cc)
  cc <- identifyOverExpressedInteractions(cc)
  cc <- computeCommunProb(cc, type = "truncatedMean", trim = 0.1)  # 两组参数必须相同
  cc <- filterCommunication(cc, min.cells = 10)
  cc <- computeCommunProbPathway(cc)
  cc <- aggregateNet(cc)
  netAnalysis_computeCentrality(cc, slot.name = "netP")
}
# 两组使用同一套标签（同一列注释、同样的写法）
seu$celltype <- factor(seu$celltype)
cc.ctrl <- run_cellchat(subset(seu, stim == "CTRL"))
cc.stim <- run_cellchat(subset(seu, stim == "STIM"))

# 某组缺少某个细胞类型时，先补齐到同一组标签
grp <- union(levels(cc.ctrl@idents), levels(cc.stim@idents))
if (!identical(levels(cc.ctrl@idents), levels(cc.stim@idents))) {
  cc.ctrl <- liftCellChat(cc.ctrl, group.new = grp)
  cc.stim <- liftCellChat(cc.stim, group.new = grp)
}
object.list <- list(CTRL = cc.ctrl, STIM = cc.stim)
cellchat <- mergeCellChat(object.list, add.names = names(object.list))

compareInteractions(cellchat, show.legend = FALSE, group = c(1, 2)) +
  compareInteractions(cellchat, show.legend = FALSE, group = c(1, 2), measure = "weight")
netVisual_diffInteraction(cellchat, weight.scale = TRUE, measure = "weight")  # 红色 = STIM 增强
rankNet(cellchat, mode = "comparison", measure = "weight", stacked = TRUE, do.stat = TRUE)
netVisual_bubble(cellchat, sources.use = "CD14 Mono", comparison = c(1, 2),
                 max.dataset = 2, title.name = "Increased in STIM", angle.x = 45, remove.isolate = TRUE)