Cellchat Cell Communication

# Install CellChat from GitHub (CRAN version may lag)
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install(c("BiocNeighbors", "ComplexHeatmap"))

install.packages("devtools")
devtools::install_github("jinworks/CellChat")

# Core dependencies
install.packages(c("NMF", "ggplot2", "ggalluvial", "igraph",
                   "dplyr", "patchwork", "circlize", "RColorBrewer"))

decisions:
  - id: D1
    param: cellGroupLabels
    kind: derived
    source: upstream
    ask: "Which annotation column defines the sender and receiver groups?"
    default: "the cell-type labels assigned upstream"

  - id: D2
    param: ligandReceptorDatabase
    kind: required
    source: user
    ask: "Which species database, and should it cover all interactions or only secreted signalling, ECM-receptor, or cell-contact?"
    default: "species-matched, all interaction categories"

  - id: D3
    param: minCellsPerGroup
    kind: required
    source: user
    ask: "How few cells may a group have before its inferred signalling stops being trustworthy?"
    default: 10

  - id: D4
    param: expressionAggregation
    kind: optional
    source: user
    ask: "How should per-group expression be summarized - stringently, or more permissively so weaker signals survive?"
    default: "triMean, the most stringent option"

  - id: D5
    param: populationSizeWeighting
    kind: optional
    source: user
    ask: "Should communication probability account for how many cells each group has?"
    default: "weighted"

  - id: D6
    param: bootstrapIterations
    kind: optional_conditional
    source: user
    ask: "More bootstrap rounds give finer p-values at proportionally more runtime - is the default enough?"
    default: 100

  - id: D7
    param: communicationPatterns
    kind: optional_conditional
    source: data
    ask: "How many latent signalling patterns should the data be decomposed into?"
    default: "chosen from the selectK elbow"

  - id: D8
    param: centralityThreshold
    kind: optional
    source: user
    ask: "Which interactions are strong enough to enter the network-centrality analysis?"
    default: "p-value 0.05"

library(CellChat)
library(Seurat)

# Assume `seurat_obj` is a processed Seurat object with cell type identities in Idents()
data.input <- GetAssayData(seurat_obj, assay = "RNA", slot = "data")  # normalized counts
meta       <- data.frame(labels = Idents(seurat_obj), row.names = names(Idents(seurat_obj)))

cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels")
cellchat@DB <- CellChatDB.human  # or CellChatDB.mouse

cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
cellchat <- computeCommunProb(cellchat, type = "triMean")
cellchat <- filterCommunication(cellchat, min.cells = 10)
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)

# Quick summary
print(cellchat)
# e.g. "An object of class CellChat created from a single dataset
#  with 8 cell groups and 312 inferred ligand-receptor pairs"

library(CellChat)
library(Seurat)

# --- Option A: from a Seurat object ---
# seurat_obj must have cell type identities set with Idents() or in meta.data
data.input <- GetAssayData(seurat_obj, assay = "RNA", slot = "data")  # log-normalized
meta <- data.frame(
  labels = Idents(seurat_obj),
  row.names = colnames(seurat_obj)
)
cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels")

# --- Option B: from a count matrix directly ---
# data.input: genes-by-cells normalized matrix (dgCMatrix or dense matrix)
# identity: named factor of cell group labels (length = ncol(data.input))
cellchat <- createCellChat(object = data.input, meta = data.frame(labels = identity),
                           group.by = "labels")

cat("Cell groups:", levels(cellchat@idents), "\n")
cat("Number of cells:", ncol(data.input), "\n")
# Cell groups: B_cell Endothelial Fibroblast Macrophage NK T_cell Tumor
# Number of cells: 12847

# Load database for the appropriate species
CellChatDB <- CellChatDB.human   # use CellChatDB.mouse for mouse data

# Inspect available signaling categories
unique(CellChatDB$interaction$annotation)
# [1] "Secreted Signaling"    "ECM-Receptor"          "Cell-Cell Contact"

# Option 1: Use all interactions (recommended for discovery)
cellchat@DB <- CellChatDB

# Option 2: Subset to secreted ligand-receptor pairs only (reduces noise)
CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling",
                           key = "annotation")
cellchat@DB <- CellChatDB.use

# Subset the CellChat data slots to only genes in the database
cellchat <- subsetData(cellchat)
cat("Genes retained after database subset:", nrow(cellchat@data.signaling), "\n")
# Genes retained after database subset: 1842

# Identify over-expressed genes per cell group (uses Seurat-style wilcoxon test)
cellchat <- identifyOverExpressedGenes(cellchat)

# Map over-expressed genes to ligand-receptor pairs in CellChatDB
cellchat <- identifyOverExpressedInteractions(cellchat)

# Inspect how many interactions were identified per group pair
df.net <- subsetCommunication(cellchat)
cat("Total inferred interactions:", nrow(df.net), "\n")
head(df.net[, c("source", "target", "ligand", "receptor", "prob")], 5)
#      source   target  ligand receptor      prob
# 1   B_cell Macrophage  CD22     PTPRC 0.0318
# 2 Fibroblast    Tumor   FN1     CD44  0.1072
# ...

# Compute pairwise communication probability
# type = "triMean": uses 25th percentile × mean × 25th percentile for robustness
# type = "truncatedMean": uses trimmed mean with threshold parameter trim
cellchat <- computeCommunProb(
  cellchat,
  type          = "triMean",   # recommended default
  trim          = 0.1,         # fraction to trim (only used if type="truncatedMean")
  nboot         = 100,         # bootstrap iterations for p-value estimation
  seed.use      = 42,
  population.size = TRUE       # weight by population size (recommended)
)

