14: Stromal CAF Scoring & IL-17 Responsiveness

CAF subtype characterisation, pro-leukemia niche scoring, and IL-17A response in CNS fibroblasts

Overview

This script characterises the stromal compartment in the context of cancer-associated fibroblast (CAF) biology and IL-17A responsiveness. The central question is whether CNS fibroblast populations adopt pro-leukemic states, and whether these states are consistent with activation downstream of Vγ6Vδ4-derived IL-17A.

CAF scoring is restricted to the three CNS fibroblast populations: - CNS Leptomeningeal Fibroblasts - CNS Meningeal Fibroblasts (ECM) - CNS Mesenchymal/Adventitial

BM Osteoblastic Stroma and Choroid Plexus Epithelium are retained for contextual visualisation but excluded from CAF and IL-17 scoring.

Analysis sections:

1. Object inspection — confirm subtype annotations and tissue distribution
2. CAF subtype scoring — myoCAF, iCAF, apCAF via UCell (fibroblasts only)
3. Pro-leukemia niche scoring — existing Niche_retention_UCell plus key ligand expression
4. IL-17 receptor expression — Il17ra / Il17rc across subtypes
5. IL-17 downstream response scoring — NF-κB/Act1 output genes via UCell
6. Integration — iCAF vs IL-17 response correlation; combined heatmap
7. Inter-subtype comparison — statistical tests across CNS fibroblast subtypes
8. Export — save scored object and summary tables

---

Setup

library(Seurat)
library(qs2)
library(UCell)
library(dplyr)
library(tidyr)
library(tibble)
library(ggplot2)
library(patchwork)
library(ggrepel)
library(RColorBrewer)
library(knitr)
library(scales)

base_path <- "/Users/alasdairduguid/Documents/Projects/BSH start up grant 2024 - scRNAseq miR128a/scRNAseq_proj/analysis_feb26"
cache_dir  <- file.path(base_path, "data")

# Colour palettes
tissue_cols <- c("BM" = "#2166AC", "CNS" = "#B2182B")

fib_cols <- c(
  "CNS Leptomeningeal Fibroblasts"  = "#E41A1C",
  "CNS Meningeal Fibroblasts (ECM)" = "#377EB8",
  "CNS Mesenchymal/Adventitial"     = "#4DAF4A"
)

all_subtype_cols <- c(fib_cols,
                      "BM Osteoblastic Stroma"    = "#FF7F00",
                      "Choroid Plexus Epithelium" = "#984EA3")

stroma <- qs_read(file.path(cache_dir, "13c_v2_stroma_annotated_clean.qs"))
DefaultAssay(stroma) <- "originalexp"

cat("Stroma object loaded:", ncol(stroma), "cells,", nrow(stroma), "features\n")

## Stroma object loaded: 993 cells, 26694 features

cat("Active assay:", DefaultAssay(stroma), "\n")

## Active assay: originalexp

---

1. Object Inspection

Interpretation: The stroma object contains 691 CNS fibroblasts across three subtypes. CNS Meningeal Fibroblasts (ECM) are the most abundant (n = 275), followed by Leptomeningeal (n = 294) and Mesenchymal/Adventitial (n = 122). All three populations are CNS-restricted, consistent with the near-complete tissue segregation observed in Script 13c. The UMAP shows clean subtype separation, validating the annotation. Identity markers distinguish populations along expected lines: ECM fibroblasts are enriched for structural matrix genes (Col1a1, Fn1, Postn); Leptomeningeal fibroblasts for solute transporter and immune-modulatory genes (Slc47a1, Dcn); Mesenchymal/Adventitial for perivascular and progenitor markers (Pi16, Cd34).

Cell Counts by Subtype and Tissue

stroma@meta.data %>%
  count(subtype, Tissue) %>%
  pivot_wider(names_from = Tissue, values_from = n, values_fill = 0) %>%
  mutate(Total = rowSums(select(., -subtype))) %>%
  arrange(desc(Total)) %>%
  kable(caption = "Cell counts per stromal subtype and tissue")

Cell counts per stromal subtype and tissue

subtype	BM	CNS	Total
CNS Leptomeningeal Fibroblasts	0	294	294
BM Osteoblastic Stroma	277	7	284
CNS Meningeal Fibroblasts (ECM)	0	275	275
CNS Mesenchymal/Adventitial	0	122	122
Choroid Plexus Epithelium	0	18	18

UMAP Overview

p1 <- DimPlot(stroma, group.by = "subtype", label = TRUE,
              label.size = 3, repel = TRUE,
              cols = all_subtype_cols) +
  ggtitle("Stromal Subtypes") +
  theme(legend.text = element_text(size = 8))

p2 <- DimPlot(stroma, group.by = "Tissue", cols = tissue_cols) +
  ggtitle("By Tissue")

p1 + p2

Split by Tissue

DimPlot(stroma, group.by = "subtype", split.by = "Tissue",
        label = TRUE, label.size = 3, repel = TRUE,
        cols = all_subtype_cols) +
  ggtitle("Stromal Subtypes by Tissue") +
  theme(legend.text = element_text(size = 8))

Identity Marker DotPlot

identity_markers <- c(
  # Leptomeningeal / meningeal fibroblasts
  "Col1a1", "Col3a1", "Dcn", "Lum", "Fgfr3", "Pdgfra",
  # Fibroblast activation
  "Acta2", "Tagln", "Fap",
  # Adventitial / Pi16+ fibroblasts
  "Pi16", "Ly6c1", "Cd34",
  # Pericytes (should be absent if clean)
  "Pdgfrb", "Rgs5",
  # Choroid plexus epithelium
  "Ttr", "Folr1", "Aqp1",
  # Osteoblastic
  "Runx2", "Bglap", "Sp7"
)

identity_present <- identity_markers[identity_markers %in% rownames(stroma)]

DotPlot(stroma, features = identity_present, group.by = "subtype",
        cols = c("lightgrey", "darkred"), dot.scale = 6) +
  RotatedAxis() +
  ggtitle("Stromal Identity Markers") +
  theme(axis.text.x = element_text(size = 9))

---

2. CAF Subtype Scoring

Interpretation: CNS Meningeal Fibroblasts (ECM) score highest for myoCAF, consistent with their heavy ECM production programme and expression of Acta2, Tagln, Mmp2, and Fap — positioning them as the structurally active population capable of physically remodelling the leukaemic niche. CNS Leptomeningeal and Mesenchymal/Adventitial fibroblasts show elevated apCAF scores, with MHC class II gene expression (H2-Ab1, H2-Eb1, Cd74) consistent with the known role of leptomeningeal fibroblasts in meningeal immune surveillance and T cell interaction. This raises the possibility that these populations influence T cell activity in the CNS niche — potentially contributing to an immunosuppressive environment through antigen-specific Treg induction rather than direct niche support. The iCAF programme is relatively low across all populations, suggesting cytokine-driven inflammatory CAF activation is not the dominant stromal state in this model.

Based on the Öhlund/Biffi framework. Scores are calculated on the three CNS fibroblast clusters only — BM Osteoblastic Stroma and Choroid Plexus Epithelium are excluded as CAF biology is not applicable to these lineages.

fibroblast_subtypes <- c(
  "CNS Leptomeningeal Fibroblasts",
  "CNS Meningeal Fibroblasts (ECM)",
  "CNS Mesenchymal/Adventitial"
)

stroma_fib <- subset(stroma, subset = subtype %in% fibroblast_subtypes)

cat("Fibroblast cells retained:", ncol(stroma_fib), "\n")

## Fibroblast cells retained: 691

cat("Excluded:", ncol(stroma) - ncol(stroma_fib),
    "cells (BM Osteoblastic Stroma + Choroid Plexus Epithelium)\n\n")

## Excluded: 302 cells (BM Osteoblastic Stroma + Choroid Plexus Epithelium)

print(table(stroma_fib$subtype, stroma_fib$Tissue))

##                                  
##                                   CNS
##   CNS Leptomeningeal Fibroblasts  294
##   CNS Meningeal Fibroblasts (ECM) 275
##   CNS Mesenchymal/Adventitial     122

myoCAF_genes <- c("Acta2", "Tagln", "Myl9", "Tpm2", "Pdgfra",
                  "Thy1", "Fap", "Postn", "Lox", "Thbs2",
                  "Col10a1", "Col11a1", "Mmp2", "Mmp14")

iCAF_genes <- c("Il6", "Cxcl1", "Cxcl2", "Cxcl12", "Cxcl16",
                "Has1", "Clec3b", "Ly6c1", "Ccl2", "Ccl7",
                "Il33", "Lcn2", "S100a8", "S100a9",
                "Pdpn", "Vcam1", "Icam1", "Tnfsf13b")

apCAF_genes <- c("H2-Ab1", "H2-Eb1", "Cd74", "Slpi",
                 "H2-Aa", "H2-DMa", "H2-DMb1", "Ciita",
                 "Saa3", "Ly6c2")

