Script 17a: CellChat — Full Niche Cell-Cell Communication
B-ALL CNS Niche Project
2026-03-18
Overview
CellChat analysis of the full B-ALL niche using the unified
annotation object (15c_immune_unified_annotated.qs). Runs
CNS and BM separately then performs comparative analysis.
Unified annotation populations (36 total, excluding artefacts): - Myeloid: cpepiMF, BAM_1, BAM_2, BAM_3, BM Macrophages, BM Macrophages (activated), BM Monocyte-derived Macrophages, dmMF, pvMF, Kolmer cell, Proliferating Myeloid, Myeloid DCs, Non-classical monocytes - Stromal: CNS Meningeal Fibroblasts (ECM), CNS Leptomeningeal Fibroblasts, CNS Mesenchymal/Adventitial, BM Stroma - T cells: CD8_Naive, CD8_EffectorMemory, CD8_Tex, CD4_CTL, Th1, Treg, CD4_other, Vγ6Vδ4, Vγ4, Other γδ, NKT-like subsets, CD8NKT_NKT1 - B cells/other: Pre-B cells, Plasma B cells, NK cells
Key biological questions: 1. What signals does the CNS niche send to CD8 T cells — particularly the retention/suppression programme (Itga1, Socs1, Cish, Bcl2)? 2. Do Vγ6Vδ4 γδ T cells interact with CD8 T cells in CNS? 3. What distinguishes CNS vs BM interaction networks? 4. Which populations are the dominant senders in each niche?
Reference: Jin et al. 2021 Nature Communications — CellChat
1. Libraries
library(CellChat)
library(Seurat)
library(ggplot2)
library(dplyr)
library(patchwork)
library(qs2)
library(tidyr)
library(RColorBrewer)
BASE_DIR <- "/Volumes/Al_seq/scratch_backup/harmony_clustering/"
# Colour palette for populations — will be built after loading
tissue_cols <- c("BM" = "#b2182b", "CNS" = "#2166ac")2. Load data and prepare
seu <- qs_read(paste0(BASE_DIR, "15c_immune_unified_annotated.qs"))
DefaultAssay(seu) <- "originalexp"
# Remove excluded populations and set identity
seu_clean <- subset(seu,
cells = colnames(seu)[
seu$annotation_unified != "exclude" &
!is.na(seu$annotation_unified)
])
Idents(seu_clean) <- "annotation_unified"
cat("=== Clean object ===\n")## === Clean object ===cat("Total cells:", ncol(seu_clean), "\n")## Total cells: 15145cat("Populations:", length(levels(Idents(seu_clean))), "\n\n")## Populations: 36print(table(seu_clean$annotation_unified, seu_clean$Tissue))##
## BM CNS
## BAM_1 (Marco+/activated) 21 544
## BAM_2 (IFN-responsive) 0 115
## BAM_3 (transitional) 0 73
## BM Macrophages 499 83
## BM Macrophages (activated) 406 27
## BM Monocyte-derived Macrophages 60 3
## BM Stroma 719 0
## CD4_CTL 162 84
## CD4_other 40 11
## CD8_EffectorMemory 307 183
## CD8_Naive 1610 127
## CD8_Tex 99 15
## CD8NKT_NKT1 156 100
## CD8NKT_other 17 15
## CNS Leptomeningeal Fibroblasts 0 326
## CNS Meningeal Fibroblasts (ECM) 0 275
## CNS Mesenchymal/Adventitial 0 122
## cpepiMF 64 1169
## dmMF (MHC-II+) 5 123
## Kolmer cell 16 10
## Myeloid DCs 76 41
## NK cells 887 92
## NKT-like_NKT1 62 28
## NKT-like_NKT10 28 65
## NKT-like_NKT17 89 30
## NKT-like_NKT2 40 46
## Non-classical monocytes 356 250
## Other γδ 253 252
## Plasma B cells 754 90
## Pre-B cells 2893 205
## Proliferating Myeloid 120 82
## pvMF 0 108
## Th1 272 78
## Treg 157 10
## Vγ4 40 32
## Vγ6Vδ4 2 121# Split into CNS and BM objects
seu_cns <- subset(seu_clean, subset = Tissue == "CNS")
seu_bm <- subset(seu_clean, subset = Tissue == "BM")
# Remove populations with <10 cells in each tissue
min_cells <- 10
cns_keep <- names(table(seu_cns$annotation_unified))[
table(seu_cns$annotation_unified) >= min_cells
]
bm_keep <- names(table(seu_bm$annotation_unified))[
table(seu_bm$annotation_unified) >= min_cells
]
seu_cns <- subset(seu_cns,
cells = colnames(seu_cns)[
seu_cns$annotation_unified %in% cns_keep
])
seu_bm <- subset(seu_bm,
cells = colnames(seu_bm)[
seu_bm$annotation_unified %in% bm_keep
])
cat("CNS populations (≥10 cells):", length(cns_keep), "\n")## CNS populations (≥10 cells): 34print(sort(table(seu_cns$annotation_unified), decreasing = TRUE))##
## cpepiMF BAM_1 (Marco+/activated)
## 1169 544
## CNS Leptomeningeal Fibroblasts CNS Meningeal Fibroblasts (ECM)
## 326 275
## Other γδ Non-classical monocytes
## 252 250
## Pre-B cells CD8_EffectorMemory
## 205 183
## CD8_Naive dmMF (MHC-II+)
## 127 123
## CNS Mesenchymal/Adventitial Vγ6Vδ4
## 122 121
## BAM_2 (IFN-responsive) pvMF
## 115 108
## CD8NKT_NKT1 NK cells
## 100 92
## Plasma B cells CD4_CTL
## 90 84
## BM Macrophages Proliferating Myeloid
## 83 82
## Th1 BAM_3 (transitional)
## 78 73
## NKT-like_NKT10 NKT-like_NKT2
## 65 46
## Myeloid DCs Vγ4
## 41 32
## NKT-like_NKT17 NKT-like_NKT1
## 30 28
## BM Macrophages (activated) CD8_Tex
## 27 15
## CD8NKT_other CD4_other
## 15 11
## Kolmer cell Treg
## 10 10cat("\nBM populations (≥10 cells):", length(bm_keep), "\n")##
## BM populations (≥10 cells): 28print(sort(table(seu_bm$annotation_unified), decreasing = TRUE))##
## Pre-B cells CD8_Naive
## 2893 1610
## NK cells Plasma B cells
## 887 754
## BM Stroma BM Macrophages
## 719 499
## BM Macrophages (activated) Non-classical monocytes
## 406 356
## CD8_EffectorMemory Th1
## 307 272
## Other γδ CD4_CTL
## 253 162
## Treg CD8NKT_NKT1
## 157 156
## Proliferating Myeloid CD8_Tex
## 120 99
## NKT-like_NKT17 Myeloid DCs
## 89 76
## cpepiMF NKT-like_NKT1
## 64 62
## BM Monocyte-derived Macrophages CD4_other
## 60 40
## NKT-like_NKT2 Vγ4
## 40 40
## NKT-like_NKT10 BAM_1 (Marco+/activated)
## 28 21
## CD8NKT_other Kolmer cell
## 17 163. Create CellChat objects
# Helper function to create CellChat object from Seurat
make_cellchat <- function(srt, label = "object") {
cat("Creating CellChat object for:", label, "\n")
# Extract normalised counts and metadata
data_input <- GetAssayData(srt, layer = "data")
meta <- srt@meta.data %>%
select(annotation_unified) %>%
rename(labels = annotation_unified)
meta$labels <- factor(meta$labels)
# Create CellChat object
cc <- createCellChat(
object = data_input,
meta = meta,
group.by = "labels"
)
# Set ligand-receptor database — mouse
CellChatDB <- CellChatDB.mouse
# Use all interactions (Secreted + Cell-Cell Contact + ECM-Receptor)
CellChatDB.use <- CellChatDB
cc@DB <- CellChatDB.use
cat(" Cells:", ncol(data_input), "\n")
cat(" Populations:", nlevels(meta$labels), "\n")
cat(" Database interactions:", nrow(CellChatDB.use$interaction), "\n")
cc
}
cc_cns <- make_cellchat(seu_cns, "CNS")## Creating CellChat object for: CNS
## [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 BAM_1 (Marco+/activated), BAM_2 (IFN-responsive), BAM_3 (transitional), BM Macrophages, BM Macrophages (activated), CD4_CTL, CD4_other, CD8_EffectorMemory, CD8_Naive, CD8_Tex, CD8NKT_NKT1, CD8NKT_other, CNS Leptomeningeal Fibroblasts, CNS Meningeal Fibroblasts (ECM), CNS Mesenchymal/Adventitial, cpepiMF, dmMF (MHC-II+), Kolmer cell, Myeloid DCs, NK cells, NKT-like_NKT1, NKT-like_NKT10, NKT-like_NKT17, NKT-like_NKT2, Non-classical monocytes, Other γδ, Plasma B cells, Pre-B cells, Proliferating Myeloid, pvMF, Th1, Treg, Vγ4, Vγ6Vδ4
## Cells: 4932
## Populations: 34
## Database interactions: 3379cc_bm <- make_cellchat(seu_bm, "BM")## Creating CellChat object for: BM
## [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 BAM_1 (Marco+/activated), BM Macrophages, BM Macrophages (activated), BM Monocyte-derived Macrophages, BM Stroma, CD4_CTL, CD4_other, CD8_EffectorMemory, CD8_Naive, CD8_Tex, CD8NKT_NKT1, CD8NKT_other, cpepiMF, Kolmer cell, Myeloid DCs, NK cells, NKT-like_NKT1, NKT-like_NKT10, NKT-like_NKT17, NKT-like_NKT2, Non-classical monocytes, Other γδ, Plasma B cells, Pre-B cells, Proliferating Myeloid, Th1, Treg, Vγ4
## Cells: 10203
## Populations: 28
## Database interactions: 33794. Run CellChat pipeline
cat("=== Running CellChat: CNS ===\n")## === Running CellChat: CNS ===# Pre-process
cc_cns <- subsetData(cc_cns)
cc_cns <- identifyOverExpressedGenes(cc_cns)
cc_cns <- identifyOverExpressedInteractions(cc_cns)## The number of highly variable ligand-receptor pairs used for signaling inference is 2554# Infer communication network
cc_cns <- computeCommunProb(
cc_cns,
type = "triMean",
trim = 0.1,
nboot = 100
)## triMean is used for calculating the average gene expression per cell group.
## [1] ">>> Run CellChat on sc/snRNA-seq data <<< [2026-03-18 08:36:04.130888]"
## [1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2026-03-18 08:39:42.431555]"# Filter low-probability interactions
cc_cns <- filterCommunication(cc_cns, min.cells = 10)## The cell-cell communication related with the following cell groups are excluded due to the few number of cells: Kolmer cell, Treg ! 15.6% interactions are removed!# Pathway-level communication
cc_cns <- computeCommunProbPathway(cc_cns)
# Aggregate network
cc_cns <- aggregateNet(cc_cns)
cat("CNS interactions computed.\n")## CNS interactions computed.cat("Significant LR pairs:",
nrow(subsetCommunication(cc_cns)), "\n")## Significant LR pairs: 14701cat("=== Running CellChat: BM ===\n")## === Running CellChat: BM ===cc_bm <- subsetData(cc_bm)
cc_bm <- identifyOverExpressedGenes(cc_bm)
cc_bm <- identifyOverExpressedInteractions(cc_bm)## The number of highly variable ligand-receptor pairs used for signaling inference is 2473cc_bm <- computeCommunProb(
cc_bm,
type = "triMean",
trim = 0.1,
nboot = 100
)## triMean is used for calculating the average gene expression per cell group.
## [1] ">>> Run CellChat on sc/snRNA-seq data <<< [2026-03-18 08:39:47.319399]"
## [1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2026-03-18 08:43:15.315278]"cc_bm <- filterCommunication(cc_bm, min.cells = 10)
cc_bm <- computeCommunProbPathway(cc_bm)
cc_bm <- aggregateNet(cc_bm)
cat("BM interactions computed.\n")## BM interactions computed.cat("Significant LR pairs:",
nrow(subsetCommunication(cc_bm)), "\n")## Significant LR pairs: 102905. Overview figures
Fig 1: Interaction number and strength — CNS vs BM
par(mfrow = c(1, 2))
netVisual_circle(
cc_cns@net$count,
vertex.weight = as.numeric(table(cc_cns@idents)),
weight.scale = TRUE,
label.edge = FALSE,
title.name = "Interaction count — CNS"
)
netVisual_circle(
cc_bm@net$count,
vertex.weight = as.numeric(table(cc_bm@idents)),
weight.scale = TRUE,
label.edge = FALSE,
title.name = "Interaction count — BM"
)
Fig 2: Dominant senders and receivers — CNS
# Compute network centrality
cc_cns <- netAnalysis_computeCentrality(cc_cns, slot.name = "netP")
cc_bm <- netAnalysis_computeCentrality(cc_bm, slot.name = "netP")
netAnalysis_signalingRole_scatter(cc_cns) +
labs(title = "Signalling roles — CNS",
subtitle = "X = outgoing strength (sender); Y = incoming strength (receiver)")
netAnalysis_signalingRole_scatter(cc_bm) +
labs(title = "Signalling roles — BM",
subtitle = "X = outgoing strength (sender); Y = incoming strength (receiver)")
Fig 3: Heatmap of outgoing signals — CNS
netAnalysis_signalingRole_heatmap(
cc_cns,
pattern = "outgoing",
signaling = NULL,
title = "Outgoing signalling — CNS",
width = 12,
height = 36
)
Fig 4: Heatmap of incoming signals — CNS
netAnalysis_signalingRole_heatmap(
cc_cns,
pattern = "incoming",
signaling = NULL,
title = "Incoming signalling — CNS",
width = 12,
height = 36
)
6. Targeted analysis — signals to CD8 T cells
# What signals does the CNS niche send TO CD8 T cells?
cd8_pops <- intersect(
c("CD8_EffectorMemory", "CD8_Tex", "CD8_Naive"),
levels(cc_cns@idents)
)
cat("CD8 populations in CNS CellChat:\n")## CD8 populations in CNS CellChat:print(cd8_pops)## [1] "CD8_EffectorMemory" "CD8_Tex" "CD8_Naive"# Extract all interactions where receiver = CD8 population
cd8_incoming_cns <- subsetCommunication(
cc_cns,
targets.use = cd8_pops
)
cat("\n=== Top pathways signalling TO CD8 T cells in CNS ===\n")##
## === Top pathways signalling TO CD8 T cells in CNS ===cd8_incoming_cns %>%
group_by(pathway_name) %>%
summarise(
n_interactions = n(),
mean_prob = round(mean(prob), 4),
senders = paste(unique(source), collapse = ", "),
.groups = "drop"
) %>%
arrange(desc(mean_prob)) %>%
head(20) %>%
print()## # A tibble: 20 × 4
## pathway_name n_interactions mean_prob senders
## <chr> <int> <dbl> <chr>
## 1 SELPLG 26 0.106 BAM_2 (IFN-responsive), BAM_3 (transit…
## 2 PTN 9 0.0955 CNS Leptomeningeal Fibroblasts, CNS Me…
## 3 MK 15 0.0839 CNS Leptomeningeal Fibroblasts, CNS Me…
## 4 LCK 26 0.0769 CD4_CTL, CD4_other, CD8_EffectorMemory…
## 5 MHC-I 960 0.0712 BAM_2 (IFN-responsive), CD4_other, CD8…
## 6 CCL 29 0.0648 CD8_EffectorMemory, CD8_Tex, CD4_CTL, …
## 7 LAIR1 6 0.0595 BAM_2 (IFN-responsive), BAM_3 (transit…
## 8 FN1 18 0.0531 CNS Leptomeningeal Fibroblasts, CNS Me…
## 9 GALECTIN 111 0.0475 BAM_1 (Marco+/activated), BAM_2 (IFN-r…
## 10 CEACAM 6 0.0423 CNS Mesenchymal/Adventitial, NK cells,…
## 11 JAM 11 0.0409 BAM_2 (IFN-responsive), BAM_3 (transit…
## 12 CXCL 33 0.0397 BAM_1 (Marco+/activated), BM Macrophag…
## 13 ADGRE 91 0.0342 CD4_CTL, CD4_other, CD8_EffectorMemory…
## 14 SIRP 12 0.0337 Non-classical monocytes, Myeloid DCs
## 15 THBS 27 0.0267 CNS Leptomeningeal Fibroblasts, CNS Me…
## 16 MHC-II 102 0.0239 BM Macrophages, CNS Leptomeningeal Fib…
## 17 COLLAGEN 40 0.0181 CNS Leptomeningeal Fibroblasts, CNS Me…
## 18 APP 21 0.0174 BAM_1 (Marco+/activated), BAM_2 (IFN-r…
## 19 PD-L1 4 0.017 Non-classical monocytes, Plasma B cells
## 20 L1CAM 4 0.0164 Myeloid DCs, Non-classical monocytes# Bubble plot: all senders → CD8 T cells in CNS
if (length(cd8_pops) > 0) {
netVisual_bubble(
cc_cns,
targets.use = cd8_pops,
remove.isolate = TRUE,
angle.x = 45
) +
labs(title = "Signals received by CD8 T cells — CNS",
subtitle = "All sender populations → CD8_EffectorMemory / CD8_Tex / CD8_Naive")
}
7. Targeted analysis — Vγ6Vδ4 as sender
vg6_present_cns <- "Vγ6Vδ4" %in% levels(cc_cns@idents)
vg6_present_bm <- "Vγ6Vδ4" %in% levels(cc_bm@idents)
cat("Vγ6Vδ4 in CNS CellChat:", vg6_present_cns, "\n")## Vγ6Vδ4 in CNS CellChat: TRUEcat("Vγ6Vδ4 in BM CellChat:", vg6_present_bm, "\n")## Vγ6Vδ4 in BM CellChat: FALSEif (vg6_present_cns) {
vg6_outgoing <- subsetCommunication(
cc_cns,
sources.use = "Vγ6Vδ4"
)
cat("\n=== Vγ6Vδ4 outgoing interactions in CNS ===\n")
vg6_outgoing %>%
arrange(desc(prob)) %>%
select(source, target, ligand, receptor, pathway_name, prob) %>%
head(30) %>%
print()
}##
## === Vγ6Vδ4 outgoing interactions in CNS ===
## source target ligand receptor pathway_name
## 1 Vγ6Vδ4 CNS Leptomeningeal Fibroblasts Ppia Bsg CypA
## 2 Vγ6Vδ4 dmMF (MHC-II+) Ptprc Mrc1 CD45
## 3 Vγ6Vδ4 CNS Meningeal Fibroblasts (ECM) Ppia Bsg CypA
## 4 Vγ6Vδ4 CNS Mesenchymal/Adventitial Ppia Bsg CypA
## 5 Vγ6Vδ4 CD8_EffectorMemory H2-K1 Cd8b1 MHC-I
## 6 Vγ6Vδ4 CD8_EffectorMemory H2-D1 Cd8b1 MHC-I
## 7 Vγ6Vδ4 pvMF Ptprc Mrc1 CD45
## 8 Vγ6Vδ4 CD8_Tex H2-K1 Cd8a MHC-I
## 9 Vγ6Vδ4 CD8_Tex H2-D1 Cd8a MHC-I
## 10 Vγ6Vδ4 Non-classical monocytes Thy1 Adgre5 ADGRE
## 11 Vγ6Vδ4 CD8_EffectorMemory H2-K1 Cd8a MHC-I
## 12 Vγ6Vδ4 CD8_EffectorMemory H2-D1 Cd8a MHC-I
## 13 Vγ6Vδ4 BAM_2 (IFN-responsive) Ppia Bsg CypA
## 14 Vγ6Vδ4 cpepiMF Ppia Bsg CypA
## 15 Vγ6Vδ4 CD8_Tex H2-K1 Cd8b1 MHC-I
## 16 Vγ6Vδ4 CD8_Tex H2-D1 Cd8b1 MHC-I
## 17 Vγ6Vδ4 NK cells Selplg Sell SELPLG
## 18 Vγ6Vδ4 Non-classical monocytes Thy1 ITGAM_ITGB2 THY1
## 19 Vγ6Vδ4 CD8_Naive Selplg Sell SELPLG
## 20 Vγ6Vδ4 BM Macrophages (activated) Ptprc Mrc1 CD45
## 21 Vγ6Vδ4 Pre-B cells Cd52 Siglecg CD52
## 22 Vγ6Vδ4 dmMF (MHC-II+) Ppia Bsg CypA
## 23 Vγ6Vδ4 CD8NKT_other Clec2d Klrb1b CLEC
## 24 Vγ6Vδ4 CD8NKT_other Thy1 Adgre5 ADGRE
## 25 Vγ6Vδ4 BAM_2 (IFN-responsive) Thy1 ITGAM_ITGB2 THY1
## 26 Vγ6Vδ4 CD8NKT_other Clec2d Klrb1c CLEC
## 27 Vγ6Vδ4 cpepiMF Thy1 ITGAM_ITGB2 THY1
## 28 Vγ6Vδ4 Non-classical monocytes Thy1 ITGAX_ITGB2 THY1
## 29 Vγ6Vδ4 Pre-B cells Selplg Sell SELPLG
## 30 Vγ6Vδ4 Myeloid DCs Thy1 ITGAM_ITGB2 THY1
## prob
## 1 0.23143812
## 2 0.22231209
## 3 0.21968154
## 4 0.21883499
## 5 0.20845291
## 6 0.20709246
## 7 0.20218566
## 8 0.18563780
## 9 0.18439157
## 10 0.17299676
## 11 0.17298234
## 12 0.17180315
## 13 0.16099346
## 14 0.15900165
## 15 0.15516272
## 16 0.15408237
## 17 0.15035416
## 18 0.14192165
## 19 0.13581211
## 20 0.13495697
## 21 0.12035322
## 22 0.11903515
## 23 0.11787538
## 24 0.11587252
## 25 0.11353825
## 26 0.10582510
## 27 0.10452176
## 28 0.10342231
## 29 0.09776311
## 30 0.09559582if (vg6_present_cns) {
netVisual_bubble(
cc_cns,
sources.use = "Vγ6Vδ4",
remove.isolate = TRUE,
angle.x = 45
) +
labs(title = "Vγ6Vδ4 outgoing signals — CNS",
subtitle = "All targets")
}
8. Targeted analysis — stromal senders to T cells
stroma_pops <- intersect(
c("CNS Meningeal Fibroblasts (ECM)",
"CNS Leptomeningeal Fibroblasts",
"CNS Mesenchymal/Adventitial"),
levels(cc_cns@idents)
)
t_cell_pops <- intersect(
c("CD8_EffectorMemory", "CD8_Tex", "CD8_Naive",
"Th1", "Treg", "CD4_CTL", "Vγ6Vδ4"),
levels(cc_cns@idents)
)
cat("Stromal populations in CNS:\n")## Stromal populations in CNS:print(stroma_pops)## [1] "CNS Meningeal Fibroblasts (ECM)" "CNS Leptomeningeal Fibroblasts"
## [3] "CNS Mesenchymal/Adventitial"cat("\nT cell populations in CNS:\n")##
## T cell populations in CNS:print(t_cell_pops)## [1] "CD8_EffectorMemory" "CD8_Tex" "CD8_Naive"
## [4] "Th1" "Treg" "CD4_CTL"
## [7] "Vγ6Vδ4"stroma_to_tcell <- subsetCommunication(
cc_cns,
sources.use = stroma_pops,
targets.use = t_cell_pops
)
cat("\n=== Top stroma → T cell pathways in CNS ===\n")##
## === Top stroma → T cell pathways in CNS ===stroma_to_tcell %>%
group_by(pathway_name, source) %>%
summarise(
n = n(),
mean_prob = round(mean(prob), 4),
.groups = "drop"
) %>%
arrange(desc(mean_prob)) %>%
head(20) %>%
print()## # A tibble: 20 × 4
## pathway_name source n mean_prob
## <chr> <fct> <int> <dbl>
## 1 PTN CNS Leptomeningeal Fibroblasts 9 0.0825
## 2 MK CNS Leptomeningeal Fibroblasts 12 0.0813
## 3 MK CNS Mesenchymal/Adventitial 12 0.0812
## 4 MK CNS Meningeal Fibroblasts (ECM) 12 0.0751
## 5 PTN CNS Mesenchymal/Adventitial 9 0.0746
## 6 PTN CNS Meningeal Fibroblasts (ECM) 9 0.0629
## 7 FN1 CNS Meningeal Fibroblasts (ECM) 17 0.061
## 8 FN1 CNS Mesenchymal/Adventitial 17 0.0529
## 9 MHC-I CNS Leptomeningeal Fibroblasts 20 0.0506
## 10 FN1 CNS Leptomeningeal Fibroblasts 17 0.0504
## 11 GALECTIN CNS Meningeal Fibroblasts (ECM) 6 0.0392
## 12 MHC-I CNS Mesenchymal/Adventitial 30 0.0389
## 13 MHC-I CNS Meningeal Fibroblasts (ECM) 30 0.0385
## 14 GALECTIN CNS Mesenchymal/Adventitial 5 0.0372
## 15 APP CNS Meningeal Fibroblasts (ECM) 5 0.0321
## 16 CXCL CNS Meningeal Fibroblasts (ECM) 7 0.0305
## 17 APP CNS Leptomeningeal Fibroblasts 5 0.0298
## 18 APP CNS Mesenchymal/Adventitial 5 0.0283
## 19 CCL CNS Meningeal Fibroblasts (ECM) 1 0.0282
## 20 CXCL CNS Mesenchymal/Adventitial 7 0.0277if (length(stroma_pops) > 0 & length(t_cell_pops) > 0) {
netVisual_bubble(
cc_cns,
sources.use = stroma_pops,
targets.use = t_cell_pops,
remove.isolate = TRUE,
angle.x = 45
) +
labs(title = "Stromal → T cell signals — CNS",
subtitle = "ECM fibroblasts / Leptomeningeal / Mesenchymal → CD8 + CD4 + γδ")
}
9. Comparative analysis — CNS vs BM
# Find shared populations between CNS and BM
shared_pops <- intersect(
levels(cc_cns@idents),
levels(cc_bm@idents)
)
cat("Shared populations (", length(shared_pops), "):\n")## Shared populations ( 26 ):print(shared_pops)## [1] "BAM_1 (Marco+/activated)" "BM Macrophages"
## [3] "BM Macrophages (activated)" "CD4_CTL"
## [5] "CD4_other" "CD8_EffectorMemory"
## [7] "CD8_Naive" "CD8_Tex"
## [9] "CD8NKT_NKT1" "CD8NKT_other"
## [11] "cpepiMF" "Kolmer cell"
## [13] "Myeloid DCs" "NK cells"
## [15] "NKT-like_NKT1" "NKT-like_NKT10"
## [17] "NKT-like_NKT17" "NKT-like_NKT2"
## [19] "Non-classical monocytes" "Other γδ"
## [21] "Plasma B cells" "Pre-B cells"
## [23] "Proliferating Myeloid" "Th1"
## [25] "Treg" "Vγ4"cat("\nCNS-only (used for targeted analyses in sections 6-8):\n")##
## CNS-only (used for targeted analyses in sections 6-8):print(setdiff(levels(cc_cns@idents), levels(cc_bm@idents)))## [1] "BAM_2 (IFN-responsive)" "BAM_3 (transitional)"
## [3] "CNS Leptomeningeal Fibroblasts" "CNS Meningeal Fibroblasts (ECM)"
## [5] "CNS Mesenchymal/Adventitial" "dmMF (MHC-II+)"
## [7] "pvMF" "Vγ6Vδ4"cat("\nBM-only:\n")##
## BM-only:print(setdiff(levels(cc_bm@idents), levels(cc_cns@idents)))## [1] "BM Monocyte-derived Macrophages" "BM Stroma"# Subset both Seurat objects to shared populations and rerun CellChat
# Note: full CNS object (cc_cns) retained for targeted sender analyses
# in sections 6-8 which require CNS-specific populations (stroma, Vγ6Vδ4 etc.)
seu_cns_shared <- subset(seu_cns,
cells = colnames(seu_cns)[
seu_cns$annotation_unified %in% shared_pops])
seu_bm_shared <- subset(seu_bm,
cells = colnames(seu_bm)[
seu_bm$annotation_unified %in% shared_pops])
cat("CNS shared cells:", ncol(seu_cns_shared), "\n")## CNS shared cells: 3669cat("BM shared cells:", ncol(seu_bm_shared), "\n")## BM shared cells: 9424# Recreate CellChat objects on shared populations
cc_cns_shared <- make_cellchat(seu_cns_shared, "CNS shared")## Creating CellChat object for: CNS shared
## [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 BAM_1 (Marco+/activated), BM Macrophages, BM Macrophages (activated), CD4_CTL, CD4_other, CD8_EffectorMemory, CD8_Naive, CD8_Tex, CD8NKT_NKT1, CD8NKT_other, cpepiMF, Kolmer cell, Myeloid DCs, NK cells, NKT-like_NKT1, NKT-like_NKT10, NKT-like_NKT17, NKT-like_NKT2, Non-classical monocytes, Other γδ, Plasma B cells, Pre-B cells, Proliferating Myeloid, Th1, Treg, Vγ4
## Cells: 3669
## Populations: 26
## Database interactions: 3379cc_bm_shared <- make_cellchat(seu_bm_shared, "BM shared")## Creating CellChat object for: BM shared
## [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 BAM_1 (Marco+/activated), BM Macrophages, BM Macrophages (activated), CD4_CTL, CD4_other, CD8_EffectorMemory, CD8_Naive, CD8_Tex, CD8NKT_NKT1, CD8NKT_other, cpepiMF, Kolmer cell, Myeloid DCs, NK cells, NKT-like_NKT1, NKT-like_NKT10, NKT-like_NKT17, NKT-like_NKT2, Non-classical monocytes, Other γδ, Plasma B cells, Pre-B cells, Proliferating Myeloid, Th1, Treg, Vγ4
## Cells: 9424
## Populations: 26
## Database interactions: 3379# Helper to run full pipeline
run_cellchat <- function(cc, label) {
cat("Running:", label, "\n")
cc <- subsetData(cc)
cc <- identifyOverExpressedGenes(cc)
cc <- identifyOverExpressedInteractions(cc)
cc <- computeCommunProb(cc, type = "triMean",
trim = 0.1, nboot = 20)
cc <- filterCommunication(cc, min.cells = 10)
cc <- computeCommunProbPathway(cc)
cc <- aggregateNet(cc)
cc <- netAnalysis_computeCentrality(cc, slot.name = "netP")
cat(" Complete —", nrow(subsetCommunication(cc)),
"interactions\n")
cc
}
cc_cns_shared <- run_cellchat(cc_cns_shared, "CNS shared")## Running: CNS shared
## The number of highly variable ligand-receptor pairs used for signaling inference is 2108
## triMean is used for calculating the average gene expression per cell group.
## [1] ">>> Run CellChat on sc/snRNA-seq data <<< [2026-03-18 08:43:41.24638]"
## [1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2026-03-18 08:44:13.006833]"
## The cell-cell communication related with the following cell groups are excluded due to the few number of cells: Kolmer cell, Treg ! 21.7% interactions are removed!
## Complete — 6081 interactionscc_bm_shared <- run_cellchat(cc_bm_shared, "BM shared")## Running: BM shared
## The number of highly variable ligand-receptor pairs used for signaling inference is 2083
## triMean is used for calculating the average gene expression per cell group.
## [1] ">>> Run CellChat on sc/snRNA-seq data <<< [2026-03-18 08:44:17.901729]"
## [1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2026-03-18 08:44:55.524508]"
## Complete — 8042 interactions# Merge for comparison
object_list <- list(BM = cc_bm_shared, CNS = cc_cns_shared)
cc_merged <- mergeCellChat(object_list,
add.names = names(object_list))
cat("Merged successfully.\n")## Merged successfully.cat("BM shared interactions:",
nrow(subsetCommunication(cc_bm_shared)), "\n")## BM shared interactions: 8042cat("CNS shared interactions:",
nrow(subsetCommunication(cc_cns_shared)), "\n")## CNS shared interactions: 6081Fig 8: Differential interaction number — CNS vs BM
compareInteractions(
cc_merged,
show.legend = FALSE,
group = c(1, 2),
color.use = c("BM" = "#b2182b", "CNS" = "#2166ac")
) +
labs(title = "Total interactions: BM vs CNS")
Fig 9: Differential outgoing signals — CNS vs BM
# Heatmaps — use individual shared objects (centrality already computed)
netAnalysis_signalingRole_heatmap(
cc_cns_shared,
pattern = "outgoing",
signaling = NULL,
title = "Outgoing signalling — CNS (shared populations)",
width = 12,
height = 8,
color.heatmap = "OrRd"
)
netAnalysis_signalingRole_heatmap(
cc_bm_shared,
pattern = "outgoing",
signaling = NULL,
title = "Outgoing signalling — BM (shared populations)",
width = 12,
height = 8,
color.heatmap = "OrRd"
)
# Pathway ranking — uses merged object, works fine
rankNet(
cc_merged,
mode = "comparison",
stacked = TRUE,
do.stat = TRUE,
color.use = c("BM" = "#b2182b", "CNS" = "#2166ac")
)
# Differential interaction count and weight — uses merged object
compareInteractions(
cc_merged,
show.legend = FALSE,
group = c(1, 2),
color.use = c("BM" = "#b2182b", "CNS" = "#2166ac")
)
Fig 10: CNS-enriched vs BM-enriched pathways
rankNet(
cc_merged,
mode = "comparison",
stacked = TRUE,
do.stat = TRUE,
color.use = c("BM" = "#b2182b", "CNS" = "#2166ac")
) +
labs(title = "Signalling pathway enrichment — CNS vs BM")
10. Key pathway deep-dives
# Check which key pathways are present
key_pathways <- c("CXCL", "MHC-II", "COLLAGEN", "LAMININ",
"FN1", "VCAM", "ICAM", "IL17", "TIGIT",
"PD-L1", "CD86", "SPP1", "PROS",
"PDGF", "NOTCH", "WNT")
cns_pathways <- cc_cns@netP$pathways
bm_pathways <- cc_bm@netP$pathways
cat("=== Key pathways detected ===\n")## === Key pathways detected ===for (p in key_pathways) {
in_cns <- p %in% cns_pathways
in_bm <- p %in% bm_pathways
if (in_cns | in_bm) {
cat(sprintf(" %-12s CNS: %-5s BM: %-5s\n",
p,
ifelse(in_cns, "YES", "no"),
ifelse(in_bm, "YES", "no")))
}
}## CXCL CNS: YES BM: YES
## MHC-II CNS: YES BM: YES
## COLLAGEN CNS: YES BM: YES
## LAMININ CNS: YES BM: YES
## FN1 CNS: YES BM: YES
## VCAM CNS: YES BM: YES
## ICAM CNS: YES BM: YES
## PD-L1 CNS: YES BM: YES
## CD86 CNS: YES BM: YES
## SPP1 CNS: YES BM: YES
## PROS CNS: YES BM: YES
## PDGF CNS: YES BM: YES
## NOTCH CNS: YES BM: no
## WNT CNS: YES BM: no# Chord diagrams for top CNS-specific pathways
# Run for first 3 CNS-specific pathways detected
cns_specific <- setdiff(cns_pathways, bm_pathways)
cat("\nCNS-specific pathways (not in BM):\n")##
## CNS-specific pathways (not in BM):print(cns_specific)## [1] "MK" "PTN" "THBS" "CX3C" "ADGRL"
## [6] "12oxoLTB4" "FLRT" "SEMA3" "CD200" "FGF"
## [11] "PTPR" "BMP" "SLIT" "CLDN" "EPHA"
## [16] "ADGRA" "HSPG" "CDH5" "UNC5" "Netrin"
## [21] "NOTCH" "PCDH" "BAFF" "NECTIN" "TENASCIN"
## [26] "ncWNT" "Desmosterol" "EPHB" "ANGPT" "WNT"
## [31] "AGRN" "Adenosine" "RELN" "CDH1"# Plot top 3 CNS-specific pathways
for (p in head(cns_specific, 3)) {
tryCatch({
netVisual_aggregate(
cc_cns,
signaling = p,
layout = "chord",
title.name = paste0(p, " — CNS only")
)
}, error = function(e) {
cat("Could not plot pathway:", p, "\n")
})
}
11. Save
qs_save(cc_cns, paste0(BASE_DIR, "17a_cellchat_cns.qs"))
qs_save(cc_bm, paste0(BASE_DIR, "17a_cellchat_bm.qs"))
qs_save(cc_cns_shared, paste0(BASE_DIR, "17a_cellchat_cns_shared.qs"))
qs_save(cc_bm_shared, paste0(BASE_DIR, "17a_cellchat_bm_shared.qs"))
qs_save(cc_merged, paste0(BASE_DIR, "17a_cellchat_merged.qs"))
# Save key interaction tables as CSVs
write.csv(
subsetCommunication(cc_cns),
paste0(BASE_DIR, "17a_interactions_CNS_all.csv")
)
write.csv(
subsetCommunication(cc_bm),
paste0(BASE_DIR, "17a_interactions_BM_all.csv")
)
write.csv(
stroma_to_tcell,
paste0(BASE_DIR, "17a_interactions_stroma_to_Tcell_CNS.csv")
)
cat("Saved:\n")## Saved:cat(" 17a_cellchat_cns.qs\n")## 17a_cellchat_cns.qscat(" 17a_cellchat_bm.qs\n")## 17a_cellchat_bm.qscat(" 17a_cellchat_merged.qs\n")## 17a_cellchat_merged.qscat(" 17a_interactions_CNS_all.csv\n")## 17a_interactions_CNS_all.csvcat(" 17a_interactions_BM_all.csv\n")## 17a_interactions_BM_all.csvcat(" 17a_interactions_stroma_to_Tcell_CNS.csv\n\n")## 17a_interactions_stroma_to_Tcell_CNS.csvcat("Use 17a_cellchat_merged.qs as input for Script 17b (NicheNet)\n")## Use 17a_cellchat_merged.qs as input for Script 17b (NicheNet)Session Info
sessionInfo()## R version 4.4.2 (2024-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS 26.3.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Europe/London
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] future_1.70.0 RColorBrewer_1.1-3 tidyr_1.3.2
## [4] qs2_0.1.7 patchwork_1.3.2 Seurat_5.4.0
## [7] SeuratObject_5.3.0 sp_2.2-1 CellChat_2.1.2
## [10] Biobase_2.66.0 BiocGenerics_0.52.0 ggplot2_4.0.2
## [13] igraph_2.2.2 dplyr_1.2.0
##
## loaded via a namespace (and not attached):
## [1] RcppAnnoy_0.0.23 splines_4.4.2 later_1.4.8
## [4] tibble_3.3.1 polyclip_1.10-7 ggnetwork_0.5.14
## [7] fastDummies_1.7.5 lifecycle_1.0.5 rstatix_0.7.3
## [10] doParallel_1.0.17 globals_0.19.1 lattice_0.22-9
## [13] MASS_7.3-65 backports_1.5.0 magrittr_2.0.4
## [16] plotly_4.12.0 sass_0.4.10 rmarkdown_2.30
## [19] jquerylib_0.1.4 yaml_2.3.12 httpuv_1.6.16
## [22] otel_0.2.0 NMF_0.28 sctransform_0.4.3
## [25] spam_2.11-3 spatstat.sparse_3.1-0 reticulate_1.45.0
## [28] cowplot_1.2.0 pbapply_1.7-4 abind_1.4-8
## [31] Rtsne_0.17 presto_1.0.0 purrr_1.2.1
## [34] circlize_0.4.17 IRanges_2.40.1 S4Vectors_0.44.0
## [37] ggrepel_0.9.8 irlba_2.3.7 listenv_0.10.1
## [40] spatstat.utils_3.2-2 goftest_1.2-3 RSpectra_0.16-2
## [43] spatstat.random_3.4-4 fitdistrplus_1.2-6 parallelly_1.46.1
## [46] svglite_2.2.2 codetools_0.2-20 tidyselect_1.2.1
## [49] shape_1.4.6.1 farver_2.1.2 matrixStats_1.5.0
## [52] stats4_4.4.2 spatstat.explore_3.7-0 jsonlite_2.0.0
## [55] GetoptLong_1.1.0 BiocNeighbors_2.0.1 progressr_0.18.0
## [58] Formula_1.2-5 ggridges_0.5.7 ggalluvial_0.12.6
## [61] survival_3.8-6 iterators_1.0.14 systemfonts_1.3.2
## [64] foreach_1.5.2 tools_4.4.2 sna_2.8
## [67] ica_1.0-3 Rcpp_1.1.1 glue_1.8.0
## [70] gridExtra_2.3 xfun_0.56 withr_3.0.2
## [73] BiocManager_1.30.27 fastmap_1.2.0 digest_0.6.39
## [76] R6_2.6.1 mime_0.13 textshaping_1.0.5
## [79] colorspace_2.1-2 Cairo_1.7-0 scattermore_1.2
## [82] tensor_1.5.1 dichromat_2.0-0.1 spatstat.data_3.1-9
## [85] utf8_1.2.6 generics_0.1.4 data.table_1.18.2.1
## [88] FNN_1.1.4.1 httr_1.4.8 htmlwidgets_1.6.4
## [91] uwot_0.2.4 pkgconfig_2.0.3 gtable_0.3.6
## [94] registry_0.5-1 ComplexHeatmap_2.25.3 lmtest_0.9-40
## [97] S7_0.2.1 htmltools_0.5.9 carData_3.0-6
## [100] dotCall64_1.2 clue_0.3-67 scales_1.4.0
## [103] png_0.1-9 spatstat.univar_3.1-6 knitr_1.51
## [106] rstudioapi_0.18.0 reshape2_1.4.5 rjson_0.2.23
## [109] coda_0.19-4.1 statnet.common_4.13.0 nlme_3.1-168
## [112] cachem_1.1.0 zoo_1.8-15 GlobalOptions_0.1.3
## [115] stringr_1.6.0 KernSmooth_2.23-26 parallel_4.4.2
## [118] miniUI_0.1.2 pillar_1.11.1 grid_4.4.2
## [121] vctrs_0.7.1 RANN_2.6.2 promises_1.5.0
## [124] ggpubr_0.6.3 stringfish_0.18.0 car_3.1-5
## [127] xtable_1.8-8 cluster_2.1.8.2 evaluate_1.0.5
## [130] cli_3.6.5 compiler_4.4.2 rlang_1.1.7
## [133] crayon_1.5.3 rngtools_1.5.2 future.apply_1.20.2
## [136] ggsignif_0.6.4 labeling_0.4.3 plyr_1.8.9
## [139] stringi_1.8.7 viridisLite_0.4.3 network_1.20.0
## [142] deldir_2.0-4 gridBase_0.4-7 lazyeval_0.2.2
## [145] spatstat.geom_3.7-0 Matrix_1.7-4 RcppHNSW_0.6.0
## [148] shiny_1.13.0 ROCR_1.0-12 broom_1.0.12
## [151] RcppParallel_5.1.11-2 bslib_0.10.0