Niche Composition: Pseudobulk Treatment Comparison

Per-brain aggregation of plaque-proximal cellular niches, Adu vs IgG

Setup

source("R/setup.R")
library(patchwork)

data <- load_slides_and_distances("full")
slide1 <- data$slide1; slide2 <- data$slide2
cell_distances_1 <- data$distances_1; cell_distances_2 <- data$distances_2
cell_distances_all <- data$distances_all

sample_lookup <- cell_distances_all |> dplyr::select(cell, sample_ID, Treatment.Group)

Motivation

The treatment comparison in notebook 11 (11_neighborhood_enrichment.qmd) pools all cells from Adu-treated brains and compares them against all cells from IgG-treated brains, treating each cell as an independent observation. This inflates effective sample size and ignores the nested structure of the experiment: 3 biological replicates per group.

Here we take a pseudobulk approach: niche proportions are first averaged within each brain, then compared across the 3 Adu brains vs 3 IgG brains. This makes the biological replication explicit and is the appropriate unit of inference.

Note on statistical power. With n = 3 replicates per group, all tests are severely underpowered. Welch’s t-test is used here because, unlike the Wilcoxon rank-sum test (minimum two-sided p-value of 0.1 at n = 3 vs n = 3), it can in principle reach p < 0.05 when between-group differences are large relative to within-group variance. Normality cannot be verified at n = 3, so results should be interpreted as exploratory and directional.

Compute per-brain niche composition

For every plaque-proximal cell (< 50 µm from the nearest plaque), the 15 nearest spatial neighbours are identified within the same tissue slide, and the proportion of each annotated cell type among those neighbours is recorded. Proportions are then averaged to the brain level, yielding 6 brain-level observations per (plaque type × neighbour cell type) combination.

compute_niche_composition <- function(seurat_obj, cell_distances, threshold = 50) {
  proximal_info <- cell_distances |>
    mutate(
      nearest_type = if_else(dist_caa < dist_paren, "CAA", "Parenchymal"),
      min_dist     = pmin(dist_caa, dist_paren)
    ) |>
    filter(min_dist < threshold)

  if (nrow(proximal_info) == 0) return(tibble())

  coords <- GetTissueCoordinates(seurat_obj) |> as.data.frame()
  meta   <- seurat_obj@meta.data |> as_tibble(rownames = "cell")

  all_mat  <- coords |> dplyr::select(x, y) |> as.matrix()
  prox_idx <- match(proximal_info$cell, coords$cell)
  prox_idx <- prox_idx[!is.na(prox_idx)]

  nn_res <- nn2(data = all_mat, query = all_mat[prox_idx, ], k = 16)

  map(seq_along(prox_idx), \(i) {
    neighbour_idx <- nn_res$nn.idx[i, -1]
    tibble(
      focal_cell     = coords$cell[prox_idx[i]],
      neighbour_type = meta$annotatedclusters[neighbour_idx]
    )
  }) |>
    bind_rows() |>
    count(focal_cell, neighbour_type) |>
    complete(focal_cell, neighbour_type, fill = list(n = 0L)) |>
    group_by(focal_cell) |>
    mutate(prop = n / sum(n)) |>
    ungroup() |>
    left_join(
      proximal_info |> dplyr::select(cell, nearest_type),
      by = c("focal_cell" = "cell")
    ) |>
    filter(!is.na(nearest_type))
}

niche_1 <- compute_niche_composition(slide1, cell_distances_1)
niche_2 <- compute_niche_composition(slide2, cell_distances_2)

niche_all <- bind_rows(
  niche_1 |> mutate(slide = 1L),
  niche_2 |> mutate(slide = 2L)
) |>
  left_join(sample_lookup, by = c("focal_cell" = "cell"))

Pseudobulk aggregation

niche_brain <- niche_all |>
  group_by(sample_ID, Treatment.Group, nearest_type, neighbour_type) |>
  summarise(mean_prop = mean(prop, na.rm = TRUE), .groups = "drop")

Treatment comparison

Dot plots: per-brain niche proportions

Each point represents one brain (n = 3 per group). Crossbars show group means. Plotting individual brains makes consistency (or inconsistency) of treatment effects across replicates directly visible.