# Filter out interactions with too few cells in sender or receiver groups
cellchat <- filterCommunication(cellchat, min.cells = 10)

# Summary of retained interactions
df.net <- subsetCommunication(cellchat)
cat("Interactions after filtering:", nrow(df.net), "\n")
cat("Significant interactions (p<0.05):", sum(df.net$pval < 0.05), "\n")
# Interactions after filtering: 247
# Significant interactions (p<0.05): 189

# Aggregate to pathway level
cellchat <- computeCommunProbPathway(cellchat)

# Build aggregate interaction count and weight networks
cellchat <- aggregateNet(cellchat)

# View significant pathways
cat("Significant signaling pathways:\n")
print(cellchat@netP$pathways)
# [1] "MHC-II"    "COLLAGEN"  "FN1"       "VEGF"      "CXCL"
# [6] "CCL"       "MIF"       "APP"       "GALECTIN"  ...

# Extract pathway-level communication probabilities between groups
df.pathways <- subsetCommunication(cellchat, slot.name = "netP")
head(df.pathways[, c("source", "target", "pathway_name", "prob")], 5)
#        source    target pathway_name     prob
# 1  Fibroblast     Tumor     COLLAGEN   0.2341
# 2  Macrophage  Fibroblast     MIF    0.1876
# ...

# Compute centrality measures for all pathways
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")

# Visualize centrality scores as a heatmap (rows=pathways, cols=cell groups)
# Each dot size: outgoing signal strength; color: incoming signal strength
netAnalysis_signalingRole_heatmap(
  cellchat,
  pattern    = "all",     # "outgoing", "incoming", or "all"
  signaling  = NULL,      # NULL = all pathways; or specify e.g. c("COLLAGEN","VEGF")
  height     = 10,
  color.heatmap = "OrRd"
)

# Identify dominant communication patterns using NMF
# outgoing patterns reveal which cell groups co-activate similar pathways
library(NMF)
selectK(cellchat, pattern = "outgoing")    # elbow plot to choose K
cellchat <- identifyCommunicationPatterns(
  cellchat,
  pattern = "outgoing",
  k       = 3,            # number of latent patterns; choose from selectK elbow
  width   = 8,
  height  = 6
)

library(ggplot2)
library(patchwork)

# --- 7a. Chord diagram: aggregate interaction count and weight ---
par(mfrow = c(1, 2))
netVisual_circle(
  cellchat@net$count,
  vertex.weight = as.numeric(table(cellchat@idents)),
  weight.scale  = TRUE,
  label.edge    = FALSE,
  title.name    = "Number of interactions"
)
netVisual_circle(
  cellchat@net$weight,
  vertex.weight = as.numeric(table(cellchat@idents)),
  weight.scale  = TRUE,
  label.edge    = FALSE,
  title.name    = "Interaction strength"
)

# --- 7b. Heatmap: cell-group × cell-group interaction matrix ---
p1 <- netVisual_heatmap(cellchat, measure = "count",  color.heatmap = "Blues")
p2 <- netVisual_heatmap(cellchat, measure = "weight", color.heatmap = "Reds")
p1 + p2

# --- 7c. Chord diagram for a specific pathway ---
netVisual_aggregate(
  cellchat,
  signaling      = "COLLAGEN",
  layout         = "chord",
  vertex.receiver = NULL   # NULL = show all groups as receivers
)

# --- 7d. Bubble plot: all significant interactions for chosen pathways ---
netVisual_bubble(
  cellchat,
  sources.use = NULL,   # NULL = all senders
  targets.use = NULL,   # NULL = all receivers
  signaling   = c("COLLAGEN", "MIF", "VEGF"),
  remove.isolate = FALSE
)
ggsave("bubble_plot_selected_pathways.pdf", width = 10, height = 8)

# Assume cellchat_ctrl and cellchat_disease are pre-computed CellChat objects
object.list <- list(Control = cellchat_ctrl, Disease = cellchat_disease)
cellchat_merged <- mergeCellChat(object.list, add.names = names(object.list))

# --- Compare total interaction count and strength ---
compareInteractions(cellchat_merged, show.legend = FALSE,
                    group = c(1, 2), measure = "count")
compareInteractions(cellchat_merged, show.legend = FALSE,
                    group = c(1, 2), measure = "weight")

# --- Differential interaction chord diagram (gained/lost connections) ---
netVisual_diffInteraction(cellchat_merged, weight.scale = TRUE)

# --- Identify signaling pathways specific to each condition ---
rankNet(cellchat_merged, mode = "comparison", stacked = TRUE, do.stat = TRUE)

# --- Scatter plot: pathways shifted in information flow ---
rankNetPairwise(
  cellchat_merged,
  comparison = c(1, 2),
  slot.name  = "netP",
  measure    = "prob"
)

P(A → B | L-R) = hill(expr_L_A) × hill(expr_R1_B) × hill(expr_R2_B) × ...

# Inspect the database structure
dim(CellChatDB.human$interaction)   # [1] 2293   17
head(CellChatDB.human$interaction[, c("interaction_name", "pathway_name",
                                       "ligand", "receptor", "annotation")], 4)
#   interaction_name pathway_name ligand receptor          annotation
# 1         TGFB1_TGFBR1_TGFBR2         TGFb  TGFB1  TGFBR1_TGFBR2  Secreted Signaling
# 2          WNT5A_FZD1_LRP5          WNT   WNT5A     FZD1_LRP5   Secreted Signaling
# 3             FN1_CD44             FN1    FN1      CD44        ECM-Receptor
# 4          NOTCH1_DLL4        NOTCH  NOTCH1       DLL4   Cell-Cell Contact