CellChat demo

1 Introduction

For detailed info, see:

• https://github.com/jinworks/CellChat.

• https://htmlpreview.github.io/?https://github.com/jinworks/CellChat/blob/master/tutorial/CellChat-vignette.html.

• https://htmlpreview.github.io/?https://github.com/jinworks/CellChat/blob/master/tutorial/Comparison_analysis_of_multiple_datasets.html.

• https://htmlpreview.github.io/?https://github.com/jinworks/CellChat/blob/master/tutorial/Comparison_analysis_of_multiple_datasets_with_different_cellular_compositions.html.

• https://figshare.com/authors/Suoqin_Jin/14184602.

2 Preliminary works

library(CellChat)
library(patchwork)

Sys.unsetenv("RETICULATE_PYTHON")
reticulate::use_python("/home/yangrui/softwares/micromamba/envs/envs/py3.12/bin/python", required = TRUE)

work_dir <- "/data/users/yangrui/lab_projs/test"
setwd(work_dir)

3 CellChat for a single scRNA-seq dataset

3.1 Infer cellular communication networks

3.1.1 Prepare the input data

Normalized gene by cell matrix and cell group info data frame.

rda_file <- "data/data_humanSkin_CellChat.rda"

load(rda_file)

# normalized data matrix
# with gene names as row names, and cell barcodes as column names
norm_data <- data_humanSkin$data
# a dataframe with cell barcodes as row names containing cell group info in columns
metadata <- data_humanSkin$meta

# subset the input data for CellChat analysis
sample_name <- "LS"
use_cells <- row.names(metadata)[metadata$condition == sample_name]
norm_data <- norm_data[, use_cells]
metadata <- metadata[use_cells, ]
# CellChat needs a column named `samples`
# if not given, CellChat will add this column and
# all cells are assumed to belong to `sample1`
metadata$samples <- factor(metadata$condition)

3.1.2 Create the CellChat object

cellchat <- createCellChat(object = norm_data, meta = metadata, group.by = "labels")

3.1.3 Set the interaction database

# for human samples, use the database `CellChatDB.human`
# for mouse samples, use the database `CellChatDB.mouse`
CellChatDB <- CellChatDB.human
# CellChatDB <- CellChatDB.mouse

showDatabaseCategory(CellChatDB)

# use a subset of CellChatDB for cell-cell communication analysis if needed
# unique(CellChatDB$interaction$annotation)
CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling")

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

3.1.4 Identify over-expressed ligands or receptors

# identify over-expressed ligands or receptors in one cell group
# and then identify over-expressed ligand-receptor interactions if either ligand or receptor is over-expressed
# this step is necessary even if the whole database is used!
# it seems that this will put the expression data containing only those genes listed in CellChatDB into the slot cellchat@data.signaling
cellchat <- subsetData(cellchat)

# do parallel
future::plan("multisession", workers = 20)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)

3.1.5 Infer cell-cell communication network

Note: the number of inferred signaling pathways depends on the method (specified by parameter type) for calculating average gene expression per cell group. Valid methods include:

• triMean: \(\frac{Q_1 + 2 Q_2 + Q_3}{4}\) where \(Q_n\) is the n-th quantile. If the number of cells expressed < 25%, then the TriMean will be 0. In contrast with arithmetic mean, it is more insensitive to outliers, but containing more info than median.

• truncatedMean: to identify weak signaling, e.g., infer L-R pairs in groups with the number of cells expressed < 25%. It needs the parameter trim (\(0-0.25\)) to be set (e.g., trim = 0.1), and then cells with expression values < \(Q_10\) or > \(Q_90\) will be discarded before calculating the arithmetic mean.

To determine a proper value of trim for signaling genes of interest, we can use function computeAveExpr(cellchat, features = c("CXCL12", "CXCR4"), type = "truncatedMean", trim = 0.1) to check their average expression values.

• thresholdedMean: if you confidently know that expression values lower than some threshold are purely technical effects, then we can reset those expression values to be 0 before calculating the arithmetic mean.

• median.

