CellChat signalling analysis of Arcuate Nucleus/Median Eminence Lxn (+/-) astrocytes dataset from Lutomska LM et al 2022

Load data and setup parameters

Code

# Load tidyverse infrastructure packages
suppressPackageStartupMessages({
  library(future)
  library(here)
  library(tidyverse)
  library(magrittr)
  library(stringr)
  library(skimr)
  library(RColorBrewer)
  library(viridis)
})


# Load packages for scRNA-seq analysis and visualisation
suppressPackageStartupMessages({
  library(ggplot2)
  library(cowplot)
  library(patchwork)
  library(ggstatsplot)
  library(ComplexHeatmap)
  library(sceasy)
  library(anndata)
  library(Seurat)
  library(SeuratDisk)
  library(SeuratWrappers)
  library(CellChat)
  library(scCustomize)
})

sc <- import("scanpy", convert = FALSE)

Set paths

Code

src_dir <- here("code")
data_dir <- here("data")
output_dir <- here("output")
plots_dir <- here(output_dir, "figures/")
tables_dir <- here(output_dir, "tables/")

Load helper functions and gene-sets

Code

source(here(src_dir, "genes.R"))
source(here(src_dir, "functions.R"))

Set fixed variables

Code

# set seed
reseed <- 42
set.seed(seed = reseed)

# Parameters for parallel execution
n_cores <- 8
plan("multisession", workers = n_cores)
options(
  future.globals.maxSize = 200000 * 1024^2,
  future.rng.onMisuse = "ignore"
)
plan()

multisession:
- args: function (..., workers = 8, envir = parent.frame())
- tweaked: TRUE
- call: plan("multisession", workers = n_cores)

Code

# ggplot2 theme
theme_set(ggmin::theme_powerpoint())

Code

bioproject <- "PRJNA847050"
project <- "lutomska2022-Arc"
cb_fpr <- 0.001
low_cutoff_gene <- 500
high_cutoff_gene <- NULL
high_cutoff_gene <- 10000
low_cutoff_umis <- NULL
low_cutoff_umis <- -Inf
high_cutoff_umis <- 45000
high_cutoff_pc_mt <- 25
high_cutoff_pc_ribo <- 20
high_cutoff_pc_hb <- 0.1
high_cutoff_doublet_score <- 0.33
high_cutoff_complexity <- 0.8

Load selected astrocytes data from Lutomska LM et al (2022)

Code

anndata <- sc$read(here(
  data_dir,
  sprintf("class_cello/%s-whole_dataset-%s-cello_annotation.h5ad", bioproject, cb_fpr)
))

astrodata <- sc$read(here(
  data_dir,
  sprintf("class_cello/%s-astrocytes_dataset-%s-initial_selection.h5ad", bioproject, cb_fpr)
))

astrocytes_idx <- py_to_r(astrodata$obs_names)

Convert adata object to R AnnDataR6 object.

Code

adata <- py_to_r(anndata)
class(adata)

[1] "AnnDataR6" "R6"       

Code

class(adata$X)