myoCAF_present <- myoCAF_genes[myoCAF_genes %in% rownames(stroma_fib)]
iCAF_present   <- iCAF_genes[iCAF_genes   %in% rownames(stroma_fib)]
apCAF_present  <- apCAF_genes[apCAF_genes %in% rownames(stroma_fib)]

cat("myoCAF genes found:", length(myoCAF_present), "/", length(myoCAF_genes), "\n")

## myoCAF genes found: 14 / 14

cat("iCAF genes found:  ", length(iCAF_present),   "/", length(iCAF_genes),   "\n")

## iCAF genes found:   18 / 18

cat("apCAF genes found: ", length(apCAF_present),   "/", length(apCAF_genes),  "\n")

## apCAF genes found:  10 / 10

cat("\nmyoCAF missing:", paste(setdiff(myoCAF_genes, rownames(stroma_fib)), collapse = ", "), "\n")

## 
## myoCAF missing:

cat("iCAF missing:  ", paste(setdiff(iCAF_genes,   rownames(stroma_fib)), collapse = ", "), "\n")

## iCAF missing:

cat("apCAF missing: ", paste(setdiff(apCAF_genes,  rownames(stroma_fib)), collapse = ", "), "\n")

## apCAF missing:

stroma_fib <- AddModuleScore_UCell(
  stroma_fib,
  features = list(myoCAF = myoCAF_present,
                  iCAF   = iCAF_present,
                  apCAF  = apCAF_present),
  name = ""
)

cat("CAF scores added. Summary:\n")

## CAF scores added. Summary:

stroma_fib@meta.data %>%
  group_by(subtype, Tissue) %>%
  summarise(
    myoCAF = round(mean(myoCAF), 3),
    iCAF   = round(mean(iCAF),   3),
    apCAF  = round(mean(apCAF),  3),
    n      = n(),
    .groups = "drop"
  ) %>%
  arrange(desc(iCAF)) %>%
  kable(caption = "Mean CAF subtype scores (CNS fibroblasts only)")

Mean CAF subtype scores (CNS fibroblasts only)

subtype	Tissue	myoCAF	iCAF	apCAF	n
CNS Meningeal Fibroblasts (ECM)	CNS	0.177	0.126	0.050	275
CNS Leptomeningeal Fibroblasts	CNS	0.102	0.108	0.114	294
CNS Mesenchymal/Adventitial	CNS	0.116	0.102	0.129	122

CAF Score UMAPs

p1 <- FeaturePlot(stroma_fib, features = "myoCAF",
                  cols = c("lightgrey", "#D6604D"), pt.size = 1) +
  ggtitle("myoCAF Score")

p2 <- FeaturePlot(stroma_fib, features = "iCAF",
                  cols = c("lightgrey", "#4393C3"), pt.size = 1) +
  ggtitle("iCAF Score")

p3 <- FeaturePlot(stroma_fib, features = "apCAF",
                  cols = c("lightgrey", "#74C476"), pt.size = 1) +
  ggtitle("apCAF Score")

p4 <- DimPlot(stroma_fib, group.by = "subtype", label = TRUE,
              label.size = 3, repel = TRUE, cols = fib_cols) +
  ggtitle("Fibroblast Subtypes") + NoLegend()

(p1 + p2) / (p3 + p4)

CAF Score Violin Plots

v1 <- VlnPlot(stroma_fib, features = "myoCAF", group.by = "subtype",
              pt.size = 0, cols = fib_cols) +
  NoLegend() + RotatedAxis() + ggtitle("myoCAF")

v2 <- VlnPlot(stroma_fib, features = "iCAF", group.by = "subtype",
              pt.size = 0, cols = fib_cols) +
  NoLegend() + RotatedAxis() + ggtitle("iCAF")

v3 <- VlnPlot(stroma_fib, features = "apCAF", group.by = "subtype",
              pt.size = 0, cols = fib_cols) +
  NoLegend() + RotatedAxis() + ggtitle("apCAF")

v1 / v2 / v3

CAF Gene DotPlot

caf_genes_plot <- c(
  # myoCAF
  "Acta2", "Tagln", "Fap", "Postn", "Lox", "Mmp2",
  # iCAF
  "Il6", "Cxcl12", "Cxcl1", "Cxcl2", "Ccl2", "Vcam1", "Il33", "Lcn2",
  # apCAF
  "Cd74", "H2-Ab1", "H2-Eb1", "Slpi"
)
caf_genes_present <- caf_genes_plot[caf_genes_plot %in% rownames(stroma_fib)]

DotPlot(stroma_fib, features = caf_genes_present, group.by = "subtype",
        cols = c("lightgrey", "darkred"), dot.scale = 7) +
  RotatedAxis() +
  ggtitle("CAF Signature Genes by Fibroblast Subtype") +
  theme(axis.text.x = element_text(size = 9))

---

3. Pro-Leukemia Niche Scoring

Interpretation: CNS Meningeal Fibroblasts (ECM) are the highest-scoring population for niche retention (mean 0.309), ahead of Leptomeningeal (0.263) and Mesenchymal/Adventitial (0.244) fibroblasts. Gene-level inspection confirms ECM fibroblasts are the primary expressers of Cxcl12 and Kitl — the two ligands most directly implicated in B-ALL niche retention through CXCR4 and c-Kit signalling respectively. Vcam1 and Fn1 enrichment in ECM fibroblasts supports adhesive interactions with leukaemia cells via VLA-4. Together, these data identify ECM fibroblasts as the dominant source of pro-leukaemic niche signals in the CNS meningeal compartment.

The object already contains Niche_retention_UCell from the prior analysis. This section visualises that score in the context of the fibroblast subtypes and supplements it with key B-ALL ligand expression.

cat("Existing UCell scores in object:\n")

## Existing UCell scores in object:

print(grep("UCell", colnames(stroma@meta.data), value = TRUE))

## [1] "ECM_remodelling_UCell"  "Growth_support_UCell"   "Niche_retention_UCell" 
## [4] "Angiogenic_UCell"       "Immunomodulatory_UCell" "CAF_programme_UCell"

existing_scores <- c("Niche_retention_UCell", "ECM_remodelling_UCell",
                     "Growth_support_UCell", "Angiogenic_UCell",
                     "Immunomodulatory_UCell", "CAF_programme_UCell")

scores_in_fib <- existing_scores[existing_scores %in% colnames(stroma_fib@meta.data)]
cat("\nExisting scores present in fibroblast object:",
    paste(scores_in_fib, collapse = ", "), "\n")

## 
## Existing scores present in fibroblast object: Niche_retention_UCell, ECM_remodelling_UCell, Growth_support_UCell, Angiogenic_UCell, Immunomodulatory_UCell, CAF_programme_UCell

Niche Retention Score UMAP

if ("Niche_retention_UCell" %in% colnames(stroma_fib@meta.data)) {
  p1 <- FeaturePlot(stroma_fib, features = "Niche_retention_UCell",
                    cols = c("lightgrey", "darkred"), pt.size = 1) +
    ggtitle("Niche Retention Score (fibroblasts)")

  p2 <- FeaturePlot(stroma, features = "Niche_retention_UCell",
                    cols = c("lightgrey", "darkred"), pt.size = 0.6) +
    ggtitle("Niche Retention Score (all stroma, context)")

  p1 + p2
}

Niche Retention by Subtype

if ("Niche_retention_UCell" %in% colnames(stroma_fib@meta.data)) {
  VlnPlot(stroma_fib, features = "Niche_retention_UCell",
          group.by = "subtype", pt.size = 0, cols = fib_cols) +
    NoLegend() + RotatedAxis() +
    ggtitle("Niche Retention Score by Fibroblast Subtype")
}

Key B-ALL Ligand Expression

key_ligands <- c("Cxcl12", "Vcam1", "Il6", "Il7", "Kitl",
                 "Spp1", "Tgfb1", "Tgfb2", "Ccl2", "Fn1",
                 "Angpt1", "Lif", "Mif")
key_ligands_present <- key_ligands[key_ligands %in% rownames(stroma_fib)]

DotPlot(stroma_fib, features = key_ligands_present, group.by = "subtype",
        cols = c("lightgrey", "darkred"), dot.scale = 7) +
  RotatedAxis() +
  ggtitle("Pro-Leukemia Ligand Expression by Fibroblast Subtype")

All Existing Scores

if (length(scores_in_fib) > 0) {
  VlnPlot(stroma_fib, features = scores_in_fib, group.by = "subtype",
          pt.size = 0, cols = fib_cols, ncol = 2) &
    NoLegend() & RotatedAxis()
}

---

4. IL-17 Receptor Expression

Interpretation: IL-17A signals through a heterodimeric receptor comprising IL-17RA and IL-17RC, both of which must be co-expressed for functional signalling. Approximately 23% of CNS Meningeal Fibroblasts (ECM) co-express both subunits — substantially higher than Mesenchymal/Adventitial fibroblasts (~11.5%) and markedly higher than Leptomeningeal fibroblasts (~2.5%). The near-absence of IL-17 receptor expression in Leptomeningeal fibroblasts is notable given their anatomical proximity to the subarachnoid space. It suggests that receptor expression — and therefore IL-17A competency — is determined by fibroblast identity rather than simply by spatial proximity to the source of ligand.