By the way, in mathematics, we have also other ways to calculate means, such as arithmetic mean; geometric mean \(\sqrt[n]{n_1 n_2 \cdots n_n}\), which will be zero if any value is zero; harmonic mean \(\frac{n}{\sum{\frac{1}{x_i}}}\), which is really sensitive to extremely low values making also the harmonic mean extremely low.

# single-cell data with shallow sequencing depth will suffer from severe dropout effects (i.e., a lot of biologically expressed L/R genes will be 0).
# to mitigate dropout effects, CellChat provides the function smoothData() to smooth gene expression based on known protein-protein interactions with highly confident experiment evidence.
# this function should be run before function computeCommunProb().
# `adj` can be either `PPI.human` for human or `PPI.mouse` for mouse
# cellchat <- smoothData(cellchat, adj = PPI.human)
# `raw.use = FALSE` should be set if using cellchat@data.project generated by function smoothData() instead of cellchat@data.signaling
cellchat <- computeCommunProb(cellchat, type = "triMean")

3.1.6 Filter the cell-cell communication

# filter cell-cell communication if there are only few number of cells in certain cell groups (`min.cells`)
# or inconsistent cell-cell communication across samples (`min.samples`)
# for keeping those interactions associated with rare populations when min.samples >= 2, set `rare.keep = TRUE`
# 即保留这些互作：如果存在至少 2 个生物学重复，对于一个 rare population 的互作，其在每一个重复中可能均无法被判定为有效的互作（如不满足 min.cells），但是把这些生物学重复 merge 后满足了保留的条件
# set `nonFilter.keep = TRUE` to keep the non-filtered cell-cell comunication in the CellChat object
# this avoids redundant computeCommunProb() calls when iteratively tuning filterCommunication() parameters
cellchat <- filterCommunication(cellchat, min.cells = 10, nonFilter.keep = TRUE)

3.1.7 Infer cell-cell communication at the signaling pathway level

Note: the inferred intercellular communication network of each ligand-receptor pair and each signaling pathway is stored in the slot net and netP, respectively.

# CellChat computes the communication probability at the signaling pathway level
# by summarizing the communication probabilities of all ligand-receptor interactions
# associated with each signaling pathway
cellchat <- computeCommunProbPathway(cellchat)

3.1.8 Calculate aggredated cell-cell communication network

By counting the number of links or summarizing the communication probability.

# you can use parameters `sources.use`, `targets.use`, `signaling`, and `pairLR.use` to only calculate the aggregated network among a subset of the cell-cell communication network
cellchat <- aggregateNet(cellchat)

3.1.9 Extract the inferred cellular communication network as a data frame

# there are a number of parameters used for filtering, see ?subsetCommunication
net_df <- subsetCommunication(cellchat)

3.1.10 Export the CellChat object together with the inferred cellular comunication networks

saveRDS(cellchat, file = gsub("\\.rda$", paste0(".", sample_name, ".0.rds"), rda_file))

3.2 Visualization of cell-cell communication network

Cell groups \(\rightarrow\) Signaling pathways \(\rightarrow\) L-R pairs.

3.2.1 Visualize the aggregated cell-cell communication network

3.2.1.1 Circle plot

CellChat 的原始推断结果是一个三维数组（来源细胞 × 靶向细胞 × 配受体对），而 cellchat@net$count 和 cellchat@net$weight 将这个三维数组压缩为二维矩阵（来源细胞 × 靶向细胞），其中的 count 和 weight 分别代表细胞群之间显著的配受体通讯数总和（by counting the number of links）以及细胞群之间显著的配受体通讯的总概率/强度（by summarizing the communication probability）。由于 CellChat 推断的细胞通讯是区分来源细胞和靶向细胞的，因此前述的 count 和 weight 均是非对称的二维矩阵。

行是来源细胞，列是靶向细胞。

# get the total number of cells for each cell group
# which will be used as the vertex weight
group_size <- as.numeric(table(cellchat@idents))

par(mfrow = c(1, 2), xpd = TRUE)

# each vertex indicates a cell group, and the edge width indicates the total number of interactions between two given cell groups
# 最宽的线对应着最大的值（`weight.scale`）
netVisual_circle(cellchat@net$count, vertex.weight = group_size, weight.scale = TRUE, label.edge = FALSE, title = "Number of interactions")

