CellChat 细胞通讯分析（可视化）

<~生~信~交~流~与~合~作~请~关~注~公~众~号@生信探索>

预处理见上一篇推文https://mp.weixin.qq.com/s/ZsUQogkqcPXkaNDIV8GhWg，本篇内容是合并两个处理好的CellChat对象，然后进行对比分析和可视化，因为有许多细节需要手动调整所以就不写成脚本了。

A.加载R包、定义函数、建立文件夹

using(data.table, tidyverse, CellChat, patchwork, reticulate, ComplexHeatmap)
reticulate::use_python("/opt/homebrew/Caskroom/mambaforge/base/envs/SC/bin/python")
output_dir <- "CellChat"
setwd(output_dir)
group_names <- c("T", "C")
dirs <- list(
    "01.Overview", "02.Diff_number_strength", "03.Counts_Compare_select", "04.Compare_pathway_strengh",
    "05.Manifold", "06.Compare_Signal", "07.Compare_Pathways/Net",
    "07.Compare_Pathways/Chord", "07.Compare_Pathways/Heatmap", "08.Compare_LR_regulated",
    "09.GeneExpression", "10.Cell"
)

walk(dirs, ~ dir.create(.x, recursive = TRUE, showWarnings = FALSE))

ps <- function(filename, plot = FALSE, w = 12, h = 6) {
    if (is.object(plot)) {
        print(plot)
    }
    plot <- recordPlot()
    pdf(file = filename, onefile = T, width = w, height = h)
    replayPlot(plot)
    dev.off()
}

B. 合并对象

cc1 <- readRDS(str_glue("cc_{group_names[1]}.rds"))
cc2 <- readRDS(str_glue("cc_{group_names[2]}.rds"))
cc_list <- list(cc1, cc2)
names(cc_list) <- group_names
cc <- mergeCellChat(cc_list, add.names = group_names)
com_pathways <- purrr::map(cc_list, ~ .x@netP$pathways) %>%
    reduce(intersect) %>%
    unique()
pathway_union <- purrr::map(cc_list, ~ .x@netP$pathways) %>%
    reduce(union) %>%
    unique()

1.barplot: Overview_number_strength

p1 <- compareInteractions(cc, show.legend = FALSE, group = c(1, 2))
p2 <- compareInteractions(cc, show.legend = FALSE, group = c(1, 2), measure = "weight")
ps("01.Overview/Overview_number_strength.pdf", plot = p1 + p2, w = 6, h = 4)

2.数量与强度差异网络图、热图

par(mfrow = c(1, 2))
netVisual_diffInteraction(cc, weight.scale = T)
netVisual_diffInteraction(cc, weight.scale = T, measure = "weight")
ps("02.Diff_number_strength/net.pdf")

p1 <- netVisual_heatmap(cc, measure = "count")
p2 <- netVisual_heatmap(cc, measure = "weight")
ps("02.Diff_number_strength/heatmap.pdf", plot = p1 + p2)

par(mfrow = c(1, 2))
weight.max <- getMaxWeight(cc_list, attribute = c("idents", "count"))
for (i in seq(length(cc_list))) {
    netVisual_circle(cc_list[[i]]@net$count,
        weight.scale = T, label.edge = F,
        edge.weight.max = weight.max[2], edge.width.max = 12,
        title.name = paste0("Number of interactions - ", names(cc_list)[i])
    )
}
ps("02.Diff_number_strength/Counts_Compare_net.pdf")

3.指定细胞互作数量对比网络图

par(mfrow = c(1, 2))
s_cell <- c("CD4+ Cells", "CD8+ Cells", "Nk Cells")
count1 <- cc_list[[1]]@net$count[s_cell, s_cell]
count2 <- cc_list[[2]]@net$count[s_cell, s_cell]
weight.max <- max(max(count1), max(count2))
netVisual_circle(count1,
    weight.scale = T, label.edge = T, edge.weight.max = weight.max, edge.width.max = 12,
    title.name = paste0("Number of interactions-", names(cc_list)[1])
)
netVisual_circle(count2,
    weight.scale = T, label.edge = T, edge.weight.max = weight.max, edge.width.max = 12,
    title.name = paste0("Number of interactions-", names(cc_list)[2])
)
ps("03.Counts_Compare_select/circle.pdf", w = 10, h = 10)

4.保守和特异性信号通路的识别与可视化

p1 <- rankNet(cc, mode = "comparison", stacked = TRUE, do.stat = TRUE)
p2 <- rankNet(cc, mode = "comparison", stacked = FALSE, do.stat = TRUE)
ps("04.Compare_pathway_strengh/net.pdf", w = 10, h = 10, plot = p1 + p2)

