Summary

Parameters

1. Data

2. Re-annotating plasma cells

3. Epithelium, April and immune cells

4. PC density by stage

5. Do PCs co-localise with April+ epithelium?

6. Sensitivity checks

7. Do PCs sit in April-rich neighbourhoods?

8. What do PCs sit next to?

9. Which cells express April?

10. Can PCs respond to April?

11. Spatial maps

12. Overall conclusions

13. Caveats and possible next steps?

Session info

Collaboration_JulieTellier_jChainPcAprilEpi_coLocalisation_xenium_20261001

2a. Gate plot

2b. Does the support-marker rule throw out real PCs?

2c. Are the new PCs really PCs, or endothelial cells / adipocytes?

2d. The supplied labels

Ashley L Weir

01 October 2026

Question. Do Jchain-expressing plasma cells (PCs) co-localise with April (Tnfsf13) expressing epithelial cells in the mouse mammary gland, across four stages: virgin, pregnancy, lactation and post-weaning.

Main findings

The plasma cells had to be re-annotated. The supplied predicted_cell_type labels contain no plasma cell type. Most real PCs had been labelled “Endothelial cell”. I called PCs from markers instead, and checked that they are dominated by immunoglobulin transcripts (IgA PCs).
PCs are rare except after weaning. About 1 PC per mm² in pregnancy and lactation, about 19 per mm² in post-weaning, and none in virgin.
PCs do not preferentially sit next to April+ epithelium, in any stage, and do not sit in April-rich neighbourhoods. This held across radii, two different random baselines and alternative cell definitions.
PCs do have a spatial preference. In lactation they sit in the stroma between alveoli. After weaning they gather around the remaining epithelial structures (about 1.5 to 3× enriched near epithelium and myoepithelium) and cluster with each other.
PCs express the April receptor BCMA (Tnfrsf17), about 25 to 45% detected against under 3% of other cells, so they are able to respond to April.

Main caveats. There is one tissue section (one mouse) per stage. Tnfsf13 is detected at low levels everywhere. Only post-weaning has enough PCs for well-powered spatial tests.

Every threshold used in the analysis is set here, so it can be changed in one place.

H5AD        <- "MG_combined_FORCT.h5ad"
OUT_DIR     <- "coloc_out"

PIXEL_SIZE  <- 0.2125  # microns per pixel (Xenium); x/y_centroid are in pixels
JCHAIN_MIN  <- 5       # min Jchain counts for a PC
SUPPORT_MIN <- 1       # min number of other PC markers detected
PC_EPI_MAX  <- 2       # max epithelial counts (Krt18+Epcam+Krt19+Cdh1) in a PC
KRT18_MIN   <- 2       # min Krt18 counts for epithelium
APRIL_MIN   <- 1       # min Tnfsf13 counts for April+
PTPRC_MIN   <- 1       # min Ptprc counts for the immune comparison set
MIN_PC      <- 10      # skip a section with fewer PCs than this
RADII       <- c(10, 15, 20, 30, 50, 75, 100)   # microns
N_PERM      <- 1000    # random shuffles per test
N_RANDOM    <- 5000    # random reference cells per section (April neighbourhood)
SEED        <- 42

STAGES <- c(MG3 = "Virgin", MG5 = "Pregnancy",
            MG4 = "Lactation", MG1 = "Post-weaning")

dir.create(OUT_DIR, showWarnings = FALSE)

The file MG_combined_FORCT.h5ad was provided by Julie Tellier. It contains 480,372 cells and 5,106 genes from four Xenium sections, one per stage. It had already been QC-filtered, log-normalised and integrated with scVI ( I think). Raw counts are stored in the counts layer, and were used throughout.

Why only a few genes are extracted: the full matrix is large. Every analysis here needs only about 30 marker genes, so I keep their raw counts in a small table (df) with one row per cell.

Why coordinates are converted: in this file x_centroid / y_centroid are in pixels, not microns. I found this when no cell appeared to have any neighbour within 15 “µm”, which is impossible. Multiplying by 0.2125 µm/pixel gives a typical spacing between neighbouring cells of about 8 µm, as expected.

sce <- readH5AD(H5AD, reader = "R", use_hdf5 = TRUE)
assayNames(sce)[assayNames(sce) == "X"] <- "logcounts"

markers <- list(
  plasma = c("Jchain", "Mzb1", "Xbp1", "Prdm1", "Sdc1", "Tnfrsf17", "Irf4",
             "Igkc", "Igha", "Ighm"),
  bcell  = c("Ms4a1", "Cd19", "Pax5", "Cd79a", "Cd79b"),
  epi    = c("Krt18", "Epcam", "Krt19", "Cdh1"),
  immune = c("Ptprc"),
  endo   = c("Pecam1", "Cdh5", "Kdr"),
  adipo  = c("Fabp4", "Adipoq", "Plin1"),
  april  = c("Tnfsf13", "Tnfrsf17", "Tnfrsf13b", "Tnfrsf13c")
)
genes <- intersect(unlist(markers), rownames(sce))
cnt   <- t(as.matrix(counts(sce)[genes, ]))

df <- data.frame(
  cell         = colnames(sce),
  sample       = sub("_.*", "", as.character(sce$sample_id)),
  x            = sce$x_centroid * PIXEL_SIZE,
  y            = sce$y_centroid * PIXEL_SIZE,
  total_counts = sce$total_counts,
  predicted    = as.character(sce$predicted_cell_type),
  cnt,
  check.names  = FALSE
)
df$stage <- factor(STAGES[df$sample], levels = STAGES)
stopifnot(!anyNA(df$stage))
missing_markers <- setdiff(unlist(markers), rownames(sce))

