Script 17b — NicheNet: CD8 T cells as receivers in the CNS leukaemia niche
Alasdair Duguid
2026-03-18
Overview
NicheNet analysis asking: which ligands from the CNS sender field best explain the CNS-specific transcriptional state of CD8_Tex and CD8_EffectorMemory cells?
Rationale: CellChat (Script 17a) identified which interactions exist. NicheNet asks which of those interactions functionally explain the receiver transcriptional state. Using CNS vs BM differential expression in CD8 T cells as the gene set of interest ensures the analysis is non-circular: CNS-specific senders must earn their ranking by explaining receiver DE rather than being assumed to matter.
Gene set strategy (informed by DE power):
- CD8_EffectorMemory (183 CNS / 307 BM cells): well-powered DE. Gene set = padj < 0.05, log2FC > 0 (15 genes). High-confidence focused set giving AUPR > 0.10 for top ligands.
- CD8_Tex (15 CNS / 99 BM cells): severely underpowered for multiple-testing correction. Gene set = p < 0.01, log2FC > 0 (238 genes). Treated as exploratory; results should be interpreted with caution given low CNS cell numbers.
1. Libraries and prior model
library(nichenetr)
library(Seurat)
library(qs2)
library(tidyverse)
library(ggplot2)
library(ggrepel)
library(patchwork)
library(RColorBrewer)
library(viridis)
library(ComplexHeatmap)
library(circlize)# ── NicheNet v2 prior model (mouse) ──────────────────────────────────────────
# Files from Zenodo: https://zenodo.org/record/7074291
# Download once to scratch; on Mac use ~/Desktop/nichenet_prior
prior_path <- "/Volumes/Al_seq/scratch_backup/nichenet_prior/"
# prior_path <- "~/Desktop/nichenet_prior" # uncomment for local run
# Download block — run once then comment out:
# dir.create(prior_path, recursive = TRUE, showWarnings = FALSE)
# options(timeout = 600)
# download.file(
# "https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final_mouse.rds",
# destfile = file.path(prior_path, "ligand_target_matrix_mouse.rds"),
# method = "curl"
# )
# download.file(
# "https://zenodo.org/record/7074291/files/lr_network_mouse_21122021.rds",
# destfile = file.path(prior_path, "lr_network_mouse.rds"),
# method = "curl"
# )
# download.file(
# "https://zenodo.org/record/7074291/files/weighted_networks_nsga2r_final_mouse.rds",
# destfile = file.path(prior_path, "weighted_networks_mouse.rds"),
# method = "curl"
# )
ligand_target_matrix <- readRDS(
file.path(prior_path, "ligand_target_matrix_mouse.rds")
)
lr_network <- readRDS(
file.path(prior_path, "lr_network_mouse.rds")
)
weighted_networks <- readRDS(
file.path(prior_path, "weighted_networks_mouse.rds")
)
cat("Prior model loaded.\n")## Prior model loaded.cat("Ligand-target matrix:", nrow(ligand_target_matrix), "targets x",
ncol(ligand_target_matrix), "ligands\n")## Ligand-target matrix: 22521 targets x 1287 ligands2. Load and prepare Seurat object
seu <- qs_read("/Volumes/Al_seq/scratch_backup/harmony_clustering/15c_immune_unified_annotated.qs")
# Tissue column is "Tissue" with capital T — create "niche" alias
seu$niche <- seu$Tissue
stopifnot("annotation_unified" %in% colnames(seu@meta.data))
stopifnot("niche" %in% colnames(seu@meta.data))
cat("Object loaded:", ncol(seu), "cells,",
length(unique(seu$annotation_unified)), "populations\n")## Object loaded: 17699 cells, 37 populationstable(seu$niche)##
## BM CNS
## 12639 5060# Assay is "originalexp" (not "RNA")
DefaultAssay(seu) <- "originalexp"
cat("Default assay:", DefaultAssay(seu), "\n")## Default assay: originalexpcat("Data slot range:", range(seu@assays$originalexp@data@x[1:1000]), "\n")## Data slot range: 1.614578 5.9699833. Define populations and helper function
cns_senders <- c(
"BAM_1 (Marco+/activated)",
"BAM_2 (IFN-responsive)",
"BAM_3 (transitional)",
"BM Macrophages",
"BM Macrophages (activated)",
"CD4_CTL",
"CD4_other",
"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"
)
cat("CNS sender populations:", length(cns_senders), "\n")## CNS sender populations: 31# Expressed genes: ≥pct_threshold of cells in population express the gene
get_expressed_genes_pop <- function(seurat_obj, population, niche_label,
pct_threshold = 0.10) {
cells <- WhichCells(
seurat_obj,
expression = annotation_unified == population & niche == niche_label
)
if (length(cells) < 5) return(character(0))
mat <- seurat_obj@assays$originalexp@data[, cells, drop = FALSE]
pct_expr <- rowMeans(mat > 0)
names(pct_expr[pct_expr >= pct_threshold])
}all_ligands <- unique(lr_network$from)
all_receptors <- unique(lr_network$to)
cat("NicheNet ligands:", length(all_ligands),
"| receptors:", length(all_receptors), "\n")## NicheNet ligands: 1287 | receptors: 10844. Expressed genes across CNS senders
cat("Computing expressed genes across CNS sender populations...\n")## Computing expressed genes across CNS sender populations...sender_expressed <- lapply(cns_senders, function(pop) {
get_expressed_genes_pop(seu, pop, "CNS", pct_threshold = 0.10)
})
names(sender_expressed) <- cns_senders
all_sender_expressed <- unique(unlist(sender_expressed))
expressed_ligands <- intersect(all_sender_expressed, all_ligands)
cat("Expressed ligands across CNS senders:", length(expressed_ligands), "\n")## Expressed ligands across CNS senders: 5815. Analysis 1 — CD8_EffectorMemory as receiver
5.1 Differential gene expression (CNS vs BM)
Idents(seu) <- "annotation_unified"
seu_em <- subset(seu, idents = "CD8_EffectorMemory")
cat("CD8_EM cells — CNS:", sum(seu_em$niche == "CNS"),
"| BM:", sum(seu_em$niche == "BM"), "\n")## CD8_EM cells — CNS: 183 | BM: 307Idents(seu_em) <- "niche"
em_de <- FindMarkers(
seu_em,
ident.1 = "CNS",
ident.2 = "BM",
min.pct = 0.10,
logfc.threshold = 0.25,
test.use = "wilcox"
) %>%
rownames_to_column("gene") %>%
arrange(p_val_adj, desc(avg_log2FC))
write.csv(em_de, "17b_CD8EM_CNS_vs_BM_DE.csv", row.names = FALSE)
cat("padj < 0.05, FC > 0:", sum(em_de$p_val_adj < 0.05 & em_de$avg_log2FC > 0), "\n")## padj < 0.05, FC > 0: 15cat("p < 0.01, FC > 0:", sum(em_de$p_val < 0.01 & em_de$avg_log2FC > 0), "\n")## p < 0.01, FC > 0: 181# Gene set: strict (padj < 0.05) — best AUPR in testing
em_geneset <- em_de %>%
filter(p_val_adj < 0.05, avg_log2FC > 0) %>%
pull(gene)
em_geneset_filtered <- intersect(em_geneset, rownames(ligand_target_matrix))
cat("Gene set (CD8_EM, strict padj < 0.05):", length(em_geneset_filtered),
"genes in NicheNet universe\n")## Gene set (CD8_EM, strict padj < 0.05): 15 genes in NicheNet universecat("\nCNS-upregulated genes in CD8_EM (padj < 0.05):\n")##
## CNS-upregulated genes in CD8_EM (padj < 0.05):print(em_de %>% filter(p_val_adj < 0.05, avg_log2FC > 0))## gene p_val avg_log2FC pct.1 pct.2 p_val_adj
## 1 Id2 4.319906e-15 0.8156465 0.945 0.752 1.153156e-10
## 2 Cxcr6 3.124094e-12 0.9242335 0.847 0.577 8.339457e-08
## 3 Fth1 3.816719e-12 0.6943322 0.973 0.902 1.018835e-07
## 4 Bcl2 1.046163e-09 1.3182757 0.563 0.303 2.792629e-05
## 5 Cish 1.781944e-09 1.3710511 0.415 0.160 4.756721e-05
## 6 Pdcd1 1.840247e-09 0.8305345 0.770 0.524 4.912355e-05
## 7 Fgl2 6.268191e-09 1.2762776 0.443 0.192 1.673231e-04
## 8 Pde3b 6.431276e-09 1.2274235 0.497 0.244 1.716765e-04
## 9 Cd3g 3.816634e-08 0.4522119 0.984 0.971 1.018812e-03
## 10 Itga1 4.338441e-08 4.3194810 0.115 0.007 1.158104e-03
## 11 Cd8a 1.165649e-07 0.9565662 0.678 0.554 3.111584e-03
## 12 Shisa5 2.212270e-07 0.4280536 1.000 0.964 5.905433e-03
## 13 Ifngr1 6.092367e-07 1.0346567 0.781 0.632 1.626296e-02
## 14 Socs1 8.065723e-07 1.3870658 0.328 0.150 2.153064e-02
## 15 Derl2 1.156763e-06 1.5204706 0.268 0.101 3.087864e-025.2 Expressed receptors and potential ligands
em_cns_expressed <- get_expressed_genes_pop(seu, "CD8_EffectorMemory", "CNS")
expressed_receptors_em <- intersect(em_cns_expressed, all_receptors)
cat("Expressed receptors in CD8_EM (CNS):", length(expressed_receptors_em), "\n")## Expressed receptors in CD8_EM (CNS): 152em_potential_ligands <- lr_network %>%
filter(from %in% expressed_ligands, to %in% expressed_receptors_em) %>%
pull(from) %>%
unique()
cat("Potential ligands (LR-filtered):", length(em_potential_ligands), "\n")## Potential ligands (LR-filtered): 288em_background_genes <- intersect(em_cns_expressed, rownames(ligand_target_matrix))
cat("Background gene set:", length(em_background_genes), "\n")## Background gene set: 52515.3 Ligand activity scoring
em_ligand_activities <- predict_ligand_activities(
geneset = em_geneset_filtered,
background_expressed_genes = em_background_genes,
ligand_target_matrix = ligand_target_matrix,
potential_ligands = em_potential_ligands
) %>%
arrange(desc(aupr_corrected)) %>%
mutate(rank = row_number())
write.csv(em_ligand_activities, "17b_CD8EM_ligand_activities.csv", row.names = FALSE)
cat("Top 20 ligands for CD8_EM:\n")## Top 20 ligands for CD8_EM:print(head(em_ligand_activities, 20))## # A tibble: 20 × 6
## test_ligand auroc aupr aupr_corrected pearson rank
## <chr> <dbl> <dbl> <dbl> <dbl> <int>
## 1 Cxcl16 0.601 0.141 0.138 0.194 1
## 2 Lamb2 0.572 0.140 0.137 0.0593 2
## 3 Wbp1 0.537 0.137 0.134 0.0325 3
## 4 Efnb2 0.590 0.112 0.109 0.0633 4
## 5 Col5a3 0.600 0.111 0.108 0.0693 5
## 6 Col27a1 0.562 0.111 0.108 0.0575 6
## 7 Pvr 0.561 0.111 0.108 0.0568 7
## 8 Col6a2 0.580 0.110 0.108 0.0570 8
## 9 AU040320 0.562 0.110 0.108 0.0607 9
## 10 Cd55 0.552 0.110 0.107 0.0528 10
## 11 Cd55b 0.552 0.110 0.107 0.0526 11
## 12 Nectin1 0.562 0.110 0.107 0.0535 12
## 13 Col18a1 0.613 0.110 0.107 0.0590 13
## 14 Sema3b 0.559 0.110 0.107 0.0567 14
## 15 Nectin4 0.527 0.110 0.107 0.0476 15
## 16 Cd96 0.544 0.110 0.107 0.0562 16
## 17 Col3a1 0.558 0.110 0.107 0.0511 17
## 18 Col1a2 0.529 0.110 0.107 0.0458 18
## 19 Gzmc 0.604 0.110 0.107 0.0605 19
## 20 Adam15 0.555 0.109 0.107 0.0552 205.4 Visualise — CD8_EffectorMemory
n_top <- 25
top_em_ligands <- em_ligand_activities %>%
top_n(n_top, aupr_corrected) %>%
pull(test_ligand)em_ligand_activities %>%
top_n(n_top, aupr_corrected) %>%
mutate(test_ligand = fct_reorder(test_ligand, aupr_corrected)) %>%
ggplot(aes(x = aupr_corrected, y = test_ligand)) +
geom_point(aes(colour = aupr_corrected), size = 4) +
scale_colour_viridis(option = "plasma", name = "AUPR\n(corrected)") +
labs(
title = "NicheNet ligand activity — CD8_EffectorMemory receiver",
subtitle = "All CNS senders → CD8_EM (CNS-upregulated gene set, padj < 0.05)",
x = "Ligand activity (AUPR corrected)",
y = NULL
) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank())
active_ligand_target_links_em <- em_ligand_activities %>%
top_n(n_top, aupr_corrected) %>%
pull(test_ligand) %>%
get_weighted_ligand_target_links(
geneset = em_geneset_filtered,
ligand_target_matrix = ligand_target_matrix,
n = 200
)
if (nrow(active_ligand_target_links_em) > 0 &&
any(!is.na(active_ligand_target_links_em$weight))) {
hm_mat_em <- active_ligand_target_links_em %>%
filter(!is.na(weight)) %>%
spread(target, weight, fill = 0) %>%
column_to_rownames("ligand") %>%
as.matrix()
ligand_order_em <- em_ligand_activities %>%
filter(test_ligand %in% rownames(hm_mat_em)) %>%
arrange(desc(aupr_corrected)) %>%
pull(test_ligand)
hm_mat_em <- hm_mat_em[ligand_order_em, , drop = FALSE]
activity_vals_em <- em_ligand_activities %>%
filter(test_ligand %in% ligand_order_em) %>%
arrange(match(test_ligand, ligand_order_em)) %>%
pull(aupr_corrected)
row_ha_em <- rowAnnotation(
AUPR = activity_vals_em,
col = list(AUPR = colorRamp2(
c(min(activity_vals_em), max(activity_vals_em)),
c("white", "#8B1A1A")
)),
annotation_name_side = "top"
)
mat_max_em <- max(hm_mat_em)
if (mat_max_em == 0) mat_max_em <- 1 # guard against all-zero matrix
print(Heatmap(
hm_mat_em,
name = "Regulatory\npotential",
col = colorRamp2(c(0, mat_max_em), c("white", "#2166AC")),
right_annotation = row_ha_em,
cluster_rows = FALSE,
cluster_columns = TRUE,
show_column_names = TRUE,
column_names_gp = gpar(fontsize = 7),
row_names_gp = gpar(fontsize = 9),
column_title = "CD8_EM: top ligand → target gene regulatory potential",
heatmap_legend_param = list(direction = "horizontal")
))
} else {
cat("Insufficient ligand-target links to draw heatmap for CD8_EM.\n")
}## Insufficient ligand-target links to draw heatmap for CD8_EM.em_lr_pairs <- lr_network %>%
filter(from %in% top_em_ligands, to %in% expressed_receptors_em) %>%
rename(ligand = from, receptor = to) %>%
left_join(
em_ligand_activities %>% select(test_ligand, aupr_corrected),
by = c("ligand" = "test_ligand")
)
em_cns_cells <- WhichCells(seu, expression = annotation_unified == "CD8_EffectorMemory" & niche == "CNS")
receptor_genes <- unique(em_lr_pairs$receptor)
receptor_genes <- receptor_genes[receptor_genes %in% rownames(seu)]
receptor_pct_em <- rowMeans(
seu@assays$originalexp@data[receptor_genes, em_cns_cells, drop = FALSE] > 0
)
em_lr_pairs <- em_lr_pairs %>%
filter(receptor %in% names(receptor_pct_em)) %>%
mutate(receptor_pct = receptor_pct_em[receptor])
em_lr_pairs %>%
mutate(
ligand = fct_reorder(ligand, aupr_corrected),
receptor = factor(receptor)
) %>%
ggplot(aes(x = receptor, y = ligand)) +
geom_point(aes(colour = aupr_corrected, size = receptor_pct)) +
scale_colour_viridis(option = "plasma", name = "Ligand\nactivity\n(AUPR)") +
scale_size_continuous(name = "Receptor\n% expressed\nin CD8_EM", range = c(2, 8)) +
labs(
title = "CD8_EM: top ligand–receptor pairs",
subtitle = "Ligand activity (colour) × receptor expression in CD8_EM CNS (size)",
x = NULL, y = "Ligand"
) +
theme_bw(base_size = 11) +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid.minor = element_blank())
sender_ligand_pct_em <- lapply(cns_senders, function(pop) {
pop_cells <- WhichCells(seu, expression = annotation_unified == pop & niche == "CNS")
if (length(pop_cells) < 5) return(NULL)
ligs <- top_em_ligands[top_em_ligands %in% rownames(seu)]
if (length(ligs) == 0) return(NULL)
mat <- seu@assays$originalexp@data[ligs, pop_cells, drop = FALSE]
pct <- rowMeans(mat > 0)
data.frame(population = pop, ligand = names(pct), pct_expr = as.numeric(pct))
}) %>% bind_rows()
sender_ligand_pct_em %>%
filter(pct_expr > 0.05) %>%
mutate(
ligand = factor(ligand, levels = rev(top_em_ligands)),
population = factor(population)
) %>%
ggplot(aes(x = population, y = ligand)) +
geom_point(aes(size = pct_expr, colour = pct_expr)) +
scale_colour_distiller(palette = "YlOrRd", direction = 1,
name = "% cells\nexpressing") +
scale_size_continuous(name = "% cells\nexpressing", range = c(1, 8)) +
labs(
title = "Sender attribution — CD8_EM top ligands",
subtitle = "Which CNS populations express each top-ranked ligand?",
x = NULL, y = "Ligand"
) +
theme_bw(base_size = 10) +
theme(axis.text.x = element_text(angle = 55, hjust = 1, size = 8),
panel.grid.minor = element_blank())
6. Analysis 2 — CD8_Tex as receiver (exploratory)
Caveat: Only 15 CD8_Tex cells in CNS. Multiple-testing correction eliminates almost all DE genes. Gene set uses uncorrected p < 0.01 (238 genes). Results are exploratory and should not be over-interpreted.
6.1 Differential gene expression (CNS vs BM)
Idents(seu) <- "annotation_unified"
seu_tex <- subset(seu, idents = "CD8_Tex")
cat("CD8_Tex cells — CNS:", sum(seu_tex$niche == "CNS"),
"| BM:", sum(seu_tex$niche == "BM"), "\n")## CD8_Tex cells — CNS: 15 | BM: 99Idents(seu_tex) <- "niche"
tex_de <- FindMarkers(
seu_tex,
ident.1 = "CNS",
ident.2 = "BM",
min.pct = 0.10,
logfc.threshold = 0.25,
test.use = "wilcox"
) %>%
rownames_to_column("gene") %>%
arrange(p_val_adj, desc(avg_log2FC))
write.csv(tex_de, "17b_CD8Tex_CNS_vs_BM_DE.csv", row.names = FALSE)
cat("padj < 0.05, FC > 0:", sum(tex_de$p_val_adj < 0.05 & tex_de$avg_log2FC > 0), "\n")## padj < 0.05, FC > 0: 1cat("p < 0.01, FC > 0:", sum(tex_de$p_val < 0.01 & tex_de$avg_log2FC > 0), "\n")## p < 0.01, FC > 0: 253# Gene set: middle ground (p < 0.01) — strict gives only 1 gene
tex_geneset <- tex_de %>%
filter(p_val < 0.01, avg_log2FC > 0) %>%
pull(gene)
tex_geneset_filtered <- intersect(tex_geneset, rownames(ligand_target_matrix))
cat("Gene set (CD8_Tex, p < 0.01):", length(tex_geneset_filtered),
"genes in NicheNet universe\n")## Gene set (CD8_Tex, p < 0.01): 238 genes in NicheNet universe6.2 Expressed receptors and potential ligands
tex_cns_expressed <- get_expressed_genes_pop(seu, "CD8_Tex", "CNS")
expressed_receptors_tex <- intersect(tex_cns_expressed, all_receptors)
cat("Expressed receptors in CD8_Tex (CNS):", length(expressed_receptors_tex), "\n")## Expressed receptors in CD8_Tex (CNS): 155tex_potential_ligands <- lr_network %>%
filter(from %in% expressed_ligands, to %in% expressed_receptors_tex) %>%
pull(from) %>%
unique()
cat("Potential ligands (LR-filtered):", length(tex_potential_ligands), "\n")## Potential ligands (LR-filtered): 309tex_background_genes <- intersect(tex_cns_expressed, rownames(ligand_target_matrix))
cat("Background gene set:", length(tex_background_genes), "\n")## Background gene set: 55406.3 Ligand activity scoring
tex_ligand_activities <- predict_ligand_activities(
geneset = tex_geneset_filtered,
background_expressed_genes = tex_background_genes,
ligand_target_matrix = ligand_target_matrix,
potential_ligands = tex_potential_ligands
) %>%
arrange(desc(aupr_corrected)) %>%
mutate(rank = row_number())
write.csv(tex_ligand_activities, "17b_CD8Tex_ligand_activities.csv", row.names = FALSE)
cat("Top 20 ligands for CD8_Tex:\n")## Top 20 ligands for CD8_Tex:print(head(tex_ligand_activities, 20))## # A tibble: 20 × 6
## test_ligand auroc aupr aupr_corrected pearson rank
## <chr> <dbl> <dbl> <dbl> <dbl> <int>
## 1 Bmp2 0.509 0.0519 0.00897 0.0422 1
## 2 Il1b 0.511 0.0493 0.00633 0.0249 2
## 3 Col11a1 0.496 0.0492 0.00628 0.0294 3
## 4 Tgfb1 0.504 0.0492 0.00624 0.0295 4
## 5 Cd40lg 0.499 0.0491 0.00616 0.00625 5
## 6 Cd5l 0.499 0.0488 0.00580 -0.000353 6
## 7 Col5a1 0.482 0.0486 0.00568 -0.00218 7
## 8 Vcan 0.483 0.0483 0.00534 -0.00426 8
## 9 Ccn3 0.488 0.0482 0.00524 0.0118 9
## 10 Cxcl12 0.493 0.0480 0.00501 0.0139 10
## 11 Thy1 0.479 0.0478 0.00483 -0.00660 11
## 12 Tslp 0.473 0.0475 0.00456 0.00578 12
## 13 Il10 0.507 0.0474 0.00442 0.0148 13
## 14 Il1a 0.502 0.0473 0.00436 0.0161 14
## 15 Col4a1 0.487 0.0471 0.00415 -0.00985 15
## 16 Tgfb2 0.492 0.0471 0.00413 0.0173 16
## 17 Mill2 0.475 0.0467 0.00374 -0.0103 17
## 18 Igf2 0.498 0.0466 0.00366 0.0101 18
## 19 Igf1 0.499 0.0465 0.00357 0.0161 19
## 20 Ccl12 0.516 0.0465 0.00355 0.00918 206.4 Visualise — CD8_Tex
n_top_tex <- 25
top_tex_ligands <- tex_ligand_activities %>%
top_n(n_top_tex, aupr_corrected) %>%
pull(test_ligand)tex_ligand_activities %>%
top_n(n_top_tex, aupr_corrected) %>%
mutate(test_ligand = fct_reorder(test_ligand, aupr_corrected)) %>%
ggplot(aes(x = aupr_corrected, y = test_ligand)) +
geom_point(aes(colour = aupr_corrected), size = 4) +
scale_colour_viridis(option = "plasma", name = "AUPR\n(corrected)") +
labs(
title = "NicheNet ligand activity — CD8_Tex receiver (exploratory)",
subtitle = "All CNS senders → CD8_Tex (p < 0.01 gene set; n=15 CNS cells)",
x = "Ligand activity (AUPR corrected)",
y = NULL
) +
theme_bw(base_size = 12) +
theme(panel.grid.minor = element_blank())
active_ligand_target_links_tex <- tex_ligand_activities %>%
top_n(n_top_tex, aupr_corrected) %>%
pull(test_ligand) %>%
get_weighted_ligand_target_links(
geneset = tex_geneset_filtered,
ligand_target_matrix = ligand_target_matrix,
n = 200
)
if (nrow(active_ligand_target_links_tex) > 0 &&
any(!is.na(active_ligand_target_links_tex$weight))) {
hm_mat_tex <- active_ligand_target_links_tex %>%
filter(!is.na(weight)) %>%
spread(target, weight, fill = 0) %>%
column_to_rownames("ligand") %>%
as.matrix()
ligand_order_tex <- tex_ligand_activities %>%
filter(test_ligand %in% rownames(hm_mat_tex)) %>%
arrange(desc(aupr_corrected)) %>%
pull(test_ligand)
hm_mat_tex <- hm_mat_tex[ligand_order_tex, , drop = FALSE]
activity_vals_tex <- tex_ligand_activities %>%
filter(test_ligand %in% ligand_order_tex) %>%
arrange(match(test_ligand, ligand_order_tex)) %>%
pull(aupr_corrected)
row_ha_tex <- rowAnnotation(
AUPR = activity_vals_tex,
col = list(AUPR = colorRamp2(
c(min(activity_vals_tex), max(activity_vals_tex)),
c("white", "#8B1A1A")
)),
annotation_name_side = "top"
)
mat_max_tex <- max(hm_mat_tex)
if (mat_max_tex == 0) mat_max_tex <- 1
print(Heatmap(
hm_mat_tex,
name = "Regulatory\npotential",
col = colorRamp2(c(0, mat_max_tex), c("white", "#2166AC")),
right_annotation = row_ha_tex,
cluster_rows = FALSE,
cluster_columns = TRUE,
show_column_names = TRUE,
column_names_gp = gpar(fontsize = 7),
row_names_gp = gpar(fontsize = 9),
column_title = "CD8_Tex: top ligand → target gene regulatory potential (exploratory)",
heatmap_legend_param = list(direction = "horizontal")
))
} else {
cat("Insufficient ligand-target links to draw heatmap for CD8_Tex.\n")
}## Insufficient ligand-target links to draw heatmap for CD8_Tex.sender_ligand_pct_tex <- lapply(cns_senders, function(pop) {
pop_cells <- WhichCells(seu, expression = annotation_unified == pop & niche == "CNS")
if (length(pop_cells) < 5) return(NULL)
ligs <- top_tex_ligands[top_tex_ligands %in% rownames(seu)]
if (length(ligs) == 0) return(NULL)
mat <- seu@assays$originalexp@data[ligs, pop_cells, drop = FALSE]
pct <- rowMeans(mat > 0)
data.frame(population = pop, ligand = names(pct), pct_expr = as.numeric(pct))
}) %>% bind_rows()
sender_ligand_pct_tex %>%
filter(pct_expr > 0.05) %>%
mutate(
ligand = factor(ligand, levels = rev(top_tex_ligands)),
population = factor(population)
) %>%
ggplot(aes(x = population, y = ligand)) +
geom_point(aes(size = pct_expr, colour = pct_expr)) +
scale_colour_distiller(palette = "YlOrRd", direction = 1,
name = "% cells\nexpressing") +
scale_size_continuous(name = "% cells\nexpressing", range = c(1, 8)) +
labs(
title = "Sender attribution — CD8_Tex top ligands (exploratory)",
subtitle = "Which CNS populations express each top-ranked ligand?",
x = NULL, y = "Ligand"
) +
theme_bw(base_size = 10) +
theme(axis.text.x = element_text(angle = 55, hjust = 1, size = 8),
panel.grid.minor = element_blank())
7. Cross-receiver comparison
all_top_ligands <- union(top_em_ligands, top_tex_ligands)
comparison_df <- full_join(
tex_ligand_activities %>%
top_n(n_top_tex, aupr_corrected) %>%
select(test_ligand, aupr_tex = aupr_corrected, rank_tex = rank),
em_ligand_activities %>%
top_n(n_top, aupr_corrected) %>%
select(test_ligand, aupr_em = aupr_corrected, rank_em = rank),
by = "test_ligand"
) %>%
mutate(category = case_when(
!is.na(aupr_tex) & !is.na(aupr_em) ~ "Shared",
!is.na(aupr_tex) & is.na(aupr_em) ~ "CD8_Tex specific",
is.na(aupr_tex) & !is.na(aupr_em) ~ "CD8_EM specific"
))
cat("Shared top ligands:", sum(comparison_df$category == "Shared"), "\n")## Shared top ligands: 0cat("CD8_Tex specific:", sum(comparison_df$category == "CD8_Tex specific"), "\n")## CD8_Tex specific: 25cat("CD8_EM specific:", sum(comparison_df$category == "CD8_EM specific"), "\n")## CD8_EM specific: 25comparison_df %>%
replace_na(list(aupr_tex = 0, aupr_em = 0)) %>%
ggplot(aes(x = aupr_tex, y = aupr_em, label = test_ligand, colour = category)) +
geom_point(size = 3) +
geom_text_repel(size = 3, max.overlaps = 20) +
scale_colour_manual(values = c(
"Shared" = "#2166AC",
"CD8_Tex specific" = "#D6604D",
"CD8_EM specific" = "#4DAC26"
)) +
geom_abline(linetype = "dashed", colour = "grey70") +
labs(
title = "Cross-receiver comparison of top ligands",
subtitle = "Note: CD8_Tex results are exploratory (n=15 CNS cells)",
x = "AUPR (CD8_Tex receiver)",
y = "AUPR (CD8_EffectorMemory receiver)",
colour = NULL
) +
theme_bw(base_size = 12) +
theme(legend.position = "bottom", panel.grid.minor = element_blank())
8. Vγ6Vδ4 ligand contribution
vg6_cells <- WhichCells(seu, expression = annotation_unified == "Vγ6Vδ4" & niche == "CNS")
vg6_expressed <- get_expressed_genes_pop(seu, "Vγ6Vδ4", "CNS", pct_threshold = 0.05)
ligs_in_seu <- all_top_ligands[all_top_ligands %in% rownames(seu)]
vg6_pct <- rowMeans(
seu@assays$originalexp@data[ligs_in_seu, vg6_cells, drop = FALSE] > 0
)
vg6_ligand_df <- data.frame(
ligand = names(vg6_pct),
pct_vg6 = as.numeric(vg6_pct),
in_em_top = names(vg6_pct) %in% top_em_ligands,
in_tex_top = names(vg6_pct) %in% top_tex_ligands
) %>%
left_join(em_ligand_activities %>% select(test_ligand, aupr_em = aupr_corrected),
by = c("ligand" = "test_ligand")) %>%
left_join(tex_ligand_activities %>% select(test_ligand, aupr_tex = aupr_corrected),
by = c("ligand" = "test_ligand")) %>%
filter(pct_vg6 >= 0.05) %>%
arrange(desc(pmax(coalesce(aupr_em, 0), coalesce(aupr_tex, 0))))
write.csv(vg6_ligand_df, "17b_Vg6Vd4_ligand_contribution.csv", row.names = FALSE)
cat("Top ligands expressed by Vγ6Vδ4 (≥5% cells):\n")## Top ligands expressed by Vγ6Vδ4 (≥5% cells):print(vg6_ligand_df)## ligand pct_vg6 in_em_top in_tex_top aupr_em aupr_tex
## 1 Wbp1 0.49586777 TRUE FALSE 0.134299120 0.0006880477
## 2 Pvr 0.08264463 TRUE FALSE 0.107718770 0.0009180495
## 3 AU040320 0.22314050 TRUE FALSE 0.107502784 -0.0002155536
## 4 Cd96 0.60330579 TRUE FALSE 0.106827673 NA
## 5 Ptdss1 0.19008264 TRUE FALSE 0.106009909 -0.0018130575
## 6 Thy1 0.90909091 FALSE TRUE 0.031237172 0.0048262134
## 7 Tgfb1 0.36363636 FALSE TRUE 0.006771057 0.0062432917
## 8 Lta 0.11570248 FALSE TRUE 0.002920690 0.0031127148vg6_ligand_df %>%
pivot_longer(c(aupr_em, aupr_tex), names_to = "receiver", values_to = "aupr") %>%
mutate(
receiver = recode(receiver, aupr_em = "CD8_EM", aupr_tex = "CD8_Tex (exploratory)"),
ligand = fct_reorder(ligand, pct_vg6)
) %>%
filter(!is.na(aupr)) %>%
ggplot(aes(x = aupr, y = ligand, fill = receiver)) +
geom_col(position = "dodge", alpha = 0.85) +
scale_fill_manual(values = c("CD8_EM" = "#4DAC26",
"CD8_Tex (exploratory)" = "#D6604D")) +
labs(
title = "Vγ6Vδ4 ligand contribution to CD8 T cell gene programmes",
subtitle = "Ligands expressed by ≥5% of Vγ6Vδ4 CNS cells",
x = "Ligand activity (AUPR corrected)",
y = "Ligand",
fill = "Receiver"
) +
theme_bw(base_size = 11) +
theme(panel.grid.minor = element_blank())
9. ECM fibroblast ligand contribution
ecm_cells <- WhichCells(
seu,
expression = annotation_unified == "CNS Meningeal Fibroblasts (ECM)" & niche == "CNS"
)
ligs_in_seu <- all_top_ligands[all_top_ligands %in% rownames(seu)]
ecm_pct <- rowMeans(
seu@assays$originalexp@data[ligs_in_seu, ecm_cells, drop = FALSE] > 0
)
ecm_ligand_df <- data.frame(
ligand = names(ecm_pct),
pct_ecm = as.numeric(ecm_pct)
) %>%
left_join(em_ligand_activities %>% select(test_ligand, aupr_em = aupr_corrected),
by = c("ligand" = "test_ligand")) %>%
left_join(tex_ligand_activities %>% select(test_ligand, aupr_tex = aupr_corrected),
by = c("ligand" = "test_ligand")) %>%
filter(pct_ecm >= 0.10) %>%
arrange(desc(pmax(coalesce(aupr_em, 0), coalesce(aupr_tex, 0))))
write.csv(ecm_ligand_df, "17b_ECM_fibroblast_ligand_contribution.csv", row.names = FALSE)
cat("Top ligands expressed by ECM fibroblasts (≥10% cells):\n")## Top ligands expressed by ECM fibroblasts (≥10% cells):print(ecm_ligand_df)## ligand pct_ecm aupr_em aupr_tex
## 1 Cxcl16 0.6945455 0.138362124 -0.0008958682
## 2 Lamb2 0.9309091 0.136858170 -0.0002703705
## 3 Wbp1 0.4181818 0.134299120 0.0006880477
## 4 Col5a3 0.1636364 0.108019781 0.0002705780
## 5 Col27a1 0.4909091 0.107841304 -0.0003524503
## 6 Pvr 0.2436364 0.107718770 0.0009180495
## 7 Col6a2 0.9890909 0.107533923 0.0006180007
## 8 AU040320 0.4145455 0.107502784 -0.0002155536
## 9 Nectin1 0.2327273 0.107095611 -0.0003509692
## 10 Col18a1 0.3454545 0.107000798 -0.0011192963
## 11 Col3a1 0.9854545 0.106826855 0.0012129037
## 12 Col1a2 1.0000000 0.106721477 0.0005813339
## 13 Adam15 0.6290909 0.106531651 0.0003936585
## 14 Ptdss1 0.5090909 0.106009909 -0.0018130575
## 15 Ccn4 0.1672727 0.105878592 -0.0004132507
## 16 Hspg2 0.6654545 0.098294135 0.0002239322
## 17 Sema3d 0.4727273 0.096835330 -0.0011694135
## 18 Ccn3 0.7600000 0.079829224 0.0052391706
## 19 Vcan 0.1345455 0.046093609 0.0053413638
## 20 Col4a1 0.9781818 0.012629320 0.0041539718
## 21 Mill2 0.6290909 0.010748723 0.0037359581
## 22 Col5a1 0.7636364 0.007496596 0.0056787984
## 23 Cxcl12 0.9890909 0.007337244 0.0050127099
## 24 Tgfb1 0.5309091 0.006771057 0.0062432917
## 25 Col11a1 0.7490909 0.005049029 0.0062805424
## 26 Tgfb2 0.7127273 0.004922591 0.0041334188
## 27 Alcam 0.5018182 0.004642349 0.0033252397
## 28 Igf2 0.4472727 0.003254202 0.0036643633
## 29 Igf1 0.6436364 0.001553831 0.003567257410. Save all results
save(
em_ligand_activities,
tex_ligand_activities,
comparison_df,
vg6_ligand_df,
ecm_ligand_df,
em_de,
tex_de,
file = "17b_nichenet_results.RData"
)
cat("Script 17b complete. Outputs:\n")## Script 17b complete. Outputs:cat(" 17b_CD8EM_CNS_vs_BM_DE.csv\n")## 17b_CD8EM_CNS_vs_BM_DE.csvcat(" 17b_CD8Tex_CNS_vs_BM_DE.csv\n")## 17b_CD8Tex_CNS_vs_BM_DE.csvcat(" 17b_CD8EM_ligand_activities.csv\n")## 17b_CD8EM_ligand_activities.csvcat(" 17b_CD8Tex_ligand_activities.csv\n")## 17b_CD8Tex_ligand_activities.csvcat(" 17b_Vg6Vd4_ligand_contribution.csv\n")## 17b_Vg6Vd4_ligand_contribution.csvcat(" 17b_ECM_fibroblast_ligand_contribution.csv\n")## 17b_ECM_fibroblast_ligand_contribution.csvcat(" 17b_nichenet_results.RData\n")## 17b_nichenet_results.RData11. 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] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] circlize_0.4.17 ComplexHeatmap_2.25.3 viridis_0.6.5
## [4] viridisLite_0.4.3 RColorBrewer_1.1-3 patchwork_1.3.2
## [7] ggrepel_0.9.8 lubridate_1.9.5 forcats_1.0.1
## [10] stringr_1.6.0 dplyr_1.2.0 purrr_1.2.1
## [13] readr_2.2.0 tidyr_1.3.2 tibble_3.3.1
## [16] ggplot2_4.0.2 tidyverse_2.0.0 qs2_0.1.7
## [19] Seurat_5.4.0 SeuratObject_5.3.0 sp_2.2-1
## [22] nichenetr_2.2.1.1
##
## loaded via a namespace (and not attached):
## [1] RcppAnnoy_0.0.23 splines_4.4.2 later_1.4.8
## [4] bitops_1.0-9 polyclip_1.10-7 hardhat_1.4.2
## [7] pROC_1.19.0.1 rpart_4.1.24 fastDummies_1.7.5
## [10] lifecycle_1.0.5 doParallel_1.0.17 globals_0.19.1
## [13] lattice_0.22-9 MASS_7.3-65 backports_1.5.0
## [16] magrittr_2.0.4 limma_3.62.2 Hmisc_5.2-5
## [19] plotly_4.12.0 sass_0.4.10 rmarkdown_2.30
## [22] jquerylib_0.1.4 yaml_2.3.12 httpuv_1.6.16
## [25] otel_0.2.0 sctransform_0.4.3 spam_2.11-3
## [28] spatstat.sparse_3.1-0 reticulate_1.45.0 cowplot_1.2.0
## [31] pbapply_1.7-4 abind_1.4-8 Rtsne_0.17
## [34] presto_1.0.0 BiocGenerics_0.52.0 nnet_7.3-20
## [37] tweenr_2.0.3 ipred_0.9-15 gdtools_0.5.0
## [40] lava_1.8.2 S4Vectors_0.44.0 IRanges_2.40.1
## [43] irlba_2.3.7 listenv_0.10.1 spatstat.utils_3.2-2
## [46] goftest_1.2-3 RSpectra_0.16-2 spatstat.random_3.4-4
## [49] fitdistrplus_1.2-6 parallelly_1.46.1 codetools_0.2-20
## [52] ggforce_0.5.0 shape_1.4.6.1 tidyselect_1.2.1
## [55] farver_2.1.2 base64enc_0.1-6 matrixStats_1.5.0
## [58] stats4_4.4.2 spatstat.explore_3.7-0 jsonlite_2.0.0
## [61] caret_7.0-1 GetoptLong_1.1.0 e1071_1.7-17
## [64] Formula_1.2-5 progressr_0.18.0 ggridges_0.5.7
## [67] survival_3.8-6 iterators_1.0.14 systemfonts_1.3.2
## [70] foreach_1.5.2 tools_4.4.2 ggnewscale_0.5.2
## [73] ica_1.0-3 Rcpp_1.1.1 glue_1.8.0
## [76] prodlim_2026.03.11 gridExtra_2.3 xfun_0.56
## [79] withr_3.0.2 fastmap_1.2.0 caTools_1.18.3
## [82] digest_0.6.39 timechange_0.4.0 R6_2.6.1
## [85] mime_0.13 colorspace_2.1-2 scattermore_1.2
## [88] tensor_1.5.1 dichromat_2.0-0.1 spatstat.data_3.1-9
## [91] DiagrammeR_1.0.11 utf8_1.2.6 generics_0.1.4
## [94] fontLiberation_0.1.0 data.table_1.18.2.1 recipes_1.3.1
## [97] class_7.3-23 httr_1.4.8 htmlwidgets_1.6.4
## [100] uwot_0.2.4 ModelMetrics_1.2.2.2 pkgconfig_2.0.3
## [103] gtable_0.3.6 timeDate_4052.112 lmtest_0.9-40
## [106] S7_0.2.1 shadowtext_0.1.6 htmltools_0.5.9
## [109] fontBitstreamVera_0.1.1 dotCall64_1.2 clue_0.3-67
## [112] scales_1.4.0 png_0.1-9 gower_1.0.2
## [115] spatstat.univar_3.1-6 knitr_1.51 rstudioapi_0.18.0
## [118] tzdb_0.5.0 reshape2_1.4.5 rjson_0.2.23
## [121] checkmate_2.3.4 visNetwork_2.1.4 nlme_3.1-168
## [124] proxy_0.4-29 cachem_1.1.0 zoo_1.8-15
## [127] GlobalOptions_0.1.3 KernSmooth_2.23-26 parallel_4.4.2
## [130] miniUI_0.1.2 foreign_0.8-91 pillar_1.11.1
## [133] vctrs_0.7.1 RANN_2.6.2 randomForest_4.7-1.2
## [136] promises_1.5.0 stringfish_0.18.0 xtable_1.8-8
## [139] cluster_2.1.8.2 htmlTable_2.4.3 evaluate_1.0.5
## [142] cli_3.6.5 compiler_4.4.2 rlang_1.1.7
## [145] crayon_1.5.3 future.apply_1.20.2 labeling_0.4.3
## [148] fdrtool_1.2.18 plyr_1.8.9 ggiraph_0.9.6
## [151] stringi_1.8.7 deldir_2.0-4 lazyeval_0.2.2
## [154] spatstat.geom_3.7-0 fontquiver_0.2.1 Matrix_1.7-4
## [157] RcppHNSW_0.6.0 hms_1.1.4 future_1.70.0
## [160] statmod_1.5.1 shiny_1.13.0 ROCR_1.0-12
## [163] igraph_2.2.2 RcppParallel_5.1.11-2 bslib_0.10.0