5.流行学习识别差异信号通路

cc <- computeNetSimilarityPairwise(cc, type = "functional")
cc <- CellChat::netEmbedding(cc, type = "functional") # Python
cc <- netClustering(cc, type = "functional", do.parallel = FALSE)

netVisual_embeddingPairwise(cc, type = "functional", label.size = 2)
ps("05.Manifold/functional.pdf", w = 12, h = 10)

netVisual_embeddingPairwiseZoomIn(cc, type = "functional", nCol = 2)
ps("05.Manifold/functional_zoom.pdf", w = 18, h = 16)

rankSimilarity(cc, type = "functional") + ggtitle("functional similarity of pathway")
ps("05.Manifold/Pathway_functional_similarity.pdf", w = 9, h = 6)

cc <- computeNetSimilarityPairwise(cc, type = "structural")
cc <- netEmbedding(cc, type = "structural")
cc <- netClustering(cc, type = "structural", do.parallel = FALSE)

netVisual_embeddingPairwise(cc, type = "structural", label.size = 3.5)
ps("05.Manifold/structural.pdf", w = 12, h = 10)

netVisual_embeddingPairwiseZoomIn(cc, type = "structural", nCol = 2)
ps("05.Manifold/structural_zoom.pdf", w = 18, h = 16)

rankSimilarity(cc, type = "structural") + ggtitle("Structural similarity of pathway")
ps("05.Manifold/Pathway_structural_similarity.pdf", w = 9, h = 6)

6.细胞信号模式对比

for (x in c("all", "outgoing", "incoming")) {
    pl <- map(seq(length(cc_list)), ~ netAnalysis_signalingRole_heatmap(cc_list[[.x]], title = group_names[.x], pattern = x, signaling = pathway_union, width = 12, height = 28))
    p <- reduce(pl, `+`)
    ps(str_glue("06.Compare_Signal/{x}.pdf"), plot = p, w = 12, h = 13)
}

7. 特定信号通路的对比

for (x in com_pathways) {
    weight.max <- getMaxWeight(cc_list, slot.name = c("netP"), attribute = x)
    par(mfrow = c(1, 2), xpd = TRUE)
    for (y in seq(length(cc_list))) {
        netVisual_aggregate(cc_list[[y]],
            signaling = x,
            layout = "circle",
            edge.weight.max = weight.max[1],
            edge.width.max = 10,
            signaling.name = paste(x, group_names[y])
        )
    }
    ps(str_glue("07.Compare_Pathways/Net/{x}.pdf"), w = 12, h = 6)
}

for (x in com_pathways) {
    weight.max <- getMaxWeight(cc_list, slot.name = c("netP"), attribute = x)
    par(mfrow = c(1, 2), xpd = TRUE)
    for (y in seq(length(cc_list))) {
        p <- CellChat::netVisual_aggregate(cc_list[[y]],
            signaling = x,
            layout = "chord",
            pt.title = 3,
            title.space = 0.05,
            vertex.label.cex = 0.6,
            edge.weight.max = weight.max[1],
            edge.width.max = 10,
            signaling.name = paste(x, group_names[y])
        )
    }
    ps(str_glue("07.Compare_Pathways/Chord/{x}.pdf"), w = 9, h = 6, plot = p)
}

for (x in pathway_union) {
    tryCatch(
        {
            pl <- map(seq(length(cc_list)), ~ netVisual_heatmap(cc_list[[.x]], signaling = x, title.name = paste(x, "signaling ", group_names[y])))
            p <- reduce(pl, `+`)
            ps(str_glue("07.Compare_Pathways/Heatmap/{x}.pdf"), plot = p)
        },
        error = function(e) {
            print(str_glue("{x} error"))
        }
    )
}

8 气泡图展示High组上调或下调的配体受体对

for (x in pathway_union)
{
    tryCatch(
        {
            p <- netVisual_bubble(cc,
                comparison = c(1, 2), signaling = x,
                max.dataset = 1, title.name = str_glue("Increased signaling in {group_names[1]}"), angle.x = 45, remove.isolate = T
            )

            ps(str_glue("08.Compare_LR_regulated/{x}_Increased.pdf"), plot = p, w = 18, h = 6)
            p <- netVisual_bubble(cc,
                comparison = c(1, 2), signaling = x,
                min.dataset = 1, title.name = str_glue("Decreased signaling in {group_names[1]}"), angle.x = 45, remove.isolate = T
            )
            ps(str_glue("08.Compare_LR_regulated/{x}_Decreased.pdf"), plot = p, w = 18, h = 6)
        },
        error = function(e) {
            print(str_glue("{x} error"))
        }
    )
}