# each vertex indicates a cell group, and the edge width indicates the total interaction weights between two given cell groups
netVisual_circle(cellchat@net$weight, vertex.weight = group_size, weight.scale = TRUE, label.edge = FALSE, title = "Interaction weights/strength")

Only examine the signaling sent from each cell group and control the parameter edge.weight.max so we can compare edge weights between different networks:

mat <- cellchat@net$weight
par(mfrow = c(3, 4), xpd = TRUE)
for (i in seq_len(ncol(mat))) {
    mat2 <- matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))
    mat2[i, ] <- mat[i, ]
    netVisual_circle(mat2, vertex.weight = group_size, weight.scale = TRUE, edge.weight.max = max(mat), title.name = row.names(mat[i]))
}

3.2.1.2 Hierarchical plot

par(mfrow = c(1, 2), xpd = TRUE)
# 通过 `vertex.receiver` 指定 target cell groups 的 numeric indices，置于中间列。
# 左侧列的 cell groups 和中间列一致，作为 source cell groups。左侧列和中间列体现了感兴趣的细胞类型之间的 autocrine 和 paracrine。
# 右侧列为剩余的 cell groups，也是作为 source cell groups。右侧列和中间列体现了剩余细胞类型对感兴趣的细胞类型的 paracrine。
# 即中间列为 target cell groups，左右两列均为 source cell groups。
netVisual_hierarchy1(cellchat@net$weight, vertex.receiver = seq(1, 4))
# 这个函数的可视化结果相当于镜像翻转了前一函数的可视化结果。
netVisual_hierarchy2(cellchat@net$weight, vertex.receiver = seq(1, 4))

# 这个调用同时绘制两个图
# 且只绘制给定的 signaling pathway（上面的两个调用绘制所有的 signaling pathway）
# 左图用与左侧列一致的 cell groups 做 targets
# 右图用与右侧列一致的 cell groups 做 targets
netVisual_aggregate(cellchat, signaling = "MIF", vertex.receiver = seq(1, 4), layout = "hierarchy")

3.2.1.3 Chord diagram

# see `cellchat@netP$pathways` for all pathways
# cellchat@netP$prob[<source>, <target>, <pathway>]
show_pathway <- "CXCL"
par(mfrow = c(1, 2), xpd = TRUE)
# each sector is a cell group
# 一个细胞类型与另一个细胞类型通过该 signaling pathway 进行 interaction
netVisual_chord_cell(cellchat, signaling = show_pathway)
# each sector is a ligand, receptor, or signaling pathway
# 一个 cell group 通过哪个 ligand (e.g., CXCL12) 与另一个 cell group 通过哪个 receptor (e.g., CXCR4) 进行 interaction
netVisual_chord_gene(cellchat, signaling = show_pathway)

# 在 chord diagram 上，将属于同一大类的细胞类型放在邻近位置
group_cell_types <- c(rep("FIB", 4), rep("DC", 4), rep("TC", 4))
names(group_cell_types) <- levels(cellchat@idents)
netVisual_chord_cell(cellchat, signaling = show_pathway, group = group_cell_types, title.name = paste0(show_pathway, " signaling network"))

3.2.1.4 Heatmap

netVisual_heatmap(cellchat, signaling = show_pathway)

3.2.2 Visualize the inferred communication network of individual L-R pairs associated with that signaling pathway

3.2.2.1 Compute the contribution of each L-R pair to the overall signaling pathway

netAnalysis_contribution(cellchat, signaling = show_pathway)

3.2.2.2 Hierarchical plot

# extract all the significant L-R pairs for a given signaling pathways
LR_pairs <- extractEnrichedLR(cellchat, signaling = show_pathway)
# 针对特定的 L-R pair 绘制 hierarchical plot
netVisual_individual(cellchat, signaling = show_pathway, pairLR.use = LR_pairs[1, ], vertex.receiver = seq(1, 4), layout = "hierarchy")

3.2.2.3 Circle plot

netVisual_individual(cellchat, signaling = show_pathway, pairLR.use = LR_pairs[1, ], layout = "circle")

3.2.2.4 Chord diagram