niche_brain |>
  ggplot(aes(x = neighbour_type, y = mean_prop, colour = Treatment.Group)) +
  geom_point(
    position = position_dodge(width = 0.6),
    size = 2.5, alpha = 0.85
  ) +
  stat_summary(
    aes(group = Treatment.Group),
    fun = mean, geom = "crossbar",
    width = 0.45,
    position = position_dodge(width = 0.6),
    linewidth = 0.5
  ) +
  scale_colour_manual(values = c(Adu = "#E05A4E", IgG = "#4E79A7")) +
  facet_wrap(~ nearest_type, ncol = 1) +
  labs(
    title = "Niche composition: Adu vs IgG (pseudobulk, n = 3 per group)",
    x = "Neighbour cell type",
    y = "Mean niche proportion",
    colour = "Treatment"
  ) +
  theme_minimal(base_size = 13) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Welch’s t-test (exploratory)

ttest_results <- niche_brain |>
  group_by(nearest_type, neighbour_type) |>
  summarise(
    mean_Adu = mean(mean_prop[Treatment.Group == "Adu"], na.rm = TRUE),
    mean_IgG = mean(mean_prop[Treatment.Group == "IgG"], na.rm = TRUE),
    delta    = mean_Adu - mean_IgG,
    log2fc   = log2((mean_Adu + 1e-6) / (mean_IgG + 1e-6)),
    p_value  = tryCatch(
      t.test(
        mean_prop[Treatment.Group == "Adu"],
        mean_prop[Treatment.Group == "IgG"]
      )$p.value,
      error = \(e) NA_real_
    ),
    .groups = "drop"
  ) |>
  mutate(p_adj = p.adjust(p_value, method = "BH")) |>
  arrange(nearest_type, p_value)

ttest_results |>
  mutate(across(c(mean_Adu, mean_IgG, delta, log2fc), \(x) round(x, 4))) |>
  mutate(across(c(p_value, p_adj), \(x) round(x, 3))) |>
  kableExtra::kbl(
    caption = "Pseudobulk Welch's t-test: Adu vs IgG niche proportions (n = 3 per group)",
    col.names = c(
      "Plaque type", "Neighbour type",
      "Mean Adu", "Mean IgG", "Δ (Adu − IgG)", "log₂FC", "p-value", "p adj (BH)"
    )
  ) |>
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = FALSE
  )

Pseudobulk Welch's t-test: Adu vs IgG niche proportions (n = 3 per group)