# perform differential expression analysis
cc <- identifyOverExpressedGenes(cc, group.dataset = "datasets", pos.dataset = group_names[1], features.name = group_names[1], only.pos = FALSE, thresh.pc = 0.1, thresh.fc = 0.1, thresh.p = 1)
#> Use the joint cell labels from the merged CellChat object
# 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(cc, features.name = group_names[1])
# extract the ligand-receptor pairs with upregulated ligands in LS
net.up <- subsetCommunication(cc, net = net, datasets = group_names[1], ligand.logFC = 0.3, receptor.logFC = 0.3)
net.down <- subsetCommunication(cc, net = net, datasets = group_names[2],ligand.logFC = -0.1, receptor.logFC = NULL)
# extract the ligand-receptor pairs with upregulated ligands and upregulated recetptors in NL, i.e.,downregulated in LS
gene.up <- extractGeneSubsetFromPair(net.up, cc)
gene.down <- extractGeneSubsetFromPair(net.down, cc)

pairLR.use.up = net.up[, "interaction_name", drop = F]
netVisual_bubble(cc, pairLR.use = pairLR.use.up, max.dataset = 1, comparison = c(1, 2),  angle.x = 90, remove.isolate = T,title.name = paste0("Up-regulated signaling in ", names(cc_list)[1]))
ps("08.Compare_LR_regulated/0up_filtered.pdf",  w = 30, h = 12)
pairLR.use.down = net.down[, "interaction_name", drop = F]
netVisual_bubble(cc, pairLR.use = pairLR.use.down, comparison = c(1, 2), min.dataset = 1,  angle.x = 90, remove.isolate = T,title.name = paste0("Down-regulated signaling in ", names(cc_list)[1]))
ps("08.Compare_LR_regulated/0down_filtered.pdf",  w = 30, h = 12)

9 比较不同数据集之间地信号基因表达分布

cc@meta$datasets <- factor(cc@meta$datasets, levels = group_names) # set factor level
for (x in pathway_union) {
    tryCatch(
        {
            p <- plotGeneExpression(cc, signaling = x, split.by = "datasets", colors.ggplot = T)
            ps(str_glue("09.GeneExpression/{x}.pdf"), w = 12, h = 9, plot = p)
        },
        error = function(e) {
            print(str_glue("{x} error"))
        }
    )
}

10 比较二维空间中的主要源和目标

for (x in unique(cc@meta$CellType)) {
    tryCatch(
        {
            p <- netAnalysis_signalingChanges_scatter(cc, idents.use = x)
            ps(str_glue("10.Cell/{x}.pdf"), w = 9, h = 6, plot = p)
        },
        error = function(e) {
            print(str_glue("{x} error"))
        }
    )
}

num.link <- sapply(cc_list, function(x) {
    rowSums(x@net$count) + colSums(x@net$count) - diag(x@net$count)
})
weight.MinMax <- c(min(num.link), max(num.link)) # control the dot size in the different datasets
pl <- purrr::map(seq(length(cc_list)), ~ netAnalysis_signalingRole_scatter(cc_list[[.x]], title = names(cc_list)[.x], weight.MinMax = weight.MinMax))
patchwork::wrap_plots(plots = pl)
ps("10.Cell/scater.pdf", w = 12, h = 6)

工作目录

├── 01.Overview
├── 02.Diff_number_strength
├── 03.Counts_Compare_select
├── 04.Compare_pathway_strengh
├── 05.Manifold
├── 06.Compare_Signal
├── 07.Compare_Pathways
├── 08.Compare_LR_regulated
├── 09.GeneExpression
├── 10.Cell
├── C.arrow
├── C.csv
├── CellChat_1.R
├── CellChat_2.R
├── T.arrow
├── T.csv
├── cc_C.rds
├── cc_T.rds
├── estimationNumCluster_C_functional_dataset_single.pdf
├── estimationNumCluster_C_structural_dataset_single.pdf
├── estimationNumCluster_T_functional_dataset_single.pdf
├── estimationNumCluster_T_structural_dataset_single.pdf
├── estimationNumCluster__functional_dataset_1-2.pdf
└── estimationNumCluster__structural_dataset_1-2.pdf

Reference

https://www.jianshu.com/p/da145cff3d41
https://www.jianshu.com/p/b3d26ac51c5a
https://cloud.tencent.com/developer/inventory/26535/article/1935670