netVisual_individual(cellchat, signaling = show_pathway, pairLR.use = LR_pairs[1, ], layout = "chord")

3.2.3 Vizualize cell-cell communication mediated by multiple L-R pairs or signaling pathways

3.2.3.1 Bubble plot

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# you can show only those significant L-R pairs associated with given signaling pathways specified by 'signaling'
# or those L-R pairs defined via 'pairLR.use'
netVisual_bubble(cellchat, sources.use = 4, targets.use = 5:11, remove.isolate = TRUE)

3.2.3.2 Chord diagram

See netVisual_chord_gene above.

netVisual_chord_gene(cellchat, sources.use = 4, targets.use = 5:8)

3.2.4 Systems analysis of cell-cell communication network

3.2.4.1 Identify signaling roles of cell groups and the major contributing signaling

# senders: out-degree metrics
# receivers: in-degree metrics
# mediators: flow betweenness metrics
# influencers: information centrality
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")

3.2.4.1.1 Identify major senders, receivers, mediators and influencers using heatmap

netAnalysis_signalingRole_network(cellchat, signaling = show_pathway)

3.2.4.1.2 Visualize dominant senders and receivers using scatter plot

# dot size is proportional to the number of inferred links (both outgoing and incoming) associated with each cell group
# dot shapes indicate different categories of cell groups if `group` is defined
p1 <- netAnalysis_signalingRole_scatter(cellchat)
p2 <- netAnalysis_signalingRole_scatter(cellchat, signaling = c("CXCL", "CCL"))
p1 + p2

3.2.4.1.3 Identify signals contributing the most to outgoing or incoming signaling of certain cell groups

p1 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing")
p2 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming")
p1 + p2

3.2.4.2 Identify global communication patterns to explore how multiple cell types and signaling pathways coordinate together

3.2.4.2.1 Identify and visualize outgoing communication pattern of secreting cells

Infer the number of patterns:

library(NMF)
library(ggalluvial)

# for a range of the number of patterns, a suitable number of patterns is the one at which Cophenetic and Silhouette values begin to drop suddenly
# in this example, it's at 6
p <- selectK(cellchat, pattern = "outgoing")
p

Heatmap:

n_patterns <- 6
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "outgoing", k = n_patterns)

River plot:

netAnalysis_river(cellchat, pattern = "outgoing")

Dot plot:

netAnalysis_dot(cellchat, pattern = "outgoing")

3.2.4.2.2 Identify and visualize incoming communication pattern of target cells

Infer the number of patterns:

p <- selectK(cellchat, pattern = "incoming")
p

Heatmap:

n_patterns <- 3
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "incoming", k = n_patterns)

River plot:

netAnalysis_river(cellchat, pattern = "incoming")

Dot plot:

netAnalysis_dot(cellchat, pattern = "incoming")

3.2.4.3 Manifold and classification learning analysis of signaling networks

Functional similarity: high degree of functional similarity indicates major senders and receivers are similar, and it can be interpreted as the two signaling pathways or two ligand-receptor pairs exhibit similar and/or redundant roles.

# Manifold and classification learning
cellchat <- computeNetSimilarity(cellchat, type = "functional")
cellchat <- netEmbedding(cellchat, type = "functional")
cellchat <- netClustering(cellchat, type = "functional")

netVisual_embedding(cellchat, type = "functional", label.size = 3.5)
netVisual_embeddingZoomIn(cellchat, type = "functional", nCol = 2)

Structural similarity: to compare their signaling network structure (singaling topology), without considering the similarity of senders and receivers.

# Manifold and classification learning
cellchat <- computeNetSimilarity(cellchat, type = "structural")
cellchat <- netEmbedding(cellchat, type = "structural")
cellchat <- netClustering(cellchat, type = "structural")

netVisual_embedding(cellchat, type = "structural", label.size = 3.5)
netVisual_embeddingZoomIn(cellchat, type = "structural", nCol = 2)

saveRDS(cellchat, file = gsub("\\.rda$", paste0(".", sample_name, ".1.rds"), rda_file))

4 CellChat for comparison analysis of multiple datasets with the same cell population compositions after joint clustering

4.1 Prepare the input data