[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"

Code

adata

AnnData object with n_obs × n_vars = 14594 × 24746
    obs: 'nCount_RAW', 'nFeature_RAW', 'nCount_RNA', 'nFeature_RNA', 'orig.ident', 'nFeature_Diff', 'nCount_Diff', 'percent_mito', 'percent_ribo', 'percent_mito_ribo', 'percent_hb', 'log10GenesPerUMI', 'cell_name', 'barcode', 'latent_RT_efficiency', 'latent_cell_probability', 'latent_scale', 'doublet_score', 'predicted_doublets', 'QC', 'var_regex', 'RNA_snn_res.0.5', 'RNA_snn_res.0.71728703859461', 'RNA_snn_res.1.18939286532015', 'RNA_snn_res.3.00000100000001', 'seurat_clusters', 'k_tree', 'comb_clstr1', 'S.Score', 'G2M.Score', 'Phase', 'nCount_SCT', 'nFeature_SCT', 'SCT_snn_res.1', 'SCT_snn_res.1.10031502002514', 'SCT_snn_res.1.2198400795032', 'SCT_snn_res.1.36467746508232', 'SCT_snn_res.1.54382016547197', 'SCT_snn_res.1.77109349840298', 'SCT_snn_res.2.0689045782524', 'SCT_snn_res.2.47612417203789', 'SCT_snn_res.3.06661007679595', 'SCT_snn_res.4.00000100000002', 'bioproject', 'project', 'model', 'tech', 'region', 'sex', 'stage', 'libname', 'expbtch', 'condit', 'ora_celltype'
    var: 'vst.mean', 'vst.variance', 'vst.variance.expected', 'vst.variance.standardized', 'vst.variable'
    uns: 'k_tree_colors', 'name', 'ora_celltype_colors'
    obsm: 'X_pacmap', 'X_pca', 'X_umap', 'ora_estimate', 'ora_pvals'

Code

srt_path <- here(
  data_dir,
  sprintf("%s-whole_dataset-%s-cello_annotation.h5ad", bioproject, cb_fpr)
)

expr_mtx <- t(as.matrix(adata$raw$X))
colnames(expr_mtx) <- rownames(adata$X)
rownames(expr_mtx) <- adata$var_names
srt <- CreateSeuratObject(
  expr_mtx,
  assay = "RNA",
  project = "individual_hypothalamic_nuclei_astrocytes_evaluation_dataset",
  meta.data = as.data.frame(adata$obs),
  row.names = colnames(adata$X)
)

Idents(srt) <- "ora_celltype"
srt <- subset(srt, idents = c("Astrocytes", "Neurons"))

Idents(srt) <- "expbtch"
srt_keep <- srt
srt <- subset(srt, idents = c("control", "hfd_05"))
srt$libname <- factor(srt$expbtch, levels = c("control" = "control", "HFHSD" = "hfd_05"))
Idents(srt) <- "libname"
Idents(srt) <- "k_tree"

print(srt)

An object of class Seurat 
24746 features across 4057 samples within 1 assay 
Active assay: RNA (24746 features, 0 variable features)

Code

rm(adata, astrodata, anndata, expr_mtx)
invisible(gc())

Code

Idents(srt, cells = WhichCells(srt, expression = (Pomc > 10 & ora_celltype == "Neurons"))) <- "POMC neurons"
Idents(srt, cells = WhichCells(srt, expression = ((Agrp > 10 | Npy > 10) & ora_celltype == "Neurons"))) <- "AGRP/NPY neurons"
Idents(srt, cells = WhichCells(srt, expression = (Lxn > 0 & cell_name %in% astrocytes_idx))) <- "Astro_LXN+"
Idents(srt, cells = WhichCells(srt, expression = (Lxn == 0 & cell_name %in% astrocytes_idx))) <- "Astro_LXN-"

srt <- subset(srt, idents = c("Astro_LXN+", "Astro_LXN-", "POMC neurons", "AGRP/NPY neurons"))
srt$goi <- Idents(srt)

Code

srt.list <- SplitObject(srt, split.by = "libname")
srt.list %<>% map(NormalizeData)
srt.list

$hfd_05
An object of class Seurat 
24746 features across 1111 samples within 1 assay 
Active assay: RNA (24746 features, 0 variable features)

$control
An object of class Seurat 
24746 features across 1201 samples within 1 assay 
Active assay: RNA (24746 features, 0 variable features)

Code

SrtCellChat <- function(seurat_object) {
  data.input <- GetAssayData(seurat_object, assay = "RNA", slot = "data") # normalized data matrix
  labels <- Idents(seurat_object)
  meta <- data.frame(labels = labels, row.names = names(labels)) # create a dataframe of the cell labels
  cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels")
  return(cellchat)
}

object.list <- lapply(X = srt.list, FUN = SrtCellChat)

[1] "Create a CellChat object from a data matrix"
Set cell identities for the new CellChat object 
The cell groups used for CellChat analysis are  Astro_LXN- Astro_LXN+ AGRP/NPY neurons POMC neurons 
[1] "Create a CellChat object from a data matrix"
Set cell identities for the new CellChat object 
The cell groups used for CellChat analysis are  Astro_LXN- Astro_LXN+ AGRP/NPY neurons POMC neurons 

Code

CellChatDB <- CellChatDB.mouse # running on mouse data
showDatabaseCategory(CellChatDB)

Code

glimpse(CellChatDB$interaction)

Rows: 2,019
Columns: 11
$ interaction_name   <chr> "TGFB1_TGFBR1_TGFBR2", "TGFB2_TGFBR1_TGFBR2", "TGFB…
$ pathway_name       <chr> "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TG…
$ ligand             <chr> "Tgfb1", "Tgfb2", "Tgfb3", "Tgfb1", "Tgfb1", "Tgfb2…
$ receptor           <chr> "TGFbR1_R2", "TGFbR1_R2", "TGFbR1_R2", "ACVR1B_TGFb…
$ agonist            <chr> "TGFb agonist", "TGFb agonist", "TGFb agonist", "TG…
$ antagonist         <chr> "TGFb antagonist", "TGFb antagonist", "TGFb antagon…
$ co_A_receptor      <chr> "", "", "", "", "", "", "", "", "", "", "", "", "",…
$ co_I_receptor      <chr> "TGFb inhibition receptor", "TGFb inhibition recept…
$ evidence           <chr> "KEGG: mmu04350", "KEGG: mmu04350", "KEGG: mmu04350…
$ annotation         <chr> "Secreted Signaling", "Secreted Signaling", "Secret…
$ interaction_name_2 <chr> "Tgfb1 - (Tgfbr1+Tgfbr2)", "Tgfb2  - (Tgfbr1+Tgfbr2…

Code

skimr::skim(CellChatDB$interaction)

Data summary

Name	CellChatDB$interaction
Number of rows	2019
Number of columns	11
_______________________
Column type frequency:
character	11
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
interaction_name	0	1	6	23	0	2019	0
pathway_name	0	1	2	12	0	229	0
ligand	0	1	2	14	0	577	0
receptor	0	1	3	18	0	510	0
agonist	0	1	0	15	1676	5	0
antagonist	0	1	0	18	1563	11	0
co_A_receptor	0	1	0	25	1637	9	0
co_I_receptor	0	1	0	27	1530	11	0
evidence	0	1	10	45	0	169	0
annotation	0	1	12	18	0	3	0
interaction_name_2	0	1	9	26	0	2019	0

Code

janitor::tabyl(CellChatDB$interaction, annotation)

Code

CellChatDB.use <- CellChatDB # simply use the default CellChatDB

Code

invisible(gc())
plan("sequential")
invisible(gc())
plan("multisession", workers = n_cores)

if (file.exists(here(data_dir, "cellchat-lutomska2022-HFHSD5_and_control-lxn.rds"))) {
    object.list <- readr::read_rds(file = here(data_dir, "cellchat-lutomska2022-HFHSD5_and_control-lxn.rds"))
} else {
  # set the used database in the object
  object.list$hfd_05@DB <- CellChatDB.use
  object.list$control@DB <- CellChatDB.use

  object.list %<>% map(subsetData)
  object.list %<>% map(identifyOverExpressedGenes)
  object.list %<>% map(identifyOverExpressedInteractions)
  object.list %<>% map(projectData, PPI.mouse)
  object.list %<>% map(computeCommunProb, raw.use = FALSE)
  object.list %<>% map(filterCommunication, min.cells = 5)
  object.list %<>% map(computeCommunProbPathway)
  object.list %<>% map(aggregateNet)
  object.list %<>% map(netAnalysis_computeCentrality)
  readr::write_rds(object.list, file = here(data_dir, "cellchat-lutomska2022-HFHSD5_and_control-lxn.rds"))
}

Plot CellChat network of HFHSD and control mice

Control

Code

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

Code

mat <- cellchat@net$weight
par(mfrow = c(2, 2), xpd = TRUE)
for (i in 1:nrow(mat)) {
  mat2 <- matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))
  mat2[i, ] <- mat[i, ]
  netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat), title.name = rownames(mat)[i])
}