Markers not on the panel (skipped automatically): .

nnd <- RANN::nn2(df[df$stage == "Lactation", c("x", "y")], k = 2)$nn.dists[, 2]
g   <- function(gene) if (gene %in% names(df)) df[[gene]] else rep(0, nrow(df))

Units check: the median distance from a lactation cell to its nearest neighbour is 8.2 µm. It should be about 5 to 15 µm.

Why: the supplied predicted_cell_type labels have no plasma cell type, so PCs had been assigned to other types (see 2d). Jchain alone is not an ideal marker… PCs make so much Ig and Jchain mRNA that some can spill into neighbouring cells in Xenium, and epithelial cells next to a PC can look Jchain+. Calling those cells PCs would partly build the co-localisation we want to test into the labels themselves.

Rule: a cell is a PC if it has

Jchain ≥ 5 counts, and
at least 1 other PC marker detected (Mzb1, Xbp1, Prdm1, Sdc1, Tnfrsf17, Irf4), and
almost no epithelial signal (Krt18 + Epcam + Krt19 + Cdh1 ≤ 2).

pc_support_genes <- intersect(c("Mzb1", "Xbp1", "Prdm1", "Sdc1",
                                "Tnfrsf17", "Irf4"), names(df))
df$epi_signal <- g("Krt18") + g("Epcam") + g("Krt19") + g("Cdh1")
df$pc_support <- rowSums(sapply(pc_support_genes, function(x) df[[x]] > 0))

df$is_PC <- df$Jchain >= JCHAIN_MIN &
            df$pc_support >= SUPPORT_MIN &
            df$epi_signal <= PC_EPI_MAX

Each cell is plotted by its Jchain and epithelial counts. “Real” PCs sit top-left (high Jchain, no keratin) and epithelium sits bottom-right. Cells with both (on the diagonal) would indicate spillover. The red lines show the thresholds.

ggplot(df, aes(log1p(epi_signal), log1p(Jchain))) +
  geom_bin2d(bins = 80) +
  scale_fill_viridis_c(trans = "log10") +
  geom_hline(yintercept = log1p(JCHAIN_MIN), linetype = 2, colour = "red") +
  geom_vline(xintercept = log1p(PC_EPI_MAX), linetype = 2, colour = "red") +
  facet_wrap(~stage) +
  labs(x = "log1p(Krt18 + Epcam + Krt19 + Cdh1 counts)",
       y = "log1p(Jchain counts)")

Conclusion: the gate separates PCs cleanly. Epithelial cells mostly have 0 to 3 Jchain counts, and the top-right corner (both signals) is nearly empty, so spillover is limited. Virgin breast has essentially no Jchain-high cells.

kable(as.data.frame.matrix(table(df$stage, df$is_PC)),
      col.names = c("Not PC", "PC"), caption = "PCs per stage")

PCs per stage
	Not PC	PC
Virgin	42451	0
Pregnancy	132554	18
Lactation	234273	47
Post-weaning	70749	280

hiJ <- df$Jchain >= JCHAIN_MIN & df$epi_signal <= PC_EPI_MAX
kable(as.data.frame.matrix(table(df$stage[hiJ], df$pc_support[hiJ])),
      caption = "High-Jchain, low-keratin cells by number of support markers detected")

zero <- hiJ & df$pc_support == 0
kable(head(sort(table(df$predicted[zero]), decreasing = TRUE), 5),
      col.names = c("Predicted label", "n"),
      caption = "Labels of high-Jchain cells with no support markers")

Labels of high-Jchain cells with no support markers
Predicted label	n
Endothelial cell_Fabp4 high	43
Endothelial cell_Aqp1 high	15
Fibroblast_Col3a1 high	3
B cell_Cd79a high	2
Fibroblast	1

Conclusion: most high-Jchain cells have 1 to 3 other PC markers. The few with none (mostly in post-weaning) are labelled endothelial or fibroblast. They are most likely picking up Jchain from nearby PCs, so I have excluded them to be cautious.

Why: most cells I have now labelled as PCs had previously been labelled “Endothelial cell”. I compared PCs with the endothelial-labelled cells on endothelial genes (Pecam1, Cdh5, Kdr), adipocyte genes (Fabp4, Adipoq, Plin1) and immunoglobulin genes.

grp <- ifelse(df$is_PC, "PC",
       ifelse(grepl("Endothelial", df$predicted), "Endo-labelled", "Other"))

check_genes <- intersect(c("Pecam1", "Cdh5", "Kdr", "Fabp4", "Adipoq", "Plin1"), names(df))
kable(round(100 * sapply(check_genes, function(x) tapply(df[[x]] > 0, grp, mean)), 1),
      caption = "% of cells detecting each gene")

% of cells detecting each gene
	Pecam1	Cdh5	Kdr	Fabp4	Adipoq	Plin1
Endo-labelled	16.0	18.8	11.5	93.9	45.4	43.1
Other	2.6	3.6	2.5	40.3	5.0	4.5
PC	6.1	6.1	2.3	87.8	33.0	29.6

level_genes <- intersect(c("Jchain", "Igkc", "Igha", "Fabp4", "Adipoq", "total_counts"), names(df))
kable(sapply(level_genes, function(x) tapply(df[[x]], grp, median)),
      caption = "Median counts per cell")

Median counts per cell
	Jchain	Igkc	Igha	Fabp4	total_counts
Endo-labelled	0	0	0	13	150
Other	0	0	0	0	289
PC	10	16	17	8	198

Conclusion: the new PCs seem to be genuine IgA plasma cells.