cellchat_NL_rds_file <- "data/cellchat_humanSkin_NL.rds"
cellchat_LS_rds_file <- "data/cellchat_humanSkin_LS.rds"

# each dataset has been analyzed using CellChat previously
cellchat_NL <- readRDS(cellchat_NL_rds_file)
cellchat_LS <- readRDS(cellchat_LS_rds_file)

# 读入的 CellChat object 来源于版本低于 1.6.0 的 CellChat
# 因此需要进行版本提升
cellchat_NL <- updateCellChat(cellchat_NL)
cellchat_LS <- updateCellChat(cellchat_LS)

# merge multiple CellChat objects
cellchat_list <- list(NL = cellchat_NL, LS = cellchat_LS)
cellchat <- mergeCellChat(cellchat_list, add.names = names(cellchat_list))

4.2 Compare the total number of interactions and interaction strength

p1 <- compareInteractions(cellchat, show.legend = FALSE, group = c(1, 2), measure = "count")
p2 <- compareInteractions(cellchat, show.legend = FALSE, group = c(1, 2), measure = "weight")
p1 + p2

4.3 Compare the number of interactions and interaction strength among different cell populations

4.3.1 Compare the number of interactions and interaction strength among different cell populations

Red (or blue) colored edges represent increased (or decreased) signaling in the second dataset compared to the first one.

par(mfrow = c(1, 2), xpd = TRUE)
netVisual_diffInteraction(cellchat, weight.scale = TRUE, measure = "count")
netVisual_diffInteraction(cellchat, weight.scale = TRUE, measure = "weight")

4.3.2 Heatmap showing differential number of interactions or interaction strength among different cell populations across two datasets

The top colored bar plot represents the sum of each column of the absolute values displayed in the heatmap (incoming signaling).

The right colored bar plot represents the sum of each row of the absolute values (outgoing signaling).

Therefore, the bar height indicates the degree of change in terms of the number of interactions or interaction strength between the two conditions.

In the colorbar, red (or blue) represents increased (or decreased) signaling in the second dataset compared to the first one.

p1 <- netVisual_heatmap(cellchat, measure = "count")
p2 <- netVisual_heatmap(cellchat, measure = "weight")
p1 + p2

The above differential network analysis only works for pairwise datasets.

4.3.3 Circle plot showing the number of interactions or interaction strength among different cell populations across multiple datasets

对于多个数据集，将 edge weight 和 vertex weight 的最大值设为一致，从而使它们之间具有可比性。

weight_max <- getMaxWeight(cellchat_list, attribute = c("idents", "count"))
par(mfrow = c(1, 2), xpd = TRUE)
for (i in seq_along(cellchat_list)) {
    group_size <- as.numeric(table(cellchat_list[[i]]@idents))
    netVisual_circle(cellchat_list[[i]]@net$count, weight.scale = TRUE, label.edge = FALSE, edge.weight.max = weight_max[2], edge.width.max = 10, vertex.weight = group_size, vertex.weight.max = weight_max[1], vertex.weight.max = 10, title.name = paste0("Number of interactions - ", names(cellchat_list)[i]))
}

4.3.4 Circle plot showing the differential number of interactions or interaction strength among coarse cell types

Merge cell groups into coarser groups:

group_cell_type <- c(rep("FIB", 4), rep("DC", 4), rep("TC", 4))
group_cell_type <- factor(group_cell_type, levels = c("FIB", "DC", "TC"))

cellchat_list <- lapply(cellchat_list, function(x) mergeInteractions(x, group_cell_type))
cellchat <- mergeCellChat(cellchat_list, add.names = names(cellchat_list))

Once they were merged, we can use all visualization methods used above to explore the merged obejct, such as using the circle plot to visualize the differential number of interactions or interaction strength:

par(mfrow = c(1, 2), xpd = TRUE)
netVisual_diffInteraction(cellchat, weight.scale = TRUE, measure = "count.merged", label.edge = TRUE)
netVisual_diffInteraction(cellchat, weight.scale = TRUE, measure = "weight.merged", label.edge = TRUE)

4.4 Compare the major sources and targets

4.4.1 Identify cell populations with significant changes in sending or receiving signals