HFHSD

Code

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

Code

mat <- cellchat@net$weight
par(mfrow = c(2, 2), xpd = TRUE)
for (i in 1:nrow(mat)) {
  mat2 <- matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))
  mat2[i, ] <- mat[i, ]
  netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat), title.name = rownames(mat)[i])
}

Compare CellChat network of HFHSD and control mice

Compare the total number of interactions and interaction strength

Code

object.list <- list(control = object.list$control, HFHSD = object.list$hfd_05)

cellchat <- mergeCellChat(object.list, add.names = names(object.list))

gg1 <- compareInteractions(cellchat, show.legend = F, group = c(1, 2))
gg2 <- compareInteractions(cellchat, show.legend = F, group = c(1, 2), measure = "weight")
gg1 + gg2

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

Differential number of interactions or interaction strength among different cell populations

Code

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

Code

gg1 <- netVisual_heatmap(cellchat)
gg2 <- netVisual_heatmap(cellchat, measure = "weight")
gg1 + gg2

Code

weight.max <- getMaxWeight(object.list, attribute = c("idents", "count"))
par(mfrow = c(1, 2), xpd = TRUE)
for (i in 1:length(object.list)) {
  netVisual_circle(object.list[[i]]@net$count, weight.scale = T, label.edge = T, edge.weight.max = weight.max[2], edge.width.max = 12, title.name = paste0("Number of interactions - ", names(object.list)[i]))
}

Compare the major sources and targets in 2D space

Code

num.link <- sapply(object.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
gg <- list()
for (i in 1:length(object.list)) {
  gg[[i]] <- netAnalysis_signalingRole_scatter(object.list[[i]], title = names(object.list)[i], weight.MinMax = weight.MinMax)
}

patchwork::wrap_plots(plots = gg)

Code

gg1 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "Astro_LXN+")
gg2 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "Astro_LXN-")
gg3 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "POMC neurons")
gg4 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "AGRP/NPY neurons")
patchwork::wrap_plots(plots = list(gg1, gg2, gg3, gg4))

Identify the conserved and context-specific signaling pathways

Identify signaling networks with larger (or less) difference as well as signaling groups based on their functional/structure 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.

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.

Identify signaling groups based on their functional similarity

Code

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

Compute signaling network similarity for datasets 1 2 

Code

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

Manifold learning of the signaling networks for datasets 1 2 

Code

cellchat <- netClustering(cellchat, type = "functional", do.parallel = FALSE)

Classification learning of the signaling networks for datasets 1 2 

Code

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

2D visualization of signaling networks from datasets 1 2 

Identify signaling groups based on structure similarity

Code

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

Compute signaling network similarity for datasets 1 2 

Code

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

Manifold learning of the signaling networks for datasets 1 2 

Code

cellchat <- netClustering(cellchat, type = "structural", do.parallel = FALSE)

Classification learning of the signaling networks for datasets 1 2 

Code

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

2D visualization of signaling networks from datasets 1 2 

Compute and visualize the pathway distance in the learned joint manifold

We can identify the signaling networks with larger (or less) 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. NB: 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 defining comparison in the function rankSimilarity.

Code

rankSimilarity(cellchat, type = "functional")

Compute the distance of signaling networks between datasets 1 2 

Identify and visualize the conserved and context-specific signaling pathways

By comparing the information flow/interaction strengh of each signaling pathway, we can identify signaling pathways, (i) turn off, (ii) decrease, (iii) turn on or (iv) increase, by change their information flow at one condition as compared to another condition.

Compare the overall information flow of each signaling pathway

We 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).

This bar graph can be plotted in a stacked mode or not. Significant signaling pathways were ranked based on differences in the overall information flow within the inferred networks between HFHSD and control mice arcuate nucleus of hypothalamus. The top signaling pathways colored red are enriched in HFHSD mice, and these colored blue were enriched in the control mice.

Code