Endothelial genes are low in PCs, about a third of the level in endothelial-labelled cells.
The median PC has around 40 Jchain + Igkc + Igha counts, roughly a fifth of its transcripts. Endothelial-labelled cells have none.
Igha is as high as Igkc, consistent with IgA PCs, the main PC type in the mammary gland.
Adipocyte genes are detected in PCs, but at lower levels than Ig genes. I think this is background from the surrounding fat pad.

kable(head(sort(table(df$predicted[df$is_PC]), decreasing = TRUE), 6),
      col.names = c("Predicted label", "n PCs"),
      caption = "Previous labels of the re-annotated PCs")

Previous labels of the re-annotated PCs
Predicted label	n PCs
Endothelial cell_Fabp4 high	247
Endothelial cell_Aqp1 high	47
Secretory alveoli cell_Csn2 high	14
Cd8 positive T cell_Cd8b1 high	11
B cell_Cd79a high	6
NK cell_Gzma high	5

Conclusion: about 85% of PCs had been labelled “Endothelial cell”. Separately, about 36% of all cells carry an endothelial label, which seems like far too many for mammary gland. Since Fabp4 is also a strong adipocyte marker, “Endothelial cell_Fabp4 high” is probably largely adipocytes. For immune and stromal cells, i don’t think predicted_cell_type should not be relied on. All main analyses below use my marker-based calls only.

Definitions (all from raw counts):

Epithelial: Krt18 ≥ 2, and not a PC (based on Julie’s knowledge)
April+: Tnfsf13 ≥ 1
Immune: Ptprc (CD45) ≥ 1 and not epithelial, plus all PCs, since I think PCs often have low CD45. This group is only used as a comparison set.

df$is_epi    <- g("Krt18") >= KRT18_MIN & !df$is_PC
df$april     <- g("Tnfsf13") >= APRIL_MIN
df$is_immune <- (g("Ptprc") >= PTPRC_MIN & !df$is_epi) | df$is_PC

epi_summary <- df |>
  group_by(stage) |>
  summarise(cells = n(), PCs = sum(is_PC), epithelial = sum(is_epi),
            april_epithelial = sum(is_epi & april),
            pct_epi_april = round(100 * mean(april[is_epi]), 1),
            immune = sum(is_immune), .groups = "drop")
kable(epi_summary, caption = "Cell counts per stage")

Cell counts per stage
stage	cells	PCs	epithelial	april_epithelial	pct_epi_april	immune
Virgin	42451	0	475	15	3.2	614
Pregnancy	132572	18	31203	6261	20.1	2967
Lactation	234320	47	33542	5533	16.5	1844
Post-weaning	71029	280	2438	45	1.8	1019

Conclusion: pregnancy and lactation have plenty of epithelium, 16 to 20% of it April+. Post-weaning has very little epithelium (involution), with only about 45 April+ epithelial cells, yet it is the stage with the most PCs. Loosening the epithelial definition only raised post-weaning to about 3,700 to 4,100 epithelial cells, so the epithelium really is sparse there.

Why: sections differ in size and composition, so raw PC counts are not comparable. Tissue area is estimated by counting the 50 × 50 µm squares that contain at least one cell.

pc_density <- df |>
  group_by(stage) |>
  summarise(PCs = sum(is_PC), cells = n(), epithelial = sum(is_epi),
            area_mm2 = n_distinct(paste(floor(x / 50), floor(y / 50))) * 2500 / 1e6,
            .groups = "drop") |>
  mutate(PC_per_mm2       = PCs / area_mm2,
         PC_per_10k_cells = 1e4 * PCs / cells,
         PC_per_1k_epi    = 1e3 * PCs / epithelial)
kable(pc_density, digits = 2, caption = "PC density")

PC density
stage	PCs	cells	epithelial	area_mm2	PC_per_mm2	PC_per_10k_cells	PC_per_1k_epi
Virgin	0	42451	475	15.36	0.00	0.00	0.00
Pregnancy	18	132572	31203	19.92	0.90	1.36	0.58
Lactation	47	234320	33542	41.35	1.14	2.01	1.40
Post-weaning	280	71029	2438	15.02	18.65	39.42	114.85

Conclusion: PCs are about 16× denser per mm² after weaning than in lactation, and absent in virgin. Per mm² is the fairest measure. I suspect per epithelial cell is inflated in post-weaning because the epithelium shrinks. This is notable as I thought mammary IgA PCs are usually thought to build up during pregnancy and lactation. Here they peak after weaning, possibly because PCs persist while the epithelium regresses?

Consequence for the analysis: only post-weaning has enough PCs for well-powered spatial tests. Lactation (47) and pregnancy (18) are underpowered, and virgin (0) is skipped.

Approach. For each section and each radius (10 to 100 µm), I take every epithelial cell within that radius of any PC and calculate the fraction that is April+. That observed value is compared against two random baselines, each built from 1000 shuffles:

Baseline 1, shuffle April among epithelium: asks whether epithelium next to PCs is more often April+ than epithelium in general.
Baseline 2, shuffle PC identity among immune cells: asks whether PCs sit next to April+ epithelium more than other immune cells do. This controls for immune cells in general favouring certain regions.

Enrichment = observed / mean of shuffles (1 = no effect). The p-value is one-sided (observed ≥ shuffled).

## epithelial neighbours within `radius` of each query cell
epi_neighbours <- function(query_xy, epi_xy, radius) {
  k  <- min(nrow(epi_xy), 1000)
  nn <- RANN::nn2(epi_xy, query_xy, k = k, searchtype = "radius", radius = radius)
  if (k < nrow(epi_xy) && any(nn$nn.idx[, k] > 0))
    warning("radius search hit k = ", k, " at radius ", radius)
  lapply(seq_len(nrow(nn$nn.idx)), function(i) { v <- nn$nn.idx[i, ]; v[v > 0] })
}