# 看每一个 cell group 的 sender or receiver 身份是否有所改变
# count the total number of associated outgoing (row) and ingoing (column) interactions (counting autocrine links only once) for each cell group
# 因为 autocrine links 被加了两次，因此减去一次
num_link <- sapply(cellchat_list, function(x) rowSums(x@net$count) + colSums(x@net$count) - diag(x@net$count))
# control the dot size in different datasets
# 统一缩放尺度
weight_min_max <- c(min(num_link), max(num_link))
gg <- list()
for (i in seq_along(cellchat_list)) {
    gg[[i]] <- netAnalysis_signalingRole_scatter(cellchat_list[[i]], title = names(cellchat_list)[i], weight.MinMax = weight_min_max)
}
patchwork::wrap_plots(plots = gg)

4.4.2 Identify the signaling changes of specific cell populations

如果一个 cell group 的 sender or receiver 身份有所改变，那是哪些 signaling pathways 发生了改变呢？

p1 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "Inflam. DC", signaling.exclude = "MIF")
p2 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "cDC1", signaling.exclude = c("MIF"))
patchwork::wrap_plots(plots = list(p1, p2))

4.5 Identify altered signaling with distinct network architecture and interaction strength

4.5.1 Identify signaling groups based on their functional similarity

Functional similarity: high degree of functional similarity indicates major senders and receivers are similar, and it can be interpreted as the two signaling pathways or two ligand-receptor pairs exhibit similar and/or redundant roles.

NB: functional similarity analysis is not applicable to multiple datsets with different cell type composition.

cellchat <- computeNetSimilarityPairwise(cellchat, type = "functional")
cellchat <- netEmbedding(cellchat, type = "functional")
cellchat <- netClustering(cellchat, type = "functional")

netVisual_embeddingPairwise(cellchat, type = "functional", label.size = 3.5)

4.5.2 Identify signaling groups based on structure similarity

Structural similarity: a structural similarity was used to compare their signaling network structure, without considering the similarity of senders and receivers.

NB: structural similarity analysis is applicable to multiple datsets with the same cell type composition or the vastly different cell type composition.

在 cell group 间进行信号转导的拓扑结构是否相似。

cellchat <- computeNetSimilarityPairwise(cellchat, type = "structural")
cellchat <- netEmbedding(cellchat, type = "structural")
cellchat <- netClustering(cellchat, type = "structural")

netVisual_embeddingPairwise(cellchat, type = "structural", label.size = 3.5)

4.5.3 Compute and visualize the pathway distance in the learned joint manifold

CellChat can identify the signaling networks with larger (or smaller) difference based on their Euclidean distance in the shared two-dimensions space.

Larger distance implies larger difference of the communication networks between two datasets in terms of either functional or structure similarity.

It should be noted that we only compute the distance of overlapped signaling pathways between two datasets. Those signaling pathways that are only identified in one dataset are not considered here.

If there are more than three datasets, one can do pairwise comparisons by modifying the parameter comparison in the function rankSimilarity.

rankSimilarity(cellchat, type = "functional")

rankSimilarity(cellchat, type = "structural")

4.6 Identify altered signaling with distinct interaction strength

By comparing the information flow/interaction strength of each signaling pathway, CellChat identifies signaling pathways that: (i) turn off, (ii) decrease, (iii) turn on, or (iv) increase, by changing their information flow at one condition as compared to another condition.

4.6.1 Compare the overall information flow of each signaling pathway or ligand-receptor pair

CellChat can identify the conserved and context-specific signaling pathways by simply comparing the information flow for each signaling pathway, which is defined by the sum of communication probability among all pairs of cell groups in the inferred network (i.e., the total weights in the network).

When setting do.stat = TRUE, a paired Wilcoxon test is performed to determine whether there is a significant difference of the signaling information flow between two conditions.

p1 <- rankNet(cellchat, mode = "comparison", measure = "weight", sources.use = NULL, targets.use = NULL, stacked = TRUE, do.stat = TRUE)
p2 <- rankNet(cellchat, mode = "comparison", measure = "weight", sources.use = NULL, targets.use = NULL, stacked = FALSE, do.stat = TRUE)

p1 + p2