gg1 <- rankNet(cellchat, mode = "comparison", stacked = T, do.stat = TRUE)
gg2 <- rankNet(cellchat, mode = "comparison", stacked = F, do.stat = TRUE)
gg1 + gg2

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

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

We can combine all the identified signaling pathways from different datasets and thus compare them side by side, including outgoing signaling, incoming signaling and overall signaling by aggregating outgoing and incoming signaling together. NB: rankNet also shows the comparison of overall signaling, but it does not show the signaling strength in specific cell populations.

Code

i <- 1
# combining all the identified signaling pathways from different datasets
pathway.union <- union(object.list[[i]]@netP$pathways, object.list[[i + 1]]@netP$pathways)
ht1 <- netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "outgoing", signaling = pathway.union, title = names(object.list)[i], width = 5, height = 14)
ht2 <- netAnalysis_signalingRole_heatmap(object.list[[i + 1]], pattern = "outgoing", signaling = pathway.union, title = names(object.list)[i + 1], width = 5, height = 14)
draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))

Code

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

Code

ht1 <- netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "all", signaling = pathway.union, title = names(object.list)[i], width = 5, height = 14, color.heatmap = "OrRd")
ht2 <- netAnalysis_signalingRole_heatmap(object.list[[i + 1]], pattern = "all", signaling = pathway.union, title = names(object.list)[i + 1], width = 5, height = 14, color.heatmap = "OrRd")
draw(ht2 + ht1, ht_gap = unit(0.5, "cm"))

Identify dysfunctional signaling by using differential expression analysis

The above method for identifying the upgulated 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. Alternative, we can identify the upgulated and down-regulated signaling ligand-receptor pairs based on the differential gene expression analysis. Specifically, we perform differential expression analysis between two biological conditions (HFHSD and control mice) for each cell group, and then obtain the upgulated and down-regulated signaling based on the fold change of ligands in the sender cells and receptors in the receiver cells. Such analysis can be done as follows.

Code

# define a positive dataset, i.e., the dataset with positive fold change against the other dataset
pos.dataset <- "HFHSD"
# define a char name used for storing the results of differential expression analysis
features.name <- pos.dataset
# perform differential expression analysis
cellchat <- identifyOverExpressedGenes(cellchat, group.dataset = "datasets", pos.dataset = pos.dataset, features.name = features.name, 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(cellchat, features.name = features.name)
# extract the ligand-receptor pairs with upregulated ligands in HFHSD
net.up <- subsetCommunication(cellchat, net = net, datasets = "HFHSD", ligand.logFC = 0.2, receptor.logFC = NULL)
# extract the ligand-receptor pairs with upregulated ligands and upregulated recetptors in control, i.e.,downregulated in HFHSD
net.down <- subsetCommunication(cellchat, net = net, datasets = "control", ligand.logFC = -0.1, receptor.logFC = -0.1)

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

Code

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

Code

pairLR.use.up <- net.up[, "interaction_name", drop = F]
gg1 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.up, sources.use = c(1, 2), targets.use = c(3, 4), comparison = c(1, 2), angle.x = 90, remove.isolate = T, title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
#> Comparing communications on a merged object
pairLR.use.down <- net.down[, "interaction_name", drop = F]
gg2 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.down, sources.use = c(1, 2), targets.use = c(3, 4), comparison = c(1, 2), angle.x = 90, remove.isolate = T, title.name = paste0("Down-regulated signaling in ", names(object.list)[2]))
#> Comparing communications on a merged object
gg1 + gg2

Code

gg1 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.up, sources.use = c(3, 4), targets.use = c(1, 2), comparison = c(1, 2), angle.x = 90, remove.isolate = T, title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
#> Comparing communications on a merged object
pairLR.use.down <- net.down[, "interaction_name", drop = F]
gg2 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.down, sources.use = c(3, 4), targets.use = c(1, 2), comparison = c(1, 2), angle.x = 90, remove.isolate = T, title.name = paste0("Down-regulated signaling in ", names(object.list)[2]))
#> Comparing communications on a merged object
gg1 + gg2

Visualize the enriched ligands, signaling,or ligand-receptor pairs in one condition compared to another condition using wordcloud

Code

computeEnrichmentScore(net.down, species = "mouse")

Visually compare cell-cell communication using Hierarchy plot, Circle plot

Code

pathways.show.all <- intersect(object.list$HFHSD@netP$pathways, object.list$control@netP$pathways)
for (i in 1:length(pathways.show.all)) {
  pathways.show <- pathways.show.all[i]
  ht <- list()
  tryCatch(
    {
        for (j in 1:length(object.list)) {
            ht[[j]] <- netVisual_heatmap(object.list[[j]], signaling = pathways.show, color.heatmap = "Reds", title.name = paste(pathways.show, "signaling ", names(object.list)[j]))
        }
        ComplexHeatmap::draw(ht[[1]] + ht[[2]], ht_gap = unit(0.5, "cm"))
    },
    error = function(e) {
      # If an error occurs, silently skip to the next iteration
      next
    }
  )
}

Code