## one stage, one radius
coloc_one <- function(d, radius, n_perm = N_PERM, seed = SEED) {
  set.seed(seed)
  epi   <- d[d$is_epi, ]
  imm   <- d[d$is_immune, ]
  april <- epi$april
  nb    <- epi_neighbours(as.matrix(imm[, c("x", "y")]),
                          as.matrix(epi[, c("x", "y")]), radius)
  n_e   <- lengths(nb)
  n_a   <- vapply(nb, function(v) sum(april[v]), numeric(1))
  isPC  <- imm$is_PC
  n_pc  <- sum(isPC)

  obs   <- sum(n_a[isPC]) / sum(n_e[isPC])
  null1 <- replicate(n_perm, mean(sample(april)[unlist(nb[isPC])]))
  null2 <- replicate(n_perm, { s <- sample.int(length(isPC), n_pc)
                               sum(n_a[s]) / sum(n_e[s]) })
  data.frame(
    radius           = radius,
    n_PC             = n_pc,
    n_epi_neighbours = sum(n_e[isPC]),
    obs_frac_april   = obs,
    enrich_vs_epi    = obs / mean(null1),
    p_vs_epi         = if (is.na(obs)) NA else (sum(null1 >= obs) + 1) / (n_perm + 1),
    enrich_vs_immune = obs / mean(null2, na.rm = TRUE),
    p_vs_immune      = if (is.na(obs)) NA else (sum(null2 >= obs, na.rm = TRUE) + 1) / (n_perm + 1)
  )
}

run_coloc <- function(df) {
  df |>
    group_split(stage) |>
    lapply(function(d) {
      if (sum(d$is_PC) < MIN_PC) return(NULL)
      cbind(stage = d$stage[1], do.call(rbind, lapply(RADII, coloc_one, d = d)))
    }) |>
    bind_rows()
}

plot_enrichment <- function(res) {
  res |>
    select(stage, radius, enrich_vs_epi, enrich_vs_immune) |>
    pivot_longer(starts_with("enrich"), names_to = "baseline", values_to = "enrichment") |>
    mutate(baseline = recode(baseline,
                             enrich_vs_epi    = "vs all epithelium (baseline 1)",
                             enrich_vs_immune = "vs other immune cells (baseline 2)")) |>
    ggplot(aes(radius, enrichment, colour = stage)) +
    geom_hline(yintercept = 1, linetype = 2) +
    geom_line() + geom_point() +
    facet_wrap(~baseline) +
    labs(x = "Radius around each PC (µm)",
         y = "April+ fraction of neighbouring epithelium\n(observed / shuffled)",
         colour = NULL)
}

res <- run_coloc(df)
write.csv(res, file.path(OUT_DIR, "colocalisation_results.csv"), row.names = FALSE)

kable(res, digits = 3, caption = "PC / April+ epithelium co-localisation")

PC / April+ epithelium co-localisation
stage	radius	n_PC	n_epi_neighbours	obs_frac_april	enrich_vs_epi	p_vs_epi	enrich_vs_immune	p_vs_immune
Pregnancy	10	18	1	0.000	0.000	1.000	0.000	0.999
Pregnancy	15	18	3	0.333	1.639	0.495	1.646	0.093
Pregnancy	20	18	5	0.200	1.000	0.671	1.001	0.508
Pregnancy	30	18	9	0.111	0.568	0.866	0.575	0.979
Pregnancy	50	18	35	0.171	0.857	0.726	0.881	0.817
Pregnancy	75	18	120	0.175	0.866	0.790	0.896	0.861
Pregnancy	100	18	239	0.155	0.767	0.974	0.788	1.000
Lactation	10	47	7	0.286	1.723	0.328	2.661	0.106
Lactation	15	47	23	0.174	1.056	0.551	1.223	0.317
Lactation	20	47	49	0.184	1.104	0.438	1.255	0.234
Lactation	30	47	121	0.165	0.990	0.554	1.064	0.368
Lactation	50	47	352	0.168	1.012	0.478	1.019	0.448
Lactation	75	47	832	0.172	1.041	0.324	1.014	0.428
Lactation	100	47	1529	0.175	1.061	0.195	1.046	0.269
Post-weaning	10	280	38	0.053	2.878	0.162	3.368	0.056
Post-weaning	15	280	120	0.025	1.388	0.380	0.977	0.498
Post-weaning	20	280	243	0.025	1.356	0.304	0.824	0.675
Post-weaning	30	280	529	0.025	1.333	0.198	1.007	0.487
Post-weaning	50	280	1102	0.026	1.433	0.069	1.047	0.412
Post-weaning	75	280	1997	0.025	1.332	0.048	1.114	0.225
Post-weaning	100	280	2899	0.020	1.069	0.344	1.019	0.439

plot_enrichment(res)

n_epi_neighbours shows how many epithelial cells each value rests on. Rows with very few (under about 30) are unreliable.

Conclusion: no convincing co-localisation in any stage.

Lactation: PCs do have epithelial neighbours (most within 30 µm), but the April+ fraction of that epithelium matches the section background. Enrichment is about 1.0 at every radius.
Post-weaning: a weak 1.3 to 1.4× enrichment against baseline 1 (best p ≈ 0.05, one of many tests), but about 1.0 against baseline 2. PCs are no closer to April+ epithelium than other immune cells are, and the test rests on about 45 April+ cells.
Pregnancy: too few PCs and neighbours to interpret.