4.6.2 Compare outgoing (or incoming) signaling patterns associated with each cell population

The above analysis summarize the information from the outgoing and incoming signaling together. CellChat can also compare the outgoing (or incoming) signaling pattern between two datasets, allowing to identify signaling pathways/ligand-receptors that exhibit different signaling patterns.

library(ComplexHeatmap)

i <- 1
union_pathways <- union(cellchat_list[[i]]@netP$pathways, cellchat_list[[i + 1]]@netP$pathways)

ht1 <- netAnalysis_signalingRole_heatmap(cellchat_list[[i]], pattern = "outgoing", signaling = union_pathways, title = names(cellchat_list)[i], width = 5, height = 6)
ht2 <- netAnalysis_signalingRole_heatmap(cellchat_list[[i + 1]], pattern = "outgoing", signaling = union_pathways, title = names(cellchat_list)[i + 1], width = 5, height = 6)
draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))

ht1 <- netAnalysis_signalingRole_heatmap(cellchat_list[[i]], pattern = "incoming", signaling = union_pathways, title = names(cellchat_list)[i], width = 5, height = 6, color.heatmap = "GnBu")
ht2 <- netAnalysis_signalingRole_heatmap(cellchat_list[[i + 1]], pattern = "incoming", signaling = union_pathways, title = names(cellchat_list)[i + 1], width = 5, height = 6, color.heatmap = "GnBu")
draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))

4.7 Identify the up-regulated and down-regulated signaling ligand-receptor pairs

4.7.1 Identify dysfunctional signaling by comparing the communication probabilities

CellChat can compare the communication probabilities mediated by L-R pairs from certain cell groups to other cell groups.

netVisual_bubble(cellchat, sources.use = 4, targets.use = 5:11, comparison = c(1, 2), angle.x = 45)

Moreover, CellChat can identify the up-regulated (increased) and down-regulated (decreased) signaling ligand-receptor pairs in one dataset compared to the other dataset.

# keeping the communications with highest probability in `max.dataset` (i.e., certrain condition)
p1 <- netVisual_bubble(cellchat, sources.use = 4, targets.use = 5:11, comparison = c(1, 2), max.dataset = 2, title.name = "Increased signaling in LS", angle.x = 45, remove.isolate = TRUE)
p2 <- netVisual_bubble(cellchat, sources.use = 4, targets.use = 5:11, comparison = c(1, 2), max.dataset = 1, title.name = "Decreased signaling in LS", angle.x = 45, remove.isolate = TRUE)
p1 + p2

4.7.2 Identify dysfunctional signaling by using differential expression analysis

The above method for identifying the up-regulated and down-regulated signaling is perfomed by comparing the communication probability between two datasets for each L-R pair and each pair of cell groups.

Alternatively, we can identify the up-regulated and down-regulated signaling ligand-receptor pairs based on the differential expression analysis (DEA). Specifically, we perform differential expression analysis between two biological conditions (i.e., NL and LS) for each cell group, and then obtain the up-regulated and down-regulated signaling based on the fold change of ligands in the sender cells and receptors in the receiver cells.

# define a positive dataset, i.e., the dataset with positive fold change against the other dataset
pos.dataset <- "LS"
# define a name used for storing the results of differential expression analysis
features.name <- paste0(pos.dataset, ".merged")

# perform differential expression analysis
# Of note, compared to CellChat version < v2, CellChat v2 now performs an ultra-fast Wilcoxon test using the presto package, which gives smaller values of logFC. Thus we here set a smaller value of thresh.fc compared to the original one (`thresh.fc = 0.1`). Users can also provide a vector and dataframe of customized DEGs by modifying the `cellchat@var.features$LS.merged` and `cellchat@var.features$LS.merged.info`.
cellchat <- identifyOverExpressedGenes(cellchat, group.dataset = "datasets", pos.dataset = pos.dataset, features.name = features.name, only.pos = FALSE, thresh.pc = 0.1, thresh.fc = 0.05, thresh.p = 0.05, group.DE.combined = FALSE)

# map the results of differential expression analysis onto the inferred cell-cell communications to easily manage/subset the ligand-receptor pairs of interest
net <- netMappingDEG(cellchat, features.name = features.name, variable.all = TRUE)