for (i in 1:length(pathways.show.all)) {
  tryCatch(
    {
        pathways.show <- pathways.show.all[i]
        weight.max <- getMaxWeight(object.list, slot.name = c("netP"), attribute = pathways.show) # control the edge weights across different datasets
        par(mfrow = c(1, 2), xpd = TRUE)
        for (j in 1:length(object.list)) {
            netVisual_aggregate(object.list[[j]], signaling = pathways.show, label.edge = T, layout = "circle", edge.weight.max = weight.max[1], edge.width.max = 10, signaling.name = paste(pathways.show, names(object.list)[j]))
        }
    },
    error = function(e) {
      # If an error occurs, silently skip to the next iteration
      next
    }
  )
}

Compare the signaling gene expression distribution between different datasets

Code

cellchat@meta$datasets <- factor(cellchat@meta$datasets, levels = c("control", "HFHSD")) # set factor level
for (i in 1:length(pathways.show.all)) {
  tryCatch(
    {
        pathways.show <- pathways.show.all[i]
        plotGeneExpression(cellchat, signaling = pathways.show, split.by = "datasets", colors.ggplot = T)
    },
    error = function(e) {
      # If an error occurs, silently skip to the next iteration
      next
    }
  )
}

Save the merged CellChat object

Code

readr::write_rds(cellchat, file = here(data_dir, "cellchat-comparison-lutomska2022-HFHSD5_vs_control-lxn.rds"))

Session information

Code