Why: to check that the result doesn’t depend on our exact definitions, this test was rerun with (a) epithelium defined from the combined signal (Krt18 + Epcam + Krt19 + Cdh1 ≥ 2) and (b) PCs called without the support-marker requirement.

df_s1 <- df
df_s1$is_epi    <- df_s1$epi_signal >= 2 & !df_s1$is_PC
df_s1$is_immune <- (g("Ptprc") >= PTPRC_MIN & !df_s1$is_epi) | df_s1$is_PC
res_s1 <- run_coloc(df_s1)

df_s2 <- df
df_s2$is_PC     <- df_s2$Jchain >= JCHAIN_MIN & df_s2$epi_signal <= PC_EPI_MAX
df_s2$is_epi    <- g("Krt18") >= KRT18_MIN & !df_s2$is_PC
df_s2$is_immune <- (g("Ptprc") >= PTPRC_MIN & !df_s2$is_epi) | df_s2$is_PC
res_s2 <- run_coloc(df_s2)

write.csv(res_s1, file.path(OUT_DIR, "sensitivity_epi_signal.csv"), row.names = FALSE)
write.csv(res_s2, file.path(OUT_DIR, "sensitivity_no_support.csv"), row.names = FALSE)

plot_enrichment(res_s1) + ggtitle("(a) Epithelium from combined signal")

plot_enrichment(res_s2) + ggtitle("(b) PCs without support-marker requirement")

Conclusion: both checks agree with the main result. With the looser epithelium, enrichment moves even closer to 1, and the weak post-weaning hint disappears. The occasional low p-value rests on about 10 cells and does not recur across runs.

Why the question was reframed: the density results show that post-weaning is the only well-powered stage, but its epithelium is sparse (about 45 April+ cells), so “near April+ epithelium” is a tiny target. April is also secreted, so its source matters less than how much is around. I therefore asked a different question:

For each PC, add up the Tnfsf13 counts of all cells within the radius, using counts rather than a yes/no cutoff to make the most of a sparse signal.
Split the total into epithelial and non-epithelial April, and also compute April per neighbouring cell (corrects for crowded areas).
Compare the PC average with random sets of other immune cells and random cells in the same section (1000 draws each).

april_niche_one <- function(d, radius, n_perm = N_PERM, seed = SEED) {
  set.seed(seed)
  xy    <- as.matrix(d[, c("x", "y")])
  a_all <- d$Tnfsf13
  a_epi <- ifelse(d$is_epi, d$Tnfsf13, 0)
  rand  <- sample(nrow(d), min(N_RANDOM, nrow(d)))
  q     <- unique(c(which(d$is_immune), rand))

  k   <- min(nrow(d), 1000)
  idx <- RANN::nn2(xy, xy[q, , drop = FALSE], k = k,
                   searchtype = "radius", radius = radius)$nn.idx
  hit <- idx > 0
  idx[!hit] <- 1
  sum_all <- rowSums(matrix(a_all[idx], nrow(idx)) * hit) - a_all[q]
  sum_epi <- rowSums(matrix(a_epi[idx], nrow(idx)) * hit) - a_epi[q]
  n_nb    <- rowSums(hit) - 1

  scores <- data.frame(total_april   = sum_all,
                       epi_april     = sum_epi,
                       non_epi_april = sum_all - sum_epi,
                       april_per_nbr = ifelse(n_nb > 0, sum_all / n_nb, 0))
  is_pc  <- d$is_PC[q]
  n_pc   <- sum(is_pc)
  imm_i  <- which(d$is_immune[q])
  rand_i <- which(q %in% rand)

  do.call(rbind, lapply(names(scores), function(sc) {
    v   <- scores[[sc]]
    obs <- mean(v[is_pc])
    n_i <- replicate(n_perm, mean(v[sample(imm_i,  n_pc)]))
    n_r <- replicate(n_perm, mean(v[sample(rand_i, n_pc)]))
    data.frame(radius = radius, score = sc, PC_mean = obs,
               ratio_vs_immune = obs / mean(n_i),
               p_vs_immune     = (sum(n_i >= obs) + 1) / (n_perm + 1),
               ratio_vs_random = obs / mean(n_r),
               p_vs_random     = (sum(n_r >= obs) + 1) / (n_perm + 1))
  }))
}

run_niche <- function(df, radii = c(15, 30, 50, 100)) {
  df |>
    group_split(stage) |>
    lapply(function(d) {
      if (sum(d$is_PC) < MIN_PC) return(NULL)
      cbind(stage = d$stage[1], do.call(rbind, lapply(radii, april_niche_one, d = d)))
    }) |>
    bind_rows()
}

niche <- run_niche(df)
write.csv(niche, file.path(OUT_DIR, "april_neighbourhood.csv"), row.names = FALSE)

kable(niche[niche$stage == "Post-weaning", -1], digits = 3, row.names = FALSE,
      caption = "April around PCs, post-weaning")