IL-17A signals through the IL-17RA/IL-17RC heterodimer. Expression of this receptor complex defines which fibroblast populations are competent to respond to Vγ6Vδ4-derived IL-17A. IL-17RD and IL-17RE are included given prior evidence of expression on immune populations in this dataset.

il17r_genes   <- c("Il17ra", "Il17rc", "Il17rb", "Il17rd", "Il17re")
il17r_present <- il17r_genes[il17r_genes %in% rownames(stroma_fib)]
cat("IL-17 receptor genes found:", paste(il17r_present, collapse = ", "), "\n")

## IL-17 receptor genes found: Il17ra, Il17rc, Il17rb, Il17rd, Il17re

cat("Missing:", paste(setdiff(il17r_genes, rownames(stroma_fib)), collapse = ", "), "\n")

## Missing:

Receptor DotPlot

if (length(il17r_present) > 0) {
  DotPlot(stroma_fib, features = il17r_present, group.by = "subtype",
          cols = c("lightgrey", "#762A83"), dot.scale = 8) +
    RotatedAxis() +
    ggtitle("IL-17 Receptor Expression by Fibroblast Subtype")
} else {
  cat("No IL-17 receptor genes detected in dataset.\n")
}

Receptor FeaturePlots

if (length(il17r_present) > 0) {
  FeaturePlot(stroma_fib, features = il17r_present,
              cols = c("lightgrey", "#762A83"), pt.size = 0.8,
              ncol = min(length(il17r_present), 3))
}

IL-17RA / IL-17RC Co-expression

if (all(c("Il17ra", "Il17rc") %in% rownames(stroma_fib))) {
  expr_mat <- GetAssayData(stroma_fib, layer = "data")[c("Il17ra", "Il17rc"), ]

  stroma_fib$Il17ra_pos   <- expr_mat["Il17ra", ] > 0
  stroma_fib$Il17rc_pos   <- expr_mat["Il17rc", ] > 0
  stroma_fib$Il17r_coexpr <- stroma_fib$Il17ra_pos & stroma_fib$Il17rc_pos

  coexpr_summary <- stroma_fib@meta.data %>%
    group_by(subtype) %>%
    summarise(
      n_cells       = n(),
      pct_Il17ra    = round(mean(Il17ra_pos)   * 100, 1),
      pct_Il17rc    = round(mean(Il17rc_pos)   * 100, 1),
      pct_coexpress = round(mean(Il17r_coexpr) * 100, 1),
      .groups = "drop"
    ) %>%
    arrange(desc(pct_coexpress))

  kable(coexpr_summary,
        caption = "IL-17RA/RC co-expression by fibroblast subtype (% cells)")

  ggplot(coexpr_summary,
         aes(x = reorder(subtype, -pct_coexpress),
             y = pct_coexpress, fill = subtype)) +
    geom_col(alpha = 0.85) +
    scale_fill_manual(values = fib_cols) +
    labs(x = NULL, y = "% cells co-expressing IL-17RA + IL-17RC",
         title = "IL-17 Receptor Co-expression by Fibroblast Subtype") +
    theme_classic() +
    theme(axis.text.x   = element_text(angle = 35, hjust = 1),
          legend.position = "none")
} else {
  cat("Il17ra and/or Il17rc not detected — individual receptor plots only.\n")
}

---

5. IL-17A Downstream Response Scoring

Interpretation: CNS Meningeal Fibroblasts (ECM) have significantly higher IL-17 response scores than the other two populations (mean 0.017 vs 0.010 and 0.003). The dominant gene driving this signal is Mmp3, expressed in approximately 30% of ECM fibroblasts — far exceeding the typical IL-17 outputs of chemokines such as Cxcl1/2. This is an important deviation from the canonical IL-17 response seen in epithelial tissues. Rather than a primarily chemotactic output, ECM fibroblasts appear to respond to IL-17A through a matrix-remodelling programme centred on Mmp3 (and to a lesser extent Mmp13). This distinction likely reflects the fibroblastic identity of the responding cells, and has direct implications for how IL-17A acts in the leukaemic niche.

Canonical NF-κB/Act1 downstream outputs of IL-17A signalling in stromal cells. High scores in IL-17 receptor-expressing clusters indicate active IL-17A responsiveness, supporting the Vγ6Vδ4 → IL-17A → fibroblast activation axis.

il17_response_genes <- c(
  # CXC chemokines — neutrophil and immune cell recruitment
  "Cxcl1", "Cxcl2", "Cxcl5",
  # CC chemokines
  "Ccl2", "Ccl20",
  # Cytokines
  "Il6", "Il1b", "Tnf",
  # Antimicrobial / acute phase
  "Lcn2", "S100a8", "S100a9",
  # Matrix metalloproteinases
  "Mmp3", "Mmp10", "Mmp13",
  # Granulopoiesis
  "Csf3", "Csf2",
  # Hyaluronan synthase (iCAF marker and IL-17 target)
  "Has1", "Has2"
)

il17_resp_present <- il17_response_genes[il17_response_genes %in% rownames(stroma_fib)]
cat("IL-17 response genes found:", length(il17_resp_present), "/",
    length(il17_response_genes), "\n")

## IL-17 response genes found: 18 / 18

cat("Missing:", paste(setdiff(il17_response_genes, rownames(stroma_fib)), collapse = ", "), "\n")

## Missing:

stroma_fib <- AddModuleScore_UCell(
  stroma_fib,
  features = list(IL17_response = il17_resp_present),
  name = ""
)

Response Score UMAP

p1 <- FeaturePlot(stroma_fib, features = "IL17_response",
                  cols = c("lightgrey", "#762A83"), pt.size = 1) +
  ggtitle("IL-17 Response Score")

p2 <- DimPlot(stroma_fib, group.by = "subtype", label = TRUE,
              label.size = 3, repel = TRUE, cols = fib_cols) +
  ggtitle("Subtypes (reference)") + NoLegend()

p1 + p2

Response Score by Subtype

VlnPlot(stroma_fib, features = "IL17_response", group.by = "subtype",
        pt.size = 0, cols = fib_cols) +
  NoLegend() + RotatedAxis() +
  ggtitle("IL-17 Response Score by Fibroblast Subtype")

Response Gene DotPlot

DotPlot(stroma_fib, features = il17_resp_present, group.by = "subtype",
        cols = c("lightgrey", "#762A83"), dot.scale = 6) +
  RotatedAxis() +
  ggtitle("IL-17 Response Gene Expression by Fibroblast Subtype")

---

6. Integration: iCAF vs IL-17 Response

Interpretation: Cell-level plotting of iCAF score against IL-17 response score confirms that in ECM fibroblasts, cells with higher IL-17 response scores do not necessarily show higher iCAF scores — the two are partially dissociated. This is consistent with the finding that the IL-17 response in these cells is MMP-driven rather than cytokine-driven, arguing against simple conflation of IL-17 responsiveness with inflammatory CAF activation. The combined score heatmap positions ECM fibroblasts as transcriptionally distinct from the other two populations across all scored axes simultaneously, reinforcing their identity as the most functionally active fibroblast population in the leukaemic CNS niche.

If iCAF score and IL-17 response score co-localise to the same fibroblast populations, this supports the model that Vγ6Vδ4-derived IL-17A drives iCAF-like activation in the CNS meningeal stroma.

# Build score summary — include existing scores if present
score_summary <- stroma_fib@meta.data %>%
  group_by(subtype) %>%
  summarise(
    myoCAF        = mean(myoCAF),
    iCAF          = mean(iCAF),
    apCAF         = mean(apCAF),
    IL17_response = mean(IL17_response),
    Niche_retention = if ("Niche_retention_UCell" %in% names(pick(everything())))
                         mean(Niche_retention_UCell) else NA_real_,
    n             = n(),
    .groups = "drop"
  )

kable(score_summary %>% mutate(across(where(is.numeric) & !n, ~round(., 3))),
      caption = "Mean scores by fibroblast subtype")

Mean scores by fibroblast subtype

subtype	myoCAF	iCAF	apCAF	IL17_response	Niche_retention	n
CNS Leptomeningeal Fibroblasts	0.102	0.108	0.114	0.010	0.263	294
CNS Meningeal Fibroblasts (ECM)	0.177	0.126	0.050	0.017	0.309	275
CNS Mesenchymal/Adventitial	0.116	0.102	0.129	0.003	0.244	122

iCAF vs IL-17 Response: Cell Level

ggplot(stroma_fib@meta.data,
       aes(x = iCAF, y = IL17_response, colour = subtype)) +
  geom_point(alpha = 0.4, size = 0.8) +
  geom_smooth(method = "lm", se = TRUE, linewidth = 0.9) +
  scale_colour_manual(values = fib_cols) +
  labs(x = "iCAF Score", y = "IL-17 Response Score",
       title = "iCAF vs IL-17 Response: Cell-Level",
       colour = NULL) +
  theme_classic() +
  theme(legend.position = "bottom",
        legend.text = element_text(size = 8))