sI <- sessioninfo::session_info()
sI$loadedOnly <- NULL
print(sI, locale = FALSE)

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.1 (2023-06-16)
 os       Ubuntu 22.04.2 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       Etc/UTC
 date     2023-08-01
 pandoc   2.19.2 @ /home/etretiakov/micromamba/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package          * version     date (UTC) lib source
 abind              1.4-5       2016-07-21 [1] RSPM (R 4.2.2)
 anndata          * 0.7.5.6     2023-03-17 [1] RSPM (R 4.3.0)
 assertthat         0.2.1       2019-03-21 [1] CRAN (R 4.2.2)
 backports          1.4.1       2021-12-13 [1] RSPM (R 4.2.2)
 base64enc          0.1-3       2015-07-28 [1] RSPM (R 4.2.2)
 bayestestR         0.13.1      2023-04-07 [1] RSPM (R 4.3.0)
 beeswarm           0.4.0       2021-06-01 [1] RSPM (R 4.2.2)
 Biobase          * 2.60.0      2023-04-25 [1] Bioconductor
 BiocGenerics     * 0.46.0      2023-04-25 [1] Bioconductor
 BiocManager        1.30.21.1   2023-07-18 [1] RSPM (R 4.3.1)
 BiocNeighbors      1.18.0      2023-04-25 [1] Bioconductor
 BiocParallel       1.34.2      2023-05-22 [1] Bioconductor
 bit                4.0.5       2022-11-15 [1] RSPM (R 4.2.2)
 bit64              4.0.5       2020-08-30 [1] RSPM (R 4.2.2)
 broom              1.0.5       2023-06-09 [1] RSPM (R 4.3.0)
 Cairo              1.6-0       2022-07-05 [1] RSPM (R 4.2.2)
 car                3.1-2       2023-03-30 [1] RSPM (R 4.3.0)
 carData            3.0-5       2022-01-06 [1] RSPM (R 4.2.2)
 CellChat         * 1.6.1       2023-07-30 [1] Github (sqjin/CellChat@7ee94c3)
 circlize           0.4.15      2022-05-10 [1] RSPM (R 4.2.2)
 cli                3.6.1       2023-03-23 [1] RSPM (R 4.3.0)
 clue               0.3-64      2023-01-31 [1] RSPM (R 4.3.0)
 cluster            2.1.4       2022-08-22 [3] CRAN (R 4.2.1)
 coda               0.19-4      2020-09-30 [1] RSPM (R 4.2.2)
 codetools          0.2-19      2023-02-01 [3] CRAN (R 4.2.2)
 colorspace         2.1-0       2023-01-23 [1] RSPM (R 4.3.0)
 ComplexHeatmap   * 2.16.0      2023-04-25 [1] Bioconductor
 correlation        0.8.4       2023-04-28 [1] Github (easystats/correlation@3d5cd1e)
 cowplot          * 1.1.1       2020-12-30 [1] RSPM (R 4.2.2)
 crayon             1.5.2       2022-09-29 [1] RSPM (R 4.2.2)
 data.table         1.14.8      2023-02-17 [1] RSPM (R 4.3.0)
 datawizard         0.8.0       2023-06-16 [1] RSPM (R 4.3.0)
 DBI                1.1.3       2022-06-18 [1] RSPM (R 4.2.2)
 deldir             1.0-9       2023-05-17 [1] RSPM (R 4.3.0)
 digest             0.6.33      2023-07-07 [1] RSPM (R 4.3.1)
 doParallel         1.0.17      2022-02-07 [1] RSPM (R 4.2.2)
 dplyr            * 1.1.2       2023-04-20 [1] RSPM (R 4.3.0)
 ellipsis           0.3.2       2021-04-29 [1] RSPM (R 4.2.2)
 emmeans            1.8.6       2023-05-11 [1] RSPM (R 4.3.0)
 estimability       1.4.1       2022-08-05 [1] RSPM (R 4.3.0)
 evaluate           0.21        2023-05-05 [1] RSPM (R 4.3.0)
 fansi              1.0.4       2023-01-22 [1] RSPM (R 4.3.0)
 farver             2.1.1       2022-07-06 [1] RSPM (R 4.2.2)
 fastmap            1.1.1       2023-02-24 [1] RSPM (R 4.3.0)
 fitdistrplus       1.1-11      2023-04-25 [1] RSPM (R 4.3.0)
 FNN                1.1.3.2     2023-03-20 [1] RSPM (R 4.3.0)
 forcats          * 1.0.0       2023-01-29 [1] RSPM (R 4.3.0)
 foreach            1.5.2       2022-02-02 [1] RSPM (R 4.2.2)
 future           * 1.33.0      2023-07-01 [1] RSPM (R 4.3.1)
 future.apply       1.11.0      2023-05-21 [1] RSPM (R 4.3.0)
 generics           0.1.3       2022-07-05 [1] RSPM (R 4.2.2)
 GetoptLong         1.0.5       2020-12-15 [1] RSPM (R 4.2.2)
 ggalluvial         0.12.5      2023-07-30 [1] Github (corybrunson/ggalluvial@0b71a4d)
 ggbeeswarm         0.7.2       2023-04-29 [1] RSPM (R 4.3.0)
 ggmin              0.0.0.9000  2022-11-24 [1] Github (sjessa/ggmin@8ada274)
 ggnetwork          0.5.12      2023-04-28 [1] Github (briatte/ggnetwork@5b6ed24)
 ggplot2          * 3.4.2       2023-04-03 [1] RSPM (R 4.3.0)
 ggprism            1.0.4       2022-11-24 [1] Github (csdaw/ggprism@0e411f4)
 ggpubr             0.6.0.999   2023-04-28 [1] Github (kassambara/ggpubr@6aeb4f7)
 ggrastr            1.0.2       2023-06-01 [1] RSPM (R 4.3.0)
 ggrepel            0.9.3       2023-02-03 [1] RSPM (R 4.3.0)
 ggridges           0.5.4       2022-09-26 [1] RSPM (R 4.2.2)
 ggsignif           0.6.4.9000  2023-07-30 [1] Github (const-ae/ggsignif@080e119)
 ggstatsplot      * 0.11.1.9000 2023-04-28 [1] Github (IndrajeetPatil/ggstatsplot@befe812)
 GlobalOptions      0.1.2       2020-06-10 [1] RSPM (R 4.2.2)
 globals            0.16.2      2022-11-21 [1] RSPM (R 4.2.2)
 glue               1.6.2       2022-02-24 [1] RSPM (R 4.2.2)
 goftest            1.2-3       2021-10-07 [1] RSPM (R 4.2.2)
 gridBase           0.4-7       2014-02-24 [1] RSPM (R 4.3.0)
 gridExtra          2.3         2017-09-09 [1] RSPM (R 4.2.2)
 gtable             0.3.3       2023-03-21 [1] RSPM (R 4.3.0)
 hdf5r              1.3.8       2023-01-21 [1] RSPM (R 4.3.0)
 here             * 1.0.1       2020-12-13 [1] RSPM (R 4.2.2)
 hms                1.1.3       2023-03-21 [1] RSPM (R 4.3.0)
 htmltools          0.5.5       2023-03-23 [1] RSPM (R 4.3.0)
 htmlwidgets        1.6.2       2023-03-17 [1] RSPM (R 4.3.0)
 httpuv             1.6.11      2023-05-11 [1] RSPM (R 4.3.0)
 httr               1.4.6       2023-05-08 [1] RSPM (R 4.3.0)
 ica                1.0-3       2022-07-08 [1] RSPM (R 4.2.2)
 igraph           * 1.5.0.1     2023-07-23 [1] RSPM (R 4.3.1)
 insight            0.19.2      2023-05-23 [1] RSPM (R 4.3.0)
 IRanges            2.34.0      2023-04-25 [1] Bioconductor
 irlba              2.3.5.1     2022-10-03 [1] RSPM (R 4.2.2)
 iterators          1.0.14      2022-02-05 [1] RSPM (R 4.2.2)
 janitor            2.2.0.9000  2023-04-28 [1] Github (sfirke/janitor@d64c8bb)
 jsonlite           1.8.7       2023-06-29 [1] RSPM (R 4.3.1)
 KernSmooth         2.23-21     2023-05-03 [3] CRAN (R 4.3.1)
 knitr              1.43        2023-05-25 [1] RSPM (R 4.3.0)
 labeling           0.4.2       2020-10-20 [1] RSPM (R 4.2.2)
 later              1.3.1       2023-05-02 [1] RSPM (R 4.3.0)
 lattice            0.21-8      2023-04-05 [1] RSPM (R 4.3.0)
 lazyeval           0.2.2       2019-03-15 [1] RSPM (R 4.2.2)
 leiden             0.4.3       2022-09-10 [1] RSPM (R 4.2.2)
 lifecycle          1.0.3       2022-10-07 [1] RSPM (R 4.2.2)
 listenv            0.9.0       2022-12-16 [1] RSPM (R 4.3.0)
 lmtest             0.9-40      2022-03-21 [1] RSPM (R 4.2.2)
 lubridate        * 1.9.2       2023-02-10 [1] RSPM (R 4.3.0)
 magick             2.7.4       2023-03-09 [1] RSPM (R 4.3.0)
 magrittr         * 2.0.3       2022-03-30 [1] RSPM (R 4.2.2)
 MASS               7.3-60      2023-05-04 [1] RSPM (R 4.3.0)
 Matrix             1.6-0       2023-07-08 [1] RSPM (R 4.3.1)
 matrixStats        1.0.0       2023-06-02 [1] RSPM (R 4.3.0)
 mime               0.12        2021-09-28 [1] RSPM (R 4.2.2)
 miniUI             0.1.1.1     2018-05-18 [1] RSPM (R 4.2.2)
 multcomp           1.4-24      2023-06-14 [1] RSPM (R 4.3.0)
 munsell            0.5.0       2018-06-12 [1] RSPM (R 4.2.2)
 mvtnorm            1.2-2       2023-06-08 [1] RSPM (R 4.3.0)
 network            1.18.1      2023-01-24 [1] RSPM (R 4.3.0)
 nlme               3.1-162     2023-01-31 [3] CRAN (R 4.2.2)
 NMF                0.26        2023-03-20 [1] RSPM (R 4.3.0)
 paletteer          1.5.0       2022-10-19 [1] RSPM (R 4.2.2)
 parallelly         1.36.0      2023-05-26 [1] RSPM (R 4.3.0)
 parameters         0.21.1      2023-05-26 [1] RSPM (R 4.3.0)
 patchwork        * 1.1.2.9000  2023-04-28 [1] Github (thomasp85/patchwork@c14c960)
 pbapply            1.7-2       2023-06-27 [1] RSPM (R 4.3.1)
 pillar             1.9.0       2023-03-22 [1] RSPM (R 4.3.0)
 pkgconfig          2.0.3       2019-09-22 [1] RSPM (R 4.2.2)
 plotly             4.10.2      2023-06-03 [1] RSPM (R 4.3.0)
 plyr               1.8.8       2022-11-11 [1] CRAN (R 4.2.2)
 png                0.1-8       2022-11-29 [1] RSPM (R 4.3.0)
 polyclip           1.10-4      2022-10-20 [1] RSPM (R 4.2.2)
 progressr          0.13.0      2023-01-10 [1] RSPM (R 4.3.0)
 promises           1.2.0.1     2021-02-11 [1] RSPM (R 4.2.2)
 purrr            * 1.0.1       2023-01-10 [1] RSPM (R 4.3.0)
 R.methodsS3        1.8.2       2022-06-13 [1] RSPM (R 4.2.2)
 R.oo               1.25.0      2022-06-12 [1] RSPM (R 4.2.2)
 R.utils            2.12.2      2022-11-11 [1] CRAN (R 4.2.2)
 R6                 2.5.1       2021-08-19 [1] RSPM (R 4.2.2)
 ragg               1.2.5       2023-01-12 [1] RSPM (R 4.3.0)
 RANN               2.6.1       2019-01-08 [1] RSPM (R 4.2.2)
 rappdirs           0.3.3       2021-01-31 [1] RSPM (R 4.2.2)
 RColorBrewer     * 1.1-3       2022-04-03 [1] RSPM (R 4.2.2)
 Rcpp               1.0.11      2023-07-06 [1] RSPM (R 4.3.1)
 RcppAnnoy          0.0.20      2022-10-27 [1] RSPM (R 4.2.2)
 readr            * 2.1.4       2023-02-10 [1] RSPM (R 4.3.0)
 registry           0.5-1       2019-03-05 [1] RSPM (R 4.2.2)
 rematch2           2.1.2       2020-05-01 [1] RSPM (R 4.2.2)
 remotes            2.4.2       2021-11-30 [1] RSPM (R 4.2.2)
 repr               1.1.6       2023-01-26 [1] RSPM (R 4.3.0)
 reshape2           1.4.4       2020-04-09 [1] RSPM (R 4.2.2)
 reticulate       * 1.30        2023-06-09 [1] RSPM (R 4.3.0)
 rjson              0.2.21      2022-01-09 [1] RSPM (R 4.2.2)
 rlang              1.1.1       2023-04-28 [1] RSPM (R 4.3.0)
 rmarkdown          2.22        2023-06-01 [1] RSPM (R 4.3.0)
 rngtools           1.5.2       2021-09-20 [1] RSPM (R 4.3.0)
 ROCR               1.0-11      2020-05-02 [1] RSPM (R 4.2.2)
 rprojroot          2.0.3       2022-04-02 [1] RSPM (R 4.2.2)
 RSpectra           0.16-1      2022-04-24 [1] RSPM (R 4.2.2)
 rstatix            0.7.2       2023-02-01 [1] RSPM (R 4.3.0)
 rstudioapi         0.14        2022-08-22 [1] RSPM (R 4.2.2)
 rsvd               1.0.5       2021-04-16 [1] RSPM (R 4.2.2)
 Rtsne              0.16        2022-04-17 [1] RSPM (R 4.2.2)
 S4Vectors          0.38.1      2023-05-02 [1] Bioconductor
 sandwich           3.0-2       2022-06-15 [1] RSPM (R 4.3.0)
 scales             1.2.1       2022-08-20 [1] RSPM (R 4.2.2)
 scattermore        1.2         2023-06-12 [1] RSPM (R 4.3.0)
 scCustomize      * 1.1.1       2023-04-28 [1] Github (samuel-marsh/scCustomize@d08268d)
 sceasy           * 0.0.7       2022-11-24 [1] Github (cellgeni/sceasy@0cfc0e3)
 sctransform        0.3.5       2022-09-21 [1] RSPM (R 4.2.2)
 sessioninfo        1.2.2       2021-12-06 [1] RSPM (R 4.2.2)
 Seurat           * 4.3.0       2022-11-18 [1] RSPM (R 4.2.2)
 SeuratDisk       * 0.0.0.9020  2022-11-12 [1] Github (mojaveazure/seurat-disk@9b89970)
 SeuratObject     * 4.1.3       2022-11-07 [1] RSPM (R 4.2.2)
 SeuratWrappers   * 0.3.1       2022-11-12 [1] Github (satijalab/seurat-wrappers@d28512f)
 shape              1.4.6       2021-05-19 [1] RSPM (R 4.2.2)
 shiny              1.7.4       2022-12-15 [1] RSPM (R 4.3.0)
 skimr            * 2.1.5       2023-04-28 [1] Github (ropensci/skimr@d5126aa)
 sna                2.7-1       2023-01-24 [1] RSPM (R 4.3.0)
 snakecase          0.11.0      2019-05-25 [1] RSPM (R 4.2.2)
 sp                 1.6-1       2023-05-31 [1] RSPM (R 4.3.0)
 spatstat.data      3.0-1       2023-03-12 [1] RSPM (R 4.3.0)
 spatstat.explore   3.2-1       2023-05-13 [1] RSPM (R 4.3.0)
 spatstat.geom      3.2-1       2023-05-09 [1] RSPM (R 4.3.0)
 spatstat.random    3.1-5       2023-05-11 [1] RSPM (R 4.3.0)
 spatstat.sparse    3.0-1       2023-03-12 [1] RSPM (R 4.3.0)
 spatstat.utils     3.0-3       2023-05-09 [1] RSPM (R 4.3.0)
 statnet.common     4.9.0       2023-05-24 [1] RSPM (R 4.3.0)
 statsExpressions   1.5.1       2023-05-23 [1] RSPM (R 4.3.0)
 stringi            1.7.12      2023-01-11 [1] RSPM (R 4.3.0)
 stringr          * 1.5.0       2022-12-02 [1] RSPM (R 4.3.0)
 survival           3.5-5       2023-03-12 [1] RSPM (R 4.3.0)
 svglite            2.1.1       2023-01-10 [1] RSPM (R 4.3.0)
 systemfonts        1.0.4       2022-02-11 [1] RSPM (R 4.2.2)
 tensor             1.5         2012-05-05 [1] RSPM (R 4.2.2)
 textshaping        0.3.6       2021-10-13 [1] RSPM (R 4.2.2)
 TH.data            1.1-2       2023-04-17 [1] RSPM (R 4.3.0)
 tibble           * 3.2.1       2023-03-20 [1] RSPM (R 4.3.0)
 tidyr            * 1.3.0       2023-01-24 [1] RSPM (R 4.3.0)
 tidyselect         1.2.0       2022-10-10 [1] RSPM (R 4.2.2)
 tidyverse        * 2.0.0.9000  2023-04-28 [1] Github (tidyverse/tidyverse@8ec2e1f)
 timechange         0.2.0       2023-01-11 [1] RSPM (R 4.3.0)
 tzdb               0.4.0       2023-05-12 [1] RSPM (R 4.3.0)
 utf8               1.2.3       2023-01-31 [1] RSPM (R 4.3.0)
 uwot               0.1.14      2022-08-22 [1] RSPM (R 4.2.2)
 vctrs              0.6.3       2023-06-14 [1] RSPM (R 4.3.0)
 vipor              0.4.5       2017-03-22 [1] RSPM (R 4.2.2)
 viridis          * 0.6.3       2023-05-03 [1] RSPM (R 4.3.0)
 viridisLite      * 0.4.2       2023-05-02 [1] RSPM (R 4.3.0)
 vroom              1.6.3       2023-04-28 [1] RSPM (R 4.3.0)
 withr              2.5.0       2022-03-03 [1] RSPM (R 4.2.2)
 wordcloud          2.6         2018-08-24 [1] RSPM (R 4.3.0)
 xfun               0.39        2023-04-20 [1] RSPM (R 4.3.0)
 xtable             1.8-4       2019-04-21 [1] RSPM (R 4.2.2)
 yaml               2.3.7       2023-01-23 [1] RSPM (R 4.3.0)
 zeallot            0.1.0       2018-01-28 [1] RSPM (R 4.2.2)
 zoo                1.8-12      2023-04-13 [1] RSPM (R 4.3.0)

 [1] /usr/local/lib/R/site-library
 [2] /usr/lib/R/site-library
 [3] /usr/lib/R/library

─ Python configuration ───────────────────────────────────────────────────────
 python:         /home/etretiakov/micromamba/bin/python3
 libpython:      /home/etretiakov/micromamba/lib/libpython3.9.so
 pythonhome:     /home/etretiakov/micromamba:/home/etretiakov/micromamba
 version:        3.9.16 (main, Jan 11 2023, 16:05:54)  [GCC 11.2.0]
 numpy:          /home/etretiakov/micromamba/lib/python3.9/site-packages/numpy
 numpy_version:  1.23.5
 scanpy:         /home/etretiakov/micromamba/lib/python3.9/site-packages/scanpy
 
 python versions found: 
  /home/etretiakov/micromamba/bin/python3
  /home/etretiakov/micromamba/bin/python

──────────────────────────────────────────────────────────────────────────────