April around PCs, post-weaning
radius	score	PC_mean	ratio_vs_immune	p_vs_immune	ratio_vs_random	p_vs_random
15	total_april	0.061	0.800	0.864	1.174	0.308
15	epi_april	0.011	1.204	0.480	1.523	0.331
15	non_epi_april	0.050	0.737	0.926	1.142	0.339
15	april_per_nbr	0.010	0.650	0.917	0.993	0.480
30	total_april	0.254	1.016	0.468	1.432	0.006
30	epi_april	0.046	1.428	0.130	2.516	0.002
30	non_epi_april	0.207	0.947	0.686	1.312	0.039
30	april_per_nbr	0.012	0.821	0.907	1.174	0.152
50	total_april	0.607	1.075	0.174	1.360	0.001
50	epi_april	0.104	1.524	0.006	2.795	0.001
50	non_epi_april	0.504	1.011	0.450	1.235	0.018
50	april_per_nbr	0.012	0.944	0.740	1.194	0.032
100	total_april	1.825	0.956	0.825	1.150	0.009
100	epi_april	0.204	1.395	0.006	1.691	0.001
100	non_epi_april	1.621	0.918	0.949	1.106	0.044
100	april_per_nbr	0.010	0.861	0.996	1.050	0.203

niche |>
  select(stage, radius, score, ratio_vs_immune, ratio_vs_random) |>
  pivot_longer(starts_with("ratio"), names_to = "vs", values_to = "ratio") |>
  mutate(vs = recode(vs, ratio_vs_immune = "vs other immune cells",
                         ratio_vs_random = "vs random cells")) |>
  ggplot(aes(radius, ratio, colour = stage)) +
  geom_hline(yintercept = 1, linetype = 2) +
  geom_line() + geom_point() +
  facet_grid(vs ~ score) +
  labs(x = "Radius around each PC (µm)",
       y = "April around PCs\n(observed / reference)", colour = NULL)

Result: in no stage are PCs in neighbourhoods where cells carry more April (april_per_nbr is at or below 1 everywhere). Lactation is about 1 throughout, and pregnancy mostly below 1. The one signal is in post-weaning: epithelial April around PCs is about 1.5× that around other immune cells.

Is that just PCs being closer to epithelium? I counted epithelial cells within 50 µm of each immune cell in post-weaning.

d   <- df[df$stage == "Post-weaning", ]
imm <- d[d$is_immune, ]
nb  <- epi_neighbours(as.matrix(imm[, c("x", "y")]),
                      as.matrix(d[d$is_epi, c("x", "y")]), 50)
n_epi <- lengths(nb)
set.seed(SEED)
obs  <- mean(n_epi[imm$is_PC])
null <- replicate(N_PERM, mean(n_epi[sample(nrow(imm), sum(imm$is_PC))]))

kable(data.frame(
  group = c("PCs", "Other immune cells"),
  mean_epithelial_neighbours_50um = round(c(obs, mean(n_epi[!imm$is_PC])), 2)))

group	mean_epithelial_neighbours_50um
PCs	3.94
Other immune cells	2.23

PCs vs random immune cells: ratio 1.46, p = 0.001.

Conclusion: PCs have about 1.8× more epithelial neighbours than other immune cells, which accounts for the 1.5× extra epithelial April. PCs are closer to epithelium, but the epithelium next to them is not more April+.

Why: having found that PCs are near epithelium after weaning, I described each PC’s immediate surroundings (within 20 µm, roughly touching) and compared them with the section as a whole. Neighbours are classified using our PC, epithelial (Krt18) and immune (Ptprc) calls first, and then the broad predicted label otherwise. The labels are only a rough guide (please see 2d).

df$nb_class <- sub("_.*", "", df$predicted)
df$nb_class[df$is_immune] <- "Immune (Ptprc+)"
df$nb_class[df$is_epi]    <- "Epithelial (Krt18+)"
df$nb_class[df$is_PC]     <- "Plasma cell"

pc_neighbours <- lapply(c("Lactation", "Post-weaning"), function(s) {
  d    <- df[df$stage == s, ]
  idx  <- RANN::nn2(d[, c("x", "y")], d[d$is_PC, c("x", "y")],
                    k = 50, searchtype = "radius", radius = 20)$nn.idx[, -1]
  near <- d$nb_class[idx[idx > 0]]
  lv   <- sort(unique(na.omit(d$nb_class)))
  data.frame(stage = s, class = lv,
             pct_near_PC = 100 * as.numeric(table(factor(near, lv))) / length(near),
             pct_section = 100 * as.numeric(table(factor(d$nb_class, lv))) / nrow(d))
}) |>
  bind_rows() |>
  mutate(enrichment = pct_near_PC / pct_section) |>
  filter(pct_near_PC > 0) |>
  arrange(stage, desc(pct_near_PC))

kable(pc_neighbours, digits = 2,
      caption = "Neighbours within 20 µm of PCs (enrichment > 1 = more common next to PCs)")

Neighbours within 20 µm of PCs (enrichment > 1 = more common next to PCs)
stage	class	pct_near_PC	pct_section	enrichment
Lactation	Secretory alveoli cell	42.45	50.49	0.84
Lactation	Endothelial cell	26.91	16.28	1.65
Lactation	Epithelial (Krt18+)	10.72	14.31	0.75
Lactation	Fibroblast	6.35	5.29	1.20
Lactation	Myoepithelial cell	5.47	4.35	1.26
Lactation	Epithelial cell	2.63	3.34	0.79
Lactation	Macrophage	2.19	2.23	0.98
Lactation	Immune (Ptprc+)	1.09	0.77	1.43
Lactation	Luminal cell	0.88	1.40	0.62
Lactation	NK cell	0.44	0.40	1.09
Lactation	Plasma cell	0.44	0.02	21.82
Lactation	Cd8 positive T cell	0.22	0.09	2.33
Lactation	Proliferating T cell	0.22	0.24	0.90
Post-weaning	Endothelial cell	69.78	82.95	0.84
Post-weaning	Epithelial (Krt18+)	9.39	3.43	2.74
Post-weaning	Fibroblast	9.16	8.24	1.11
Post-weaning	Myoepithelial cell	4.21	1.23	3.43
Post-weaning	Plasma cell	3.09	0.39	7.84
Post-weaning	Macrophage	2.16	2.32	0.93
Post-weaning	Immune (Ptprc+)	1.35	1.04	1.30
Post-weaning	Secretory alveoli cell	0.81	0.38	2.16
Post-weaning	Luminal cell	0.04	0.02	2.50