iCAF vs IL-17 Response: Subtype Means

ggplot(score_summary,
       aes(x = iCAF, y = IL17_response, colour = subtype,
           size = n, label = subtype)) +
  geom_point(alpha = 0.9) +
  geom_text_repel(size = 3.5, max.overlaps = 10) +
  scale_colour_manual(values = fib_cols) +
  scale_size_continuous(range = c(4, 10)) +
  labs(x = "Mean iCAF Score", y = "Mean IL-17 Response Score",
       title = "iCAF vs IL-17 Response: Subtype Means",
       subtitle = "Size = cell count") +
  theme_classic() +
  theme(legend.position = "none")

Combined Score Heatmap

hm_base_cols <- c("myoCAF", "iCAF", "apCAF", "IL17_response")
optional_cols <- c("Niche_retention_UCell", "ECM_remodelling_UCell",
                   "Angiogenic_UCell", "Immunomodulatory_UCell")
hm_cols <- c(hm_base_cols,
             optional_cols[optional_cols %in% colnames(stroma_fib@meta.data)])

hm_data <- stroma_fib@meta.data %>%
  group_by(subtype) %>%
  summarise(across(all_of(hm_cols), mean), .groups = "drop") %>%
  column_to_rownames("subtype") %>%
  as.matrix()

hm_scaled <- scale(hm_data)

heatmap(hm_scaled,
        col     = colorRampPalette(c("#2166AC", "white", "#B2182B"))(50),
        margins = c(14, 22),
        main    = "Stromal Scores (z-scaled) by Fibroblast Subtype",
        cexRow  = 0.9, cexCol  = 0.85)

---

7. Inter-Subtype Comparison: CNS Fibroblasts

Interpretation: Formal statistical comparison (Kruskal-Wallis with pairwise Wilcoxon post-hoc tests) confirms that score differences between subtypes are significant across myoCAF, niche retention, and IL-17 response axes. ECM fibroblasts are statistically distinguishable from both other populations on these scores. Leptomeningeal and Mesenchymal/Adventitial fibroblasts are more similar to each other than either is to ECM fibroblasts, particularly on myoCAF and niche retention axes.

All three fibroblast populations are CNS-restricted, so the relevant comparison is between subtypes rather than across tissues. This section tests whether iCAF character and IL-17 responsiveness differ significantly between Leptomeningeal, Meningeal (ECM), and Mesenchymal/Adventitial populations.

scores_to_test <- c("myoCAF", "iCAF", "apCAF", "IL17_response")
scores_present <- scores_to_test[scores_to_test %in% colnames(stroma_fib@meta.data)]

kw_results <- lapply(scores_present, function(sc) {
  vals <- stroma_fib@meta.data[[sc]]
  grps <- stroma_fib$subtype
  kw   <- kruskal.test(vals ~ grps)
  data.frame(
    Score   = sc,
    H       = round(kw$statistic, 2),
    df      = kw$parameter,
    p_value = signif(kw$p.value, 3)
  )
}) %>% bind_rows()

kw_results$p_adj <- p.adjust(kw_results$p_value, method = "BH")
kw_results$sig   <- case_when(
  kw_results$p_adj < 0.001 ~ "***",
  kw_results$p_adj < 0.01  ~ "**",
  kw_results$p_adj < 0.05  ~ "*",
  TRUE                     ~ "ns"
)

kable(kw_results, caption = "Kruskal-Wallis test across CNS fibroblast subtypes")

Kruskal-Wallis test across CNS fibroblast subtypes

	Score	H	df	p_value	p_adj	sig
Kruskal-Wallis chi-squared…1	myoCAF	154.78	2	0.0e+00	0.0e+00	***
Kruskal-Wallis chi-squared…2	iCAF	23.11	2	9.6e-06	1.1e-05	***
Kruskal-Wallis chi-squared…3	apCAF	74.53	2	0.0e+00	0.0e+00	***
Kruskal-Wallis chi-squared…4	IL17_response	22.84	2	1.1e-05	1.1e-05	***

sig_scores <- kw_results %>% filter(p_adj < 0.05) %>% pull(Score)

if (length(sig_scores) > 0) {
  pw_results <- lapply(sig_scores, function(sc) {
    pw <- pairwise.wilcox.test(
      stroma_fib@meta.data[[sc]],
      stroma_fib$subtype,
      p.adjust.method = "BH"
    )
    as.data.frame(pw$p.value) %>%
      rownames_to_column("Group1") %>%
      pivot_longer(-Group1, names_to = "Group2", values_to = "p_adj") %>%
      filter(!is.na(p_adj)) %>%
      mutate(Score = sc,
             p_adj = signif(p_adj, 3),
             sig   = case_when(p_adj < 0.001 ~ "***",
                               p_adj < 0.01  ~ "**",
                               p_adj < 0.05  ~ "*",
                               TRUE          ~ "ns"))
  }) %>% bind_rows()

  kable(pw_results %>% arrange(Score, p_adj),
        caption = "Pairwise Wilcoxon (BH-adjusted) for significant scores")
} else {
  cat("No scores significant after BH correction.\n")
}

Pairwise Wilcoxon (BH-adjusted) for significant scores

Group1	Group2	p_adj	Score	sig
CNS Mesenchymal/Adventitial	CNS Meningeal Fibroblasts (ECM)	1.12e-05	IL17_response	***
CNS Mesenchymal/Adventitial	CNS Leptomeningeal Fibroblasts	1.86e-03	IL17_response	**
CNS Meningeal Fibroblasts (ECM)	CNS Leptomeningeal Fibroblasts	1.77e-02	IL17_response	*
CNS Meningeal Fibroblasts (ECM)	CNS Leptomeningeal Fibroblasts	0.00e+00	apCAF	***
CNS Mesenchymal/Adventitial	CNS Meningeal Fibroblasts (ECM)	0.00e+00	apCAF	***
CNS Mesenchymal/Adventitial	CNS Leptomeningeal Fibroblasts	4.78e-01	apCAF	ns
CNS Mesenchymal/Adventitial	CNS Meningeal Fibroblasts (ECM)	4.13e-05	iCAF	***
CNS Meningeal Fibroblasts (ECM)	CNS Leptomeningeal Fibroblasts	4.44e-04	iCAF	***
CNS Mesenchymal/Adventitial	CNS Leptomeningeal Fibroblasts	1.21e-01	iCAF	ns
CNS Meningeal Fibroblasts (ECM)	CNS Leptomeningeal Fibroblasts	0.00e+00	myoCAF	***
CNS Mesenchymal/Adventitial	CNS Meningeal Fibroblasts (ECM)	0.00e+00	myoCAF	***
CNS Mesenchymal/Adventitial	CNS Leptomeningeal Fibroblasts	1.47e-02	myoCAF	*

Score Box Plots

plot_list <- lapply(scores_present, function(sc) {
  ggplot(stroma_fib@meta.data,
         aes(x = subtype, y = .data[[sc]], fill = subtype)) +
    geom_boxplot(outlier.size = 0.4, alpha = 0.85) +
    scale_fill_manual(values = fib_cols) +
    labs(x = NULL, y = NULL, title = sc) +
    theme_classic() +
    theme(axis.text.x   = element_text(angle = 35, hjust = 1, size = 8),
          legend.position = "none")
})

wrap_plots(plot_list, ncol = 2)

iCAF vs IL-17 Response Scatter

ggplot(stroma_fib@meta.data,
       aes(x = iCAF, y = IL17_response, colour = subtype)) +
  geom_point(alpha = 0.5, size = 0.8) +
  geom_smooth(method = "lm", se = TRUE, linewidth = 0.8) +
  scale_colour_manual(values = fib_cols) +
  labs(x = "iCAF Score", y = "IL-17 Response Score",
       title = "iCAF vs IL-17 Response: Cell-Level by Subtype",
       colour = NULL) +
  theme_classic() +
  theme(legend.position = "bottom",
        legend.text = element_text(size = 8))

---

8. IL-17-High ECM Fibroblast Subset

Interpretation: Within ECM fibroblasts, the distribution of IL-17 response scores is bimodal. Defining the top-quartile as IL17-high and comparing to the majority reveals that the defining transcriptional differences are almost entirely accounted for by Mmp3 and Mmp13 — no significant chemokine or cytokine upregulation is detectable in IL17-high cells. The UMAP localisation of IL17-high cells is spatially coherent, concentrated in the upper-right region of the ECM cluster, consistent with a genuine sub-state rather than stochastic variation. Genes depressed in IL17-high cells (Slc38a2, Slc47a1, Slc23a2 — amino acid and vitamin C transporters) suggest the shift toward an MMP-high state involves partial loss of a homeostatic metabolic programme. Co-expression analysis between Mmp3 status and niche ligand expression is more nuanced than a simple enrichment story. Cxcl12 is near-universally expressed across ECM fibroblasts regardless of Mmp3 status (100% vs 98.4%), and Kitl is actually slightly higher in Mmp3- cells (76.9% vs 87.5%). Vcam1 shows the most genuine enrichment in Mmp3+ cells (75.8% vs 60.9%), and Il6 shows a modest fold-change though at low absolute frequencies. The important implication is that Cxcl12 and Kitl are constitutive features of ECM fibroblast identity rather than being specifically linked to the IL-17-responding subset. Mmp3 may therefore act on constitutively expressed niche ligands by liberating matrix-bound forms — releasing sequestered Cxcl12 and membrane-associated Kitl through proteolytic ECM remodelling — rather than by transcriptionally upregulating their expression in the same cells.

The bimodal IL-17 response score distribution in CNS Meningeal Fibroblasts (ECM) suggests a transcriptionally distinct subset is driving the signal. This section defines high vs low responders within that cluster and asks what distinguishes them.

# Work within ECM fibroblasts only
ecm_fib <- subset(stroma_fib,
                  subset = subtype == "CNS Meningeal Fibroblasts (ECM)")

cat("ECM fibroblasts:", ncol(ecm_fib), "cells\n")

## ECM fibroblasts: 275 cells

# Define threshold — top quartile of IL17_response as "high"
thresh <- quantile(ecm_fib$IL17_response, 0.75)
cat("IL17_response 75th percentile threshold:", round(thresh, 4), "\n")

## IL17_response 75th percentile threshold: 0.0275

ecm_fib$IL17_group <- ifelse(ecm_fib$IL17_response >= thresh,
                              "IL17_high", "IL17_low")

cat("\nIL17 group counts:\n")

## 
## IL17 group counts:

print(table(ecm_fib$IL17_group))

## 
## IL17_high  IL17_low 
##        69       206

# Propagate back to stroma_fib so Section 9 can filter from its metadata
stroma_fib$IL17_group <- NA_character_
stroma_fib$IL17_group[colnames(ecm_fib)] <- ecm_fib$IL17_group

Score Distribution with Threshold

ggplot(ecm_fib@meta.data, aes(x = IL17_response, fill = IL17_group)) +
  geom_histogram(bins = 40, alpha = 0.85, colour = "white") +
  geom_vline(xintercept = thresh, linetype = "dashed", colour = "black") +
  scale_fill_manual(values = c("IL17_high" = "#762A83", "IL17_low" = "grey70")) +
  labs(x = "IL-17 Response Score", y = "Cell count",
       title = "IL-17 Response Score Distribution — ECM Fibroblasts",
       subtitle = paste0("Threshold (Q75) = ", round(thresh, 4)),
       fill = NULL) +
  theme_classic()

UMAP: IL17-High vs Low

DimPlot(ecm_fib, group.by = "IL17_group",
        cols = c("IL17_high" = "#762A83", "IL17_low" = "grey80"),
        pt.size = 1.2) +
  ggtitle("IL-17-High vs Low Cells within ECM Fibroblasts")

Co-score Comparison: High vs Low

score_cols <- c("myoCAF", "iCAF", "apCAF", "IL17_response",
                "Niche_retention_UCell")[
  c("myoCAF", "iCAF", "apCAF", "IL17_response",
    "Niche_retention_UCell") %in% colnames(ecm_fib@meta.data)]

score_long <- ecm_fib@meta.data %>%
  select(IL17_group, all_of(score_cols)) %>%
  pivot_longer(-IL17_group, names_to = "Score", values_to = "Value")

ggplot(score_long, aes(x = IL17_group, y = Value, fill = IL17_group)) +
  geom_boxplot(outlier.size = 0.4, alpha = 0.85) +
  scale_fill_manual(values = c("IL17_high" = "#762A83", "IL17_low" = "grey70")) +
  facet_wrap(~Score, scales = "free_y", ncol = 3) +
  labs(x = NULL, y = "Score", title = "Score Comparison: IL-17-High vs Low ECM Fibroblasts") +
  theme_classic() +
  theme(legend.position = "none",
        strip.text = element_text(size = 9))

Differential Expression: High vs Low

Idents(ecm_fib) <- "IL17_group"

# Require ≥10 cells per group — already confirmed above
de_il17 <- FindMarkers(ecm_fib,
                       ident.1    = "IL17_high",
                       ident.2    = "IL17_low",
                       test.use   = "wilcox",
                       min.pct    = 0.05,
                       logfc.threshold = 0.2,
                       verbose    = FALSE)

de_il17 <- de_il17 %>%
  rownames_to_column("gene") %>%
  arrange(p_val_adj, desc(abs(avg_log2FC)))

cat("DE genes (IL17_high vs IL17_low):", nrow(de_il17), "\n")

## DE genes (IL17_high vs IL17_low): 7346

cat("Significant (p_adj < 0.05):", sum(de_il17$p_val_adj < 0.05, na.rm = TRUE), "\n\n")

## Significant (p_adj < 0.05): 27

de_il17 %>%
  filter(p_val_adj < 0.05) %>%
  slice_head(n = 30) %>%
  mutate(across(where(is.numeric), ~round(., 4))) %>%
  kable(caption = "Top DE genes: IL-17-high vs IL-17-low ECM fibroblasts")

Top DE genes: IL-17-high vs IL-17-low ECM fibroblasts

gene	p_val	avg_log2FC	pct.1	pct.2	p_val_adj
Mmp3	0	7.4275	0.826	0.165	0.0000
Slc38a2	0	-0.9033	0.986	0.995	0.0000
Slc47a1	0	-1.3987	0.696	0.942	0.0000
Slc16a3	0	2.3833	0.406	0.092	0.0001
Prmt8	0	-1.2537	0.725	0.883	0.0001
Slc23a2	0	-1.2446	0.652	0.879	0.0001
Ifitm2	0	0.4479	1.000	0.995	0.0014
Foxd1	0	-0.7307	0.942	0.961	0.0042
Fxyd5	0	-0.5142	1.000	0.995	0.0050
Ccl19-ps2	0	1.1676	0.812	0.524	0.0053
Mmp13	0	7.3005	0.145	0.005	0.0075
Ctsl	0	0.7835	1.000	0.990	0.0083
Mif	0	1.0697	1.000	0.956	0.0086
C3	0	1.2482	0.855	0.709	0.0116
Tspan11	0	-0.5901	0.986	0.995	0.0118
Ddx5	0	-0.4371	1.000	0.995	0.0129
Pkm	0	0.6526	1.000	0.956	0.0155
Cp	0	0.6103	1.000	0.990	0.0161
Lmo4	0	-0.6691	0.913	0.951	0.0210
Sdc4	0	1.3391	0.928	0.879	0.0265
Slc4a10	0	-0.6289	0.971	0.990	0.0273
Add3	0	-0.6706	0.928	0.966	0.0299
Gpx3	0	0.8479	0.957	0.951	0.0300
Six2	0	-1.3614	0.348	0.636	0.0313
Fam124a	0	2.2244	0.304	0.078	0.0331
Ecrg4	0	-0.7766	0.942	0.981	0.0385
Ptn	0	-0.7835	0.768	0.932	0.0463

Volcano Plot

de_il17_plot <- de_il17 %>%
  mutate(
    sig      = p_val_adj < 0.05 & abs(avg_log2FC) > 0.2,
    direction = case_when(
      sig & avg_log2FC > 0 ~ "Up in IL17-high",
      sig & avg_log2FC < 0 ~ "Up in IL17-low",
      TRUE                 ~ "ns"
    ),
    label = ifelse(sig & abs(avg_log2FC) > 0.5, gene, NA)
  )

ggplot(de_il17_plot,
       aes(x = avg_log2FC, y = -log10(p_val_adj + 1e-300),
           colour = direction, label = label)) +
  geom_point(alpha = 0.6, size = 0.9) +
  geom_text_repel(size = 3, max.overlaps = 20, na.rm = TRUE) +
  scale_colour_manual(values = c("Up in IL17-high" = "#762A83",
                                 "Up in IL17-low"  = "#4393C3",
                                 "ns"              = "grey70")) +
  geom_vline(xintercept = c(-0.2, 0.2), linetype = "dashed", alpha = 0.5) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", alpha = 0.5) +
  labs(x = "log2 Fold Change (IL17-high vs IL17-low)",
       y = "-log10(adjusted p-value)",
       title = "DE: IL-17-High vs IL-17-Low ECM Fibroblasts",
       colour = NULL) +
  theme_classic()

Mmp3/Mmp13 Co-expression and Niche Ligand Context

The volcano identifies Mmp3 and Mmp13 as the dominant outputs of IL-17 signalling in ECM fibroblasts. This tab asks whether the same cells co-express canonical niche retention ligands (Kitl, Cxcl12), directly testing the hypothesis that MMP-mediated ECM remodelling and growth factor presentation occur in the same cells.

# Genes of interest
mmp_genes   <- c("Mmp3", "Mmp13")
niche_genes <- c("Cxcl12", "Kitl", "Vcam1", "Il6", "Tgfb1")
query_genes <- c(mmp_genes, niche_genes)
query_genes <- query_genes[query_genes %in% rownames(ecm_fib)]

# Extract normalised expression
expr <- GetAssayData(ecm_fib, layer = "data")[query_genes, , drop = FALSE]

# Binary expression flags
mmp3_pos  <- expr["Mmp3",  ] > 0
mmp13_pos <- if ("Mmp13" %in% rownames(expr)) expr["Mmp13", ] > 0 else rep(FALSE, ncol(ecm_fib))

ecm_fib$Mmp3_pos   <- as.logical(mmp3_pos)
ecm_fib$Mmp13_pos  <- as.logical(mmp13_pos)
ecm_fib$Mmp_either <- ecm_fib$Mmp3_pos | ecm_fib$Mmp13_pos
ecm_fib$Mmp_both   <- ecm_fib$Mmp3_pos & ecm_fib$Mmp13_pos

cat("ECM fibroblast MMP expression:\n")

## ECM fibroblast MMP expression:

cat("  Mmp3+:            ", sum(ecm_fib$Mmp3_pos),  "/", ncol(ecm_fib),
    "(", round(mean(ecm_fib$Mmp3_pos)*100, 1), "%)\n")

##   Mmp3+:             91 / 275 ( 33.1 %)

cat("  Mmp13+:           ", sum(ecm_fib$Mmp13_pos), "/", ncol(ecm_fib),
    "(", round(mean(ecm_fib$Mmp13_pos)*100, 1), "%)\n")

##   Mmp13+:            11 / 275 ( 4 %)

cat("  Mmp3+ OR Mmp13+:  ", sum(ecm_fib$Mmp_either),"/", ncol(ecm_fib),
    "(", round(mean(ecm_fib$Mmp_either)*100, 1), "%)\n")

##   Mmp3+ OR Mmp13+:   96 / 275 ( 34.9 %)

cat("  Mmp3+ AND Mmp13+: ", sum(ecm_fib$Mmp_both),  "/", ncol(ecm_fib),
    "(", round(mean(ecm_fib$Mmp_both)*100, 1), "%)\n\n")

##   Mmp3+ AND Mmp13+:  6 / 275 ( 2.2 %)

# Co-expression with niche ligands in Mmp3+ vs Mmp3- cells
coexpr_tbl <- data.frame(
  gene    = niche_genes[niche_genes %in% rownames(expr)],
  pct_mmp3pos  = sapply(niche_genes[niche_genes %in% rownames(expr)], function(g)
    round(mean(expr[g, ecm_fib$Mmp3_pos] > 0) * 100, 1)),
  pct_mmp3neg  = sapply(niche_genes[niche_genes %in% rownames(expr)], function(g)
    round(mean(expr[g, !ecm_fib$Mmp3_pos] > 0) * 100, 1))
)
colnames(coexpr_tbl) <- c("Gene", "% expressing (Mmp3+)", "% expressing (Mmp3-)")

cat("Niche ligand co-expression with Mmp3:\n")

## Niche ligand co-expression with Mmp3:

print(coexpr_tbl, row.names = FALSE)

##    Gene % expressing (Mmp3+) % expressing (Mmp3-)
##  Cxcl12                100.0                 98.4
##    Kitl                 76.9                 87.5
##   Vcam1                 75.8                 60.9
##     Il6                  7.7                  2.7
##   Tgfb1                 56.0                 51.6

# DotPlot split by Mmp3 status
ecm_fib$Mmp3_group <- ifelse(ecm_fib$Mmp3_pos, "Mmp3+", "Mmp3-")

DotPlot(ecm_fib,
        features   = query_genes,
        group.by   = "Mmp3_group",
        cols       = c("lightgrey", "#762A83"),
        dot.scale  = 8) +
  RotatedAxis() +
  labs(title = "MMP and Niche Ligand Expression: Mmp3+ vs Mmp3- ECM Fibroblasts",
       x = NULL, y = NULL)

# UMAP: Mmp3 status | Mmp3 expression | Cxcl12 expression | Kitl expression
p_group  <- DimPlot(ecm_fib, group.by = "Mmp3_group",
                    cols = c("Mmp3+" = "#762A83", "Mmp3-" = "grey80"),
                    pt.size = 1.2, order = "Mmp3+") +
  ggtitle("Mmp3 Expression Status")

p_mmp3   <- FeaturePlot(ecm_fib, features = "Mmp3",
                         cols = c("grey90", "#762A83"), pt.size = 1) +
  ggtitle("Mmp3")

p_cxcl12 <- FeaturePlot(ecm_fib, features = "Cxcl12",
                          cols = c("grey90", "#1B7837"), pt.size = 1) +
  ggtitle("Cxcl12")

p_kitl   <- FeaturePlot(ecm_fib, features = "Kitl",
                          cols = c("grey90", "#D6604D"), pt.size = 1) +
  ggtitle("Kitl")

p_group + p_mmp3 + p_cxcl12 + p_kitl + plot_layout(ncol = 4)

---

9. IL-17 Response vs Niche Retention: Cell-Level Correlation

Interpretation: Cell-level correlation between IL-17 response score and niche retention score within ECM fibroblasts yields a significant positive association (Pearson r = 0.163, p = 0.0066). The modest effect size is expected — niche retention reflects constitutive ligand expression high across ECM fibroblasts generally, while IL-17 response captures a dynamic, stimulus-dependent layer on top of that baseline. The directional signal is clear: cells with greater IL-17A transcriptional output are measurably more pro-leukaemic. IL17-high cells (purple) are distributed across the middle and upper range of niche retention scores, with few falling below ~0.2 — suggesting the responding subset is drawn preferentially from the already niche-supportive portion of the ECM fibroblast population. Together with Section 10, this closes the loop between γδ T cell-derived IL-17A and stromal pro-leukaemic function at single-cell resolution.

Within ECM fibroblasts, does the IL-17-responding subset also show higher niche retention capacity? If yes, this directly links the IL-17 signalling axis to pro-leukemic function in the same cells.

if (!"Niche_retention_UCell" %in% colnames(stroma_fib@meta.data)) {
  cat("Niche_retention_UCell not present — skipping section.\n")
  knitr::knit_exit()
}

Cell-Level Scatter: ECM Fibroblasts

ecm_data <- ecm_fib@meta.data

# Pearson correlation
cor_ecm <- cor.test(ecm_data$IL17_response,
                    ecm_data$Niche_retention_UCell,
                    method = "pearson")

ggplot(ecm_data,
       aes(x = IL17_response, y = Niche_retention_UCell,
           colour = IL17_group)) +
  geom_point(alpha = 0.5, size = 0.9) +
  geom_smooth(method = "lm", se = TRUE, colour = "black", linewidth = 0.8) +
  scale_colour_manual(values = c("IL17_high" = "#762A83", "IL17_low" = "grey60")) +
  annotate("text", x = Inf, y = Inf, hjust = 1.1, vjust = 1.5,
           label = paste0("r = ", round(cor_ecm$estimate, 3),
                          "\np = ", signif(cor_ecm$p.value, 3)),
           size = 4) +
  labs(x = "IL-17 Response Score",
       y = "Niche Retention Score",
       title = "IL-17 Response vs Niche Retention — ECM Fibroblasts",
       colour = "IL-17 Group") +
  theme_classic()

Cell-Level Scatter: All Fibroblasts

ggplot(stroma_fib@meta.data,
       aes(x = IL17_response, y = Niche_retention_UCell,
           colour = subtype)) +
  geom_point(alpha = 0.4, size = 0.8) +
  geom_smooth(method = "lm", se = TRUE, linewidth = 0.8) +
  scale_colour_manual(values = fib_cols) +
  labs(x = "IL-17 Response Score",
       y = "Niche Retention Score",
       title = "IL-17 Response vs Niche Retention — All Fibroblasts",
       colour = NULL) +
  theme_classic() +
  theme(legend.position = "bottom",
        legend.text = element_text(size = 8))

Niche Retention: High vs Low IL-17 Responders

# Restrict to ECM fibroblasts
wt <- wilcox.test(
  ecm_fib@meta.data %>% filter(IL17_group == "IL17_high") %>% pull(Niche_retention_UCell),
  ecm_fib@meta.data %>% filter(IL17_group == "IL17_low")  %>% pull(Niche_retention_UCell)
)

cat("Wilcoxon test — Niche Retention: IL17-high vs IL17-low ECM fibroblasts\n")

## Wilcoxon test — Niche Retention: IL17-high vs IL17-low ECM fibroblasts

cat("W =", wt$statistic, "| p =", signif(wt$p.value, 3), "\n")

## W = 7595 | p = 0.394

ggplot(ecm_fib@meta.data,
       aes(x = IL17_group, y = Niche_retention_UCell, fill = IL17_group)) +
  geom_boxplot(outlier.size = 0.5, alpha = 0.85) +
  scale_fill_manual(values = c("IL17_high" = "#762A83", "IL17_low" = "grey70")) +
  annotate("text", x = 1.5, y = max(ecm_fib$Niche_retention_UCell) * 1.05,
           label = paste0("p = ", signif(wt$p.value, 3)), size = 4) +
  labs(x = NULL, y = "Niche Retention Score",
       title = "Niche Retention: IL-17-High vs Low ECM Fibroblasts") +
  theme_classic() +
  theme(legend.position = "none")

---

10. IL-17 Receptor-Positive Cell Characterisation

Interpretation: IL17R+ cells are concentrated in the ECM and Mesenchymal/Adventitial compartments and essentially absent from Leptomeningeal fibroblasts, consistent with co-expression percentages in Section 4. The DotPlot comparing IL17R+ and IL17R- ECM fibroblasts reveals that receptor-positive cells express substantially higher levels of Cxcl12, Vcam1, Kitl, Fn1, Tgfb2, and Col1a1. The fibroblasts structurally equipped to receive IL-17A signalling are already the most niche-supportive cells at baseline — IL-17A therefore does not convert a neutral fibroblast into a pro-leukaemic one, but amplifies an existing pro-leukaemic state in a pre-committed subpopulation. Differential expression between IL17R+ and IL17R- ECM fibroblasts returns only three significant genes: Il17rc, Il17ra (expected), and Ppp1r13b (iASPP) — a p53 inhibitor and pro-survival regulator. Its enrichment in IL17R+ cells suggests this population may be intrinsically more resistant to apoptotic stress, a feature that would help maintain the niche-supportive fibroblast population in the inflammatory leukaemic CNS environment. The absence of broader transcriptional divergence between R+ and R- cells is itself informative: receptor expression is not associated with wholesale transcriptional reprogramming. The key distinction is competency to respond when IL-17A is present.

Rather than thresholding on the response score, this section identifies cells that express the functional IL-17RA/RC receptor heterodimer and asks what they look like transcriptionally — providing a receptor-based view of IL-17 competency independent of downstream score.

# Requires Il17ra and Il17rc to be in the dataset
if (!all(c("Il17ra", "Il17rc") %in% rownames(stroma_fib))) {
  cat("Il17ra and/or Il17rc not detected — skipping section.\n")
  knitr::knit_exit()
}

# Cells already labelled from Section 4
if (!"Il17r_coexpr" %in% colnames(stroma_fib@meta.data)) {
  expr_mat <- GetAssayData(stroma_fib, layer = "data")[c("Il17ra", "Il17rc"), ]
  stroma_fib$Il17ra_pos   <- expr_mat["Il17ra", ] > 0
  stroma_fib$Il17rc_pos   <- expr_mat["Il17rc", ] > 0
  stroma_fib$Il17r_coexpr <- stroma_fib$Il17ra_pos & stroma_fib$Il17rc_pos
}

stroma_fib$IL17R_group <- ifelse(stroma_fib$Il17r_coexpr,
                                  "IL17R_positive", "IL17R_negative")

cat("IL-17R+ cells (RA+RC co-expressing):\n")

## IL-17R+ cells (RA+RC co-expressing):

print(table(stroma_fib$IL17R_group, stroma_fib$subtype))

##                 
##                  CNS Leptomeningeal Fibroblasts CNS Meningeal Fibroblasts (ECM)
##   IL17R_negative                            287                             210
##   IL17R_positive                              7                              65
##                 
##                  CNS Mesenchymal/Adventitial
##   IL17R_negative                         108
##   IL17R_positive                          14

UMAP: Receptor-Positive Cells

p1 <- DimPlot(stroma_fib, group.by = "IL17R_group",
              cols = c("IL17R_positive" = "#762A83", "IL17R_negative" = "grey85"),
              pt.size = 1, order = "IL17R_positive") +
  ggtitle("IL-17RA/RC Co-expressing Cells")

p2 <- DimPlot(stroma_fib, group.by = "subtype",
              cols = fib_cols, pt.size = 0.8, label = TRUE,
              label.size = 3, repel = TRUE) +
  ggtitle("Subtypes (reference)") + NoLegend()

p1 + p2

Score Comparison: Receptor+ vs Receptor-

score_cols_r <- c("myoCAF", "iCAF", "apCAF", "IL17_response",
                  "Niche_retention_UCell")[
  c("myoCAF", "iCAF", "apCAF", "IL17_response",
    "Niche_retention_UCell") %in% colnames(stroma_fib@meta.data)]

score_long_r <- stroma_fib@meta.data %>%
  select(IL17R_group, subtype, all_of(score_cols_r)) %>%
  pivot_longer(all_of(score_cols_r), names_to = "Score", values_to = "Value")

ggplot(score_long_r,
       aes(x = IL17R_group, y = Value, fill = IL17R_group)) +
  geom_boxplot(outlier.size = 0.3, alpha = 0.85) +
  scale_fill_manual(values = c("IL17R_positive" = "#762A83",
                                "IL17R_negative" = "grey70")) +
  facet_grid(Score ~ subtype, scales = "free_y") +
  labs(x = NULL, y = "Score",
       title = "Score Comparison: IL-17R+ vs IL-17R- Cells by Subtype") +
  theme_classic() +
  theme(legend.position  = "none",
        axis.text.x      = element_text(angle = 35, hjust = 1, size = 7),
        strip.text       = element_text(size = 8))

DE: Receptor+ vs Receptor- (ECM Fibroblasts Only)

ecm_fib$IL17R_group <- stroma_fib$IL17R_group[colnames(ecm_fib)]

n_pos <- sum(ecm_fib$IL17R_group == "IL17R_positive", na.rm = TRUE)
n_neg <- sum(ecm_fib$IL17R_group == "IL17R_negative", na.rm = TRUE)
cat("ECM fibroblasts — IL17R+:", n_pos, "| IL17R-:", n_neg, "\n")

## ECM fibroblasts — IL17R+: 65 | IL17R-: 210

if (n_pos >= 10 & n_neg >= 10) {
  Idents(ecm_fib) <- "IL17R_group"

  de_il17r <- FindMarkers(ecm_fib,
                          ident.1          = "IL17R_positive",
                          ident.2          = "IL17R_negative",
                          test.use         = "wilcox",
                          min.pct          = 0.05,
                          logfc.threshold  = 0.2,
                          verbose          = FALSE) %>%
    rownames_to_column("gene") %>%
    arrange(p_val_adj, desc(abs(avg_log2FC)))

  cat("DE genes:", nrow(de_il17r), "\n")
  cat("Significant (p_adj < 0.05):", sum(de_il17r$p_val_adj < 0.05, na.rm = TRUE), "\n\n")

  de_il17r %>%
    filter(p_val_adj < 0.05) %>%
    slice_head(n = 30) %>%
    mutate(across(where(is.numeric), ~round(., 4))) %>%
    kable(caption = "Top DE genes: IL-17R+ vs IL-17R- ECM fibroblasts")
} else {
  cat("Insufficient cells in one group for DE — skipping.\n")
}

## DE genes: 6185 
## Significant (p_adj < 0.05): 3

Top DE genes: IL-17R+ vs IL-17R- ECM fibroblasts

gene	p_val	avg_log2FC	pct.1	pct.2	p_val_adj
Il17rc	0	1.7631	1.000	0.248	0.0000
Il17ra	0	1.0796	1.000	0.367	0.0000
Ppp1r13b	0	1.4316	0.446	0.143	0.0195

Receptor+ Cells: Key Marker Expression

# Expression of key functional genes in receptor+ vs receptor- cells
key_markers <- c("Il17ra", "Il17rc",
                 "Cxcl12", "Vcam1", "Il6", "Kitl", "Fn1",   # niche support
                 "Mmp3", "Cxcl1", "Ccl2",                    # IL-17 outputs
                 "Acta2", "Fap", "Col1a1",                   # myoCAF
                 "Tgfb1", "Tgfb2")                           # immunosuppression
key_present <- key_markers[key_markers %in% rownames(ecm_fib)]

DotPlot(ecm_fib, features = key_present, group.by = "IL17R_group",
        cols = c("lightgrey", "#762A83"), dot.scale = 7) +
  RotatedAxis() +
  ggtitle("Key Gene Expression: IL-17R+ vs IL-17R- ECM Fibroblasts")

---

11. Export

# Save fibroblast object with all new scores and groupings
qs_save(stroma_fib, file.path(cache_dir, "14_stroma_fib_CAF_IL17_scored.qs"))

# Score summary
write.csv(
  score_summary %>% mutate(across(where(is.numeric) & !n, ~round(., 4))),
  file.path(cache_dir, "14_fibroblast_score_summary.csv"),
  row.names = FALSE
)

# Inter-subtype stats
write.csv(kw_results,
          file.path(cache_dir, "14_KW_intersubtype_stats.csv"),
          row.names = FALSE)

if (exists("pw_results") && nrow(pw_results) > 0) {
  write.csv(pw_results,
            file.path(cache_dir, "14_pairwise_intersubtype_stats.csv"),
            row.names = FALSE)
}

# IL-17 receptor co-expression
if (exists("coexpr_summary")) {
  write.csv(coexpr_summary,
            file.path(cache_dir, "14_IL17R_coexpression.csv"),
            row.names = FALSE)
}

# IL-17-high vs low DE (ECM fibroblasts)
if (exists("de_il17") && nrow(de_il17) > 0) {
  write.csv(de_il17,
            file.path(cache_dir, "14_DE_IL17_high_vs_low_ECM.csv"),
            row.names = FALSE)
}

# IL-17R+ vs R- DE (ECM fibroblasts)
if (exists("de_il17r") && nrow(de_il17r) > 0) {
  write.csv(de_il17r,
            file.path(cache_dir, "14_DE_IL17R_pos_vs_neg_ECM.csv"),
            row.names = FALSE)
}

cat("Outputs saved to:", cache_dir, "\n")

## Outputs saved to: /Users/alasdairduguid/Documents/Projects/BSH start up grant 2024 - scRNAseq miR128a/scRNAseq_proj/analysis_feb26/data

cat("Files written:\n")

## Files written:

cat("  14_stroma_fib_CAF_IL17_scored.qs\n")

##   14_stroma_fib_CAF_IL17_scored.qs

cat("  14_fibroblast_score_summary.csv\n")

##   14_fibroblast_score_summary.csv

cat("  14_KW_intersubtype_stats.csv\n")

##   14_KW_intersubtype_stats.csv

cat("  14_pairwise_intersubtype_stats.csv\n")

##   14_pairwise_intersubtype_stats.csv

cat("  14_IL17R_coexpression.csv\n")

##   14_IL17R_coexpression.csv

cat("  14_DE_IL17_high_vs_low_ECM.csv\n")

##   14_DE_IL17_high_vs_low_ECM.csv

cat("  14_DE_IL17R_pos_vs_neg_ECM.csv\n")

##   14_DE_IL17R_pos_vs_neg_ECM.csv

---

Synthesis

Proposed axis: Vγ6Vδ4 γδ T cells → IL-17A → CNS Meningeal Fibroblasts (ECM) → Mmp3/Mmp13 upregulation + constitutive Cxcl12/Kitl/Vcam1 expression → ECM remodelling, growth factor release, and leukaemia retention in the CNS niche.

Key supporting observations:

• IL-17 receptor (RA/RC heterodimer) is selectively expressed in ECM fibroblasts (~23% co-expression) compared to other CNS fibroblast populations
• IL17R+ ECM fibroblasts already express high levels of pro-leukaemic ligands at baseline (Cxcl12, Kitl, Vcam1, Fn1)
• IL-17A response in ECM fibroblasts is MMP-centric (Mmp3, Mmp13) rather than chemokine-centric, consistent with a matrix remodelling output
• Cxcl12 and Kitl are constitutively expressed across ECM fibroblasts regardless of Mmp3 status; Vcam1 shows genuine enrichment in Mmp3+ cells (~15pp difference). Mmp3 may amplify niche retention by liberating matrix-bound ligand forms rather than by transcriptional co-regulation
• IL17-high ECM fibroblasts are spatially coherent within the UMAP, representing a genuine sub-state
• IL-17 response score correlates positively with niche retention score within ECM fibroblasts (r = 0.163, p = 0.0066)
• IL17R+ cells additionally express Ppp1r13b (iASPP), suggesting pro-survival properties that may maintain this population in the hostile leukaemic niche

Outstanding questions:

1. Is Mmp3/Kitl co-expression functionally linked at the protein level? Spatial transcriptomics or multiplex IF could confirm co-localisation.
2. What determines IL-17 receptor expression within the ECM fibroblast population? A further sub-state may exist below current clustering resolution.
3. Is iASPP (Ppp1r13b) enrichment in IL17R+ fibroblasts a cause or consequence of their niche-supportive identity?
4. NicheNet (Script 15) will provide a complementary, data-driven view of sender-receiver interactions to validate the γδ T cell → ECM fibroblast → leukaemia axis identified here.

Session Information

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] scales_1.4.0       knitr_1.51         RColorBrewer_1.1-3 ggrepel_0.9.6     
##  [5] patchwork_1.3.2    ggplot2_4.0.2      tibble_3.3.1       tidyr_1.3.2       
##  [9] dplyr_1.2.0        UCell_2.10.1       qs2_0.1.7          Seurat_5.4.0      
## [13] SeuratObject_5.3.0 sp_2.2-1          
## 
## loaded via a namespace (and not attached):
##   [1] rstudioapi_0.18.0           jsonlite_2.0.0             
##   [3] magrittr_2.0.4              ggbeeswarm_0.7.3           
##   [5] spatstat.utils_3.2-2        farver_2.1.2               
##   [7] rmarkdown_2.30              zlibbioc_1.52.0            
##   [9] vctrs_0.7.1                 ROCR_1.0-12                
##  [11] spatstat.explore_3.7-0      S4Arrays_1.6.0             
##  [13] htmltools_0.5.9             BiocNeighbors_2.0.1        
##  [15] SparseArray_1.6.2           sass_0.4.10                
##  [17] sctransform_0.4.3           parallelly_1.46.1          
##  [19] KernSmooth_2.23-26          bslib_0.10.0               
##  [21] htmlwidgets_1.6.4           ica_1.0-3                  
##  [23] plyr_1.8.9                  plotly_4.12.0              
##  [25] zoo_1.8-15                  cachem_1.1.0               
##  [27] igraph_2.2.2                mime_0.13                  
##  [29] lifecycle_1.0.5             pkgconfig_2.0.3            
##  [31] Matrix_1.7-4                R6_2.6.1                   
##  [33] fastmap_1.2.0               GenomeInfoDbData_1.2.13    
##  [35] MatrixGenerics_1.18.1       fitdistrplus_1.2-6         
##  [37] future_1.69.0               shiny_1.13.0               
##  [39] digest_0.6.39               S4Vectors_0.44.0           
##  [41] tensor_1.5.1                RSpectra_0.16-2            
##  [43] irlba_2.3.7                 GenomicRanges_1.58.0       
##  [45] labeling_0.4.3              progressr_0.18.0           
##  [47] spatstat.sparse_3.1-0       mgcv_1.9-4                 
##  [49] httr_1.4.8                  polyclip_1.10-7            
##  [51] abind_1.4-8                 compiler_4.4.2             
##  [53] withr_3.0.2                 S7_0.2.1                   
##  [55] BiocParallel_1.40.2         fastDummies_1.7.5          
##  [57] MASS_7.3-65                 DelayedArray_0.32.0        
##  [59] tools_4.4.2                 vipor_0.4.7                
##  [61] lmtest_0.9-40               otel_0.2.0                 
##  [63] beeswarm_0.4.0              httpuv_1.6.16              
##  [65] future.apply_1.20.2         goftest_1.2-3              
##  [67] glue_1.8.0                  nlme_3.1-168               
##  [69] promises_1.5.0              grid_4.4.2                 
##  [71] Rtsne_0.17                  cluster_2.1.8.2            
##  [73] reshape2_1.4.5              generics_0.1.4             
##  [75] gtable_0.3.6                spatstat.data_3.1-9        
##  [77] data.table_1.18.2.1         XVector_0.46.0             
##  [79] stringfish_0.18.0           BiocGenerics_0.52.0        
##  [81] spatstat.geom_3.7-0         RcppAnnoy_0.0.23           
##  [83] RANN_2.6.2                  pillar_1.11.1              
##  [85] stringr_1.6.0               limma_3.62.2               
##  [87] spam_2.11-3                 RcppHNSW_0.6.0             
##  [89] later_1.4.8                 splines_4.4.2              
##  [91] lattice_0.22-9              survival_3.8-6             
##  [93] deldir_2.0-4                tidyselect_1.2.1           
##  [95] SingleCellExperiment_1.28.1 miniUI_0.1.2               
##  [97] pbapply_1.7-4               gridExtra_2.3              
##  [99] IRanges_2.40.1              SummarizedExperiment_1.36.0
## [101] scattermore_1.2             stats4_4.4.2               
## [103] xfun_0.56                   Biobase_2.66.0             
## [105] statmod_1.5.1               matrixStats_1.5.0          
## [107] UCSC.utils_1.2.0            stringi_1.8.7              
## [109] lazyeval_0.2.2              yaml_2.3.12                
## [111] evaluate_1.0.5              codetools_0.2-20           
## [113] cli_3.6.5                   uwot_0.2.4                 
## [115] RcppParallel_5.1.11-2       xtable_1.8-8               
## [117] reticulate_1.45.0           jquerylib_0.1.4            
## [119] GenomeInfoDb_1.42.3         dichromat_2.0-0.1          
## [121] Rcpp_1.1.1                  globals_0.19.0             
## [123] spatstat.random_3.4-4       png_0.1-8                  
## [125] ggrastr_1.0.2               spatstat.univar_3.1-6      
## [127] parallel_4.4.2              presto_1.0.0               
## [129] dotCall64_1.2               listenv_0.10.1             
## [131] viridisLite_0.4.3           ggridges_0.5.7             
## [133] crayon_1.5.3                purrr_1.2.1                
## [135] rlang_1.1.7                 cowplot_1.2.0