Plaque type	Neighbour type	Mean Adu	Mean IgG	Δ (Adu − IgG)	log₂FC	p-value	p adj (BH)
CAA	GABAergic Neuron 3	0.0193	0.0256	-0.0063	-0.4095	0.036	0.504
CAA	Microglia	0.0825	0.0598	0.0227	0.4650	0.042	0.504
CAA	T Cell	0.0024	0.0008	0.0016	1.5659	0.056	0.504
CAA	BAM	0.0342	0.0219	0.0123	0.6437	0.078	0.504
CAA	Glutamatergic Neuron 2	0.0971	0.1665	-0.0694	-0.7778	0.088	0.504
CAA	GABAergic Neuron 1	0.0477	0.0196	0.0281	1.2821	0.150	0.567
CAA	OPC	0.0210	0.0186	0.0024	0.1761	0.195	0.641
CAA	Ependymal	0.0010	0.0090	-0.0080	-3.2282	0.329	0.765
CAA	Oligodendrocyte	0.1020	0.0845	0.0175	0.2708	0.333	0.765
CAA	Astrocyte	0.1824	0.1703	0.0121	0.0994	0.369	0.777
CAA	Glutamatergic Neuron 3	0.0014	0.0086	-0.0072	-2.6434	0.423	0.777
CAA	VSMC	0.0693	0.0618	0.0074	0.1639	0.518	0.851
CAA	GABAergic Neuron 4	0.0067	0.0078	-0.0011	-0.2109	0.589	0.934
CAA	NA	0.0260	0.0275	-0.0015	-0.0812	0.696	0.941
CAA	Endothelial	0.0932	0.0970	-0.0038	-0.0574	0.722	0.941
CAA	VLMC	0.0117	0.0169	-0.0053	-0.5385	0.733	0.941
CAA	Glutamatergic Neuron 5	0.0005	0.0006	-0.0001	-0.2186	0.803	0.941
CAA	Pericyte	0.0115	0.0112	0.0003	0.0442	0.816	0.941
CAA	Glutamatergic Neuron 1	0.0382	0.0400	-0.0019	-0.0689	0.876	0.941
CAA	GABAergic Neuron 2	0.0145	0.0152	-0.0006	-0.0631	0.876	0.941
CAA	Glutamatergic Neuron 4	0.0044	0.0042	0.0002	0.0784	0.879	0.941
CAA	CP	0.0013	0.0015	-0.0002	-0.2276	0.902	0.943
CAA	Fibroblast	0.1318	0.1311	0.0007	0.0074	0.940	0.961
Parenchymal	Microglia	0.1543	0.1223	0.0320	0.3356	0.014	0.504
Parenchymal	T Cell	0.0032	0.0016	0.0016	1.0243	0.048	0.504
Parenchymal	OPC	0.0306	0.0276	0.0030	0.1476	0.071	0.504
Parenchymal	NA	0.0266	0.0343	-0.0077	-0.3656	0.127	0.567
Parenchymal	Fibroblast	0.0099	0.0134	-0.0035	-0.4414	0.144	0.567
Parenchymal	VSMC	0.0057	0.0071	-0.0014	-0.3102	0.157	0.567
Parenchymal	Pericyte	0.0187	0.0166	0.0021	0.1705	0.160	0.567
Parenchymal	BAM	0.0037	0.0031	0.0006	0.2432	0.217	0.665
Parenchymal	Oligodendrocyte	0.2117	0.2067	0.0050	0.0343	0.270	0.730
Parenchymal	GABAergic Neuron 3	0.0101	0.0112	-0.0011	-0.1498	0.280	0.730
Parenchymal	CP	0.0001	0.0010	-0.0009	-3.3708	0.286	0.730
Parenchymal	Glutamatergic Neuron 1	0.1052	0.1169	-0.0117	-0.1518	0.383	0.777
Parenchymal	Ependymal	0.0025	0.0041	-0.0016	-0.7175	0.394	0.777
Parenchymal	GABAergic Neuron 4	0.0097	0.0114	-0.0017	-0.2291	0.415	0.777
Parenchymal	Glutamatergic Neuron 3	0.0225	0.0175	0.0050	0.3653	0.451	0.798
Parenchymal	Glutamatergic Neuron 2	0.0676	0.0814	-0.0139	-0.2694	0.482	0.822
Parenchymal	Endothelial	0.0723	0.0708	0.0015	0.0302	0.701	0.941
Parenchymal	VLMC	0.0006	0.0007	-0.0001	-0.3299	0.763	0.941
Parenchymal	GABAergic Neuron 1	0.0461	0.0508	-0.0046	-0.1382	0.794	0.941
Parenchymal	GABAergic Neuron 2	0.0406	0.0428	-0.0022	-0.0770	0.799	0.941
Parenchymal	Astrocyte	0.1340	0.1355	-0.0015	-0.0158	0.827	0.941
Parenchymal	Glutamatergic Neuron 4	0.0206	0.0194	0.0012	0.0874	0.850	0.941
Parenchymal	Glutamatergic Neuron 5	0.0038	0.0039	-0.0001	-0.0433	0.969	0.969

Difference heatmap (Adu − IgG)

Red tiles indicate cell types proportionally more abundant in Adu niche; blue tiles indicate depletion relative to IgG. The colour scale shows the raw difference in mean niche proportion.

ttest_results |>
  ggplot(aes(x = neighbour_type, y = nearest_type, fill = delta)) +
  geom_tile(colour = "white") +
  scale_fill_gradient2(
    low = "#4E79A7", mid = "white", high = "#E05A4E", midpoint = 0,
    name = "Δ proportion\n(Adu − IgG)"
  ) +
  labs(
    title = "Treatment difference in niche composition (pseudobulk)",
    subtitle = "Δ = mean(Adu) − mean(IgG) per brain; n = 3 per group",
    x = "Neighbour cell type",
    y = "Nearest plaque type"
  ) +
  theme_minimal(base_size = 15) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))