Conclusion: PCs may move from the stroma to the epithelium after weaning.

Lactation: PCs are in the stroma between alveoli. They are enriched next to endothelial/adipocyte cells (about 1.7×), myoepithelium and fibroblasts (about 1.2×), and slightly depleted next to alveolar and Krt18+ epithelium (about 0.75 to 0.85×).
Post-weaning: PCs gather around the remaining ducts and alveoli. They are enriched next to myoepithelium (about 3.4×), Krt18+ epithelium (about 2.7×) and leftover alveolar cells (about 2×), and depleted next to endothelial/adipocyte cells.
In both stages PCs primarily cluster with each other (about 8 to 20×).

df$broad <- sub("_.*", "", df$predicted)
df$broad[df$is_PC] <- "Plasma cell (ours)"
april_by_type <- df |>
  group_by(stage, broad) |>
  summarise(n = n(), pct_april = 100 * mean(Tnfsf13 > 0),
            mean_april = mean(Tnfsf13), .groups = "drop") |>
  filter(n >= 50) |>
  group_by(stage) |>
  slice_max(pct_april, n = 5) |>
  ungroup()
kable(april_by_type, digits = 2, caption = "Top 5 Tnfsf13-expressing cell types per stage")

Top 5 Tnfsf13-expressing cell types per stage
stage	broad	n	pct_april	mean_april
Virgin	Macrophage	1448	10.43	0.12
Virgin	Myoepithelial cell	315	3.81	0.05
Virgin	Mast cell	54	3.70	0.04
Virgin	Fibroblast	4029	3.70	0.04
Virgin	Endothelial cell	35841	1.10	0.01
Pregnancy	Epithelial cell	54344	18.43	0.21
Pregnancy	Secretory alveoli cell	10470	17.25	0.20
Pregnancy	Macrophage	4519	10.78	0.12
Pregnancy	Luminal cell	1102	6.44	0.07
Pregnancy	Myoepithelial cell	11598	5.98	0.06
Lactation	Secretory alveoli cell	149577	11.77	0.13
Lactation	Epithelial cell	9262	9.61	0.10
Lactation	Macrophage	6157	8.01	0.09
Lactation	Cd8 positive T cell	339	5.60	0.06
Lactation	Luminal cell	3669	3.82	0.04
Post-weaning	Macrophage	1971	7.05	0.07
Post-weaning	Luminal cell	120	2.50	0.03
Post-weaning	Secretory alveoli cell	1186	2.28	0.02
Post-weaning	Fibroblast	5918	2.11	0.02
Post-weaning	Myoepithelial cell	1778	1.86	0.02

Conclusion: Tnfsf13 is low everywhere. At most about 18% of any cell type has a single transcript, so any April readout from this panel is likely noisy. Epithelium is the main source in pregnancy and lactation. Macrophages are the main source in virgin and post-weaning (about 7 to 10%, against about 2% for epithelium).

Why: if PCs don’t express April receptors, co-localisation with April would matter less. I checked BCMA (Tnfrsf17, the main April receptor on PCs), TACI (Tnfrsf13b, which also binds April) and BAFF-R (Tnfrsf13c, which binds BAFF only and is mainly on B cells, used as a control).

rec <- intersect(c("Tnfrsf17", "Tnfrsf13b", "Tnfrsf13c"), names(df))
df$rec_group <- "Other"
df$rec_group[df$is_immune] <- "Other immune"
df$rec_group[(g("Ms4a1") > 0 | g("Cd19") > 0) & !df$is_PC] <- "B cell"
df$rec_group[df$is_PC] <- "PC"

kable(df |>
        filter(stage != "Virgin") |>
        group_by(stage, rec_group) |>
        summarise(n = n(), across(all_of(rec), ~ round(100 * mean(.x > 0), 1)),
                  .groups = "drop"),
      caption = "% of cells detecting each receptor")

% of cells detecting each receptor
stage	rec_group	n	Tnfrsf17	Tnfrsf13b	Tnfrsf13c
Pregnancy	B cell	948	0.6	15.4	15.9
Pregnancy	Other	128955	0.3	0.3	0.2
Pregnancy	Other immune	2651	2.6	5.2	0.9
Pregnancy	PC	18	38.9	11.1	11.1
Lactation	B cell	226	0.4	6.6	6.2
Lactation	Other	232278	0.2	1.8	0.1
Lactation	Other immune	1769	2.9	5.0	0.4
Lactation	PC	47	44.7	6.4	2.1
Post-weaning	B cell	19	0.0	10.5	0.0
Post-weaning	Other	69994	0.3	0.2	0.1
Post-weaning	Other immune	736	2.7	2.3	0.1
Post-weaning	PC	280	28.2	8.9	1.4

BCMA is one of the markers used to call PCs, which could inflate its rate in PCs. To check, I looked only at PCs that qualify through a different support marker:

other_support <- setdiff(pc_support_genes, "Tnfrsf17")
pc_other <- df$is_PC & rowSums(sapply(other_support, function(x) df[[x]] > 0)) >= 1
kable(data.frame(
  n_PC     = as.vector(table(droplevels(df$stage[pc_other]))),
  pct_BCMA = round(100 * tapply(df$Tnfrsf17[pc_other] > 0,
                                droplevels(df$stage[pc_other]), mean), 1)),
  caption = "BCMA in PCs called without relying on BCMA")