# extract the ligand-receptor pairs with up-regulated ligands in LS
net.up <- subsetCommunication(cellchat, net = net, datasets = "LS", ligand.logFC = 0.05, receptor.logFC = NULL)

# extract the ligand-receptor pairs with up-regulated ligands and up-regulated receptors in NL, i.e., down-regulated in LS
net.down <- subsetCommunication(cellchat, net = net, datasets = "NL", ligand.logFC = -0.05, receptor.logFC = -0.05)

Since the signaling genes in the net.up and net.down might be complex with multi-subunits, we can do further deconvolution to obtain the individual signaling genes.

gene.up <- extractGeneSubsetFromPair(net.up, cellchat)
gene.down <- extractGeneSubsetFromPair(net.down, cellchat)

# 查询特定的 features 在特定的 cell groups 中是否富集
# 此处仅看 outgoing
df <- findEnrichedSignaling(cellchat_list[[2]], features = c("CCL19", "CXCL12"), idents = c("Inflam. FIB", "COL11A1+ FIB"), pattern = "outgoing")

4.7.2.1 Visualize the identified up-regulated and down-regulated signaling ligand-receptor pairs

4.7.2.1.1 Bubble plot

pairLR.use.up <- unique(net.up[, "interaction_name", drop = FALSE])

p1 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.up, sources.use = 4, targets.use = 5:11, comparison = c(1, 2), angle.x = 90, remove.isolate = TRUE, title.name = paste0("Up-regulated signaling in ", names(cellchat_list)[2]))

pairLR.use.down <- unique(net.down[, "interaction_name", drop = FALSE])

p2 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.down, sources.use = c(4, 7), targets.use = 1, comparison = c(1, 2), angle.x = 90, remove.isolate = T, title.name = paste0("Down-regulated signaling in ", names(cellchat_list)[2]))

p1 + p2

4.7.2.1.2 Chord diagram

par(mfrow = c(1, 2), xpd = TRUE)
netVisual_chord_gene(cellchat_list[[2]], sources.use = 4, targets.use = 5:11, slot.name = "net", net = net.up, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(cellchat_list)[2]))
netVisual_chord_gene(cellchat_list[[1]], sources.use = c(4, 7), targets.use = 1, slot.name = "net", net = net.down, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Down-regulated signaling in ", names(cellchat_list)[2]))

4.7.2.1.3 Wordcloud plot

par(mfrow = c(1, 2), xpd = TRUE)
computeEnrichmentScore(net.up, species = "human", variable.both = TRUE)
computeEnrichmentScore(net.down, species = "human", variable.both = TRUE)

4.8 Visually compare cell-cell communication using Hierarchy plot, Circle plot or Chord diagram

Similar to the CellChat analysis of individual dataset, CellChat can visually compare cell-cell communication networks using hierarchy plot, circle plot, chord diagram, or heatmap. More details on the visualization can be found in the CellChat analysis of individual dataset.

Edge color/weight, node color/size/shape: in all visualization plots, edge colors are consistent with the sources as sender, and edge weights are proportional to the interaction strength. Thicker edge line indicates a stronger signal. In the Hierarchy plot and Circle plot, circle sizes are proportional to the number of cells in each cell group. In the hierarchy plot, solid and open circles represent source and target, respectively. In the Chord diagram, the inner thinner bar colors represent the targets that receive signal from the corresponding outer bar. The inner bar size is proportional to the signal strength received by the targets. Such inner bar is helpful for interpreting the complex chord diagram. Note that there exist some inner bars without any chord for some cell groups, please just igore it because this is an issue that has not been addressed by circlize package.

5 Comparison analysis of multiple datasets with different cell type compositions

5.1 Comparative analysis of multiple datasets with slightly different cell type compositions

For the datasets with differing cell type (group) compositions, CellChat lifts up the cell groups to the same cell type compositions across all datasets using the function liftCellChat, and then performs comparative analysis as the joint analysis of datasets with the same cell type compositions.

5.2 Comparison analysis of multiple datasets with vastly different cell type compositions

除了需要细胞类型对应的分析不能做，其它均可以做。