BCMA in PCs called without relying on BCMA
	n_PC	pct_BCMA
Pregnancy	18	38.9
Lactation	47	44.7
Post-weaning	269	25.3

Conclusion: BCMA is clearly PC-specific: about 25 to 45% of PCs against under 3% of other cells. This holds when BCMA is not used to call the PC. BAFF-R is high in B cells and low in PCs, the expected switch as B cells become PCs, which further supports the PC calls. Xenium captures only part of each cell’s mRNA, so the true BCMA+ fraction is likely higher. PCs can respond to April in every stage.

for (s in levels(df$stage)) {
  d <- df[df$stage == s, ]
  p <- ggplot() +
    geom_point(data = d, aes(x, y), colour = "grey92", size = 0.02) +          # all cells (tissue)
    geom_point(data = d[d$is_epi & !d$april, ], aes(x, y), colour = "grey45", size = 0.05) +
    geom_point(data = d[d$is_epi &  d$april, ], aes(x, y), colour = "#E69F00", size = 0.1) +
    geom_point(data = d[d$is_PC, ],             aes(x, y), colour = "#0072B2", size = 0.4) +
    coord_fixed() +
    labs(title = s,
         subtitle = "Light grey: all cells   Dark grey: April- epithelium   Orange: April+ epithelium   Blue: plasma cells") +
    theme_void()
  print(p)
}

PC annotation: marker-based calls (Jchain plus a second PC marker, no epithelial signal) give clean IgA PCs: 0 in virgin, 18 in pregnancy, 47 in lactation, 280 in post-weaning.
PC abundance: PCs peak after weaning (about 19 per mm², against about 1 per mm² in pregnancy and lactation).
No co-localisation with April+ epithelium: in any stage, by either baseline, at any radius, or with alternative definitions.
No April-rich niche: PCs are not surrounded by more April than other immune cells or random cells, once cell density is accounted for.
PCs do favour the epithelium after weaning: about 1.5 to 3× enriched next to epithelium and myoepithelium, and clustered together. During lactation they sit in the stroma.
PCs express BCMA, so they could respond to April, wherever it comes from.

Caveats

One section (one mouse) per stage. The p-values describe patterns within a section, not differences between animals. Stage comparisons need replicates.
Tnfsf13 is detected at low levels (1 count = April+), so April+ calls are noisy and absence of a signal is not strong evidence of absence.
Only post-weaning is well powered. Lactation and pregnancy results are indicative only.
Transcript spillover in Xenium can mix signals between touching cells. The PC rule was designed to limit this.

Possible next steps

Replicate sections per stage, to test stage differences properly.
If the 20 protein markers in obsm/IFQUANT include CD138 or IgA, use them to confirm PCs.
Ask why PCs seem to gather around regressing epithelium after weaning, independently of April.

sessionInfo()

## R version 4.5.3 (2026-03-11)
## Platform: x86_64-pc-linux-gnu
## Running under: Red Hat Enterprise Linux 9.6 (Plow)
## 
## Matrix products: default
## BLAS:   /stornext/System/data/software/rhel/9/base/tools/R/4.5.3/lib64/R/lib/libRblas.so 
## LAPACK: /stornext/System/data/software/rhel/9/base/tools/R/4.5.3/lib64/R/lib/libRlapack.so;  LAPACK version 3.12.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] knitr_1.52                  tidyr_1.3.2                
##  [3] dplyr_1.2.1                 ggplot2_4.0.3              
##  [5] RANN_2.6.3                  zellkonverter_1.20.1       
##  [7] SingleCellExperiment_1.32.0 SummarizedExperiment_1.40.0
##  [9] Biobase_2.70.0              GenomicRanges_1.62.1       
## [11] Seqinfo_1.0.0               IRanges_2.44.0             
## [13] S4Vectors_0.48.1            BiocGenerics_0.56.0        
## [15] generics_0.1.4              MatrixGenerics_1.22.0      
## [17] matrixStats_1.5.0          
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10         SparseArray_1.10.10 lattice_0.23-1     
##  [4] magrittr_2.0.5      digest_0.6.39       RColorBrewer_1.1-3 
##  [7] evaluate_1.0.5      grid_4.5.3          fastmap_1.2.0      
## [10] jsonlite_2.0.0      Matrix_1.7-6        purrr_1.2.2        
## [13] viridisLite_0.4.3   scales_1.4.0        jquerylib_0.1.4    
## [16] abind_1.4-8         cli_3.6.6           rlang_1.3.0        
## [19] XVector_0.50.0      withr_3.0.3         cachem_1.1.0       
## [22] DelayedArray_0.36.1 yaml_2.3.12         otel_0.2.0         
## [25] S4Arrays_1.10.1     tools_4.5.3         dir.expiry_1.18.0  
## [28] parallel_4.5.3      filelock_1.0.3      basilisk_1.22.0    
## [31] reticulate_1.47.0   vctrs_0.7.3         R6_2.6.1           
## [34] png_0.1-9           lifecycle_1.0.5     pkgconfig_2.0.3    
## [37] pillar_1.11.1       bslib_0.12.0        gtable_0.3.6       
## [40] glue_1.8.1          Rcpp_1.1.2          tidyselect_1.2.1   
## [43] tibble_3.3.1        xfun_0.61           dichromat_2.0-1    
## [46] rstudioapi_0.19.0   farver_2.1.2        htmltools_0.5.9    
## [49] labeling_0.4.3      rmarkdown_2.32      compiler_4.5.3     
## [52] S7_0.2.2

Post-weaning