Summary

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

  1. The plasma cells had to be re-annotated. The supplied predicted_cell_type labels contain no plasma cell type, and 85% of the PCs I called had been labelled “Endothelial cell”. I called PCs from markers instead, and checked that they are dominated by immunoglobulin transcripts (IgA PCs).
  2. The supplied endothelial labels are largely adipocytes. After weaning, the apparent 80% endothelial cells are mostly adipocytes by marker genes (Lipe, Adipoq, Plin1, Fabp4).
  3. 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.
  4. PCs do not preferentially sit next to April+ epithelium, and are not in April-richer neighbourhoods than other immune cells, in any stage. This held across radii, baselines and alternative cell definitions.
  5. PCs do have a spatial preference. In lactation they sit in the stroma, next to adipocytes and endothelium. After weaning they sit next to the remaining epithelial structures (2.8× enriched next to epithelium, 3.7× next to myoepithelium; median distance to epithelium 29 µm against 76 µm for random cells), often in small PC clusters.
  6. PCs express the April receptor BCMA (Tnfrsf17), 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. April is detected at low levels everywhere. Only post-weaning has enough PCs for well-powered spatial tests.

Parameters

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_update_20261028"

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)

1. Data

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: none (all on the panel).

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.

2. Re-annotating plasma cells

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

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

2a. Gate plot

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 epithelial threshold separates the two populations well. Most epithelial cells have no Jchain counts, and a few have 1 to 4. A handful of cells per stage have high Jchain together with low-moderate epithelial counts (3 to 5); the gate excludes these, so spillover appears limited. On the Jchain axis there is no clear gap: counts run continuously from 1 upwards, so the cutoff of 5 is a judgement call. Virgin has no cell with more than 1 Jchain count.

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

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

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")
High-Jchain, low-keratin cells by number of support markers detected
0 1 2 3 4 5
Virgin 0 0 0 0 0 0
Pregnancy 5 4 8 3 2 1
Lactation 2 9 13 14 9 2
Post-weaning 57 125 91 45 15 4 
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. Those with none (64 cells, 57 of them in post-weaning) are mostly labelled endothelial (58 of 64), with a few fibroblasts and B cells. They may be picking up Jchain from nearby PCs, or be PCs whose other markers weren’t captured. I excluded them to be cautious; section 6b shows that including them does not change the result.

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

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 Adipoq total_counts
Endo-labelled 0 0 0 13 0 150
Other 0 0 0 0 0 289
PC 10 16 17 8 0 198

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

  • Endothelial genes (Pecam1, Cdh5, Kdr) are detected in 2 to 6% of PCs, against 12 to 19% of endothelial-labelled cells. Notably, most endothelial-labelled cells don’t detect these genes either, while 94% detect Fabp4 (see 2d).
  • The median PC has 43 Jchain + Igkc + Igha counts, about a fifth of its transcripts. The median endothelial-labelled cell has none.
  • Igha is as high as Igkc, consistent with IgA PCs, the main PC type in the mammary gland.
  • Adipocyte genes are also common in PCs: Fabp4 in 88% (median 8 counts), Adipoq and Plin1 in about 30%. This is probably signal from neighbouring adipocytes in the fat pad, which can bleed into small cells in Xenium, but I have not tested this directly.

2d. The supplied labels

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: 85% of PCs had been labelled “Endothelial cell”, mostly “Endothelial cell_Fabp4 high”. Across the whole dataset, 36% of cells carry an endothelial label, which seems too many for mammary gland. In 2c, most endothelial-labelled cells detected Fabp4 (94%) and many detected Adipoq or Plin1, but few detected endothelial genes, so this class is probably largely adipocytes (tested directly in section 8). I don’t think predicted_cell_type should be relied on for immune and stromal cells. PCs, epithelium and immune cells use my marker-based calls throughout; sections 8 and 9 use the supplied labels only for cell types I did not re-annotate (e.g. fibroblasts, myoepithelium, macrophages), and say so.

3. Epithelium, April and immune cells

Definitions (all from raw counts):

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 (about 31,000 to 34,000 cells), 16 to 20% of it April+. Post-weaning has very little epithelium (2,438 cells, about 3% of the section, consistent with involution), and only 1.8% of it is April+, giving just 45 April+ epithelial cells, yet this is the stage with the most PCs. Using a looser epithelial definition (Krt18 + Epcam + Krt19 + Cdh1 ≥ 2) only raises post-weaning to 3,734 epithelial cells, so the epithelium really is sparse there. Virgin epithelium is also sparse (475 cells, about 1% of the section), but virgin has no PCs.

4. PC density by stage

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 (18.7 vs 1.1 per mm²; 0.9 in pregnancy), and absent in virgin. Per 10,000 cells gives a similar picture (about 20×). Per mm² is the fairest measure: per epithelial cell is inflated in post-weaning because the epithelium shrinks. This was unexpected, as mammary IgA PCs are usually thought to build up during pregnancy and lactation. Here they are most abundant after weaning, possibly because PCs persist while the epithelium regresses, although with one section per stage this needs confirming.

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

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

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:

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)

Conclusion: no convincing co-localisation in any stage.

6. Sensitivity checks

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.

So the weak post-weaning signal depends on the epithelial definition, and in every run PCs are no closer to April+ epithelium than other immune cells are.

7. Do PCs sit in April-rich neighbourhoods?

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:

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: once neighbourhood density is accounted for (april_per_nbr), PCs are not in April-richer surroundings than other immune cells in any stage (0.3 to 1.0×). There are a few isolated signals, each from one radius or one comparison:

The “other immune cells” reference set includes PCs (as defined in section 3), which makes these comparisons conservative.

P-values are uncorrected (96 comparisons).

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: in post-weaning, PCs have about 1.8× as many epithelial neighbours within 50 µm as other immune cells, which is more than enough to explain the 1.5× extra epithelial April. PCs are closer to epithelium, but the epithelium next to them is not more April+ (consistent with section 5, baseline 2).

8. What do PCs sit next to?

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. Following Julie’s feedback, the neighbour classes are now marker-based rather than taken from predicted_cell_type:

8a. Marker check by original label

extra_genes <- c("Lipe", "Fabp4", "Fasn", "Acly", "Adipoq", "Plin1",
                 "Pecam1", "Vwf", "Egfl7", "Itgb1", "Plvap", "Cdh5",
                 "Csn2", "Wap", "Lalba", "Glycam1")
not_on_panel <- setdiff(extra_genes, rownames(sce))
add <- setdiff(intersect(extra_genes, rownames(sce)), names(df))
if (length(add) > 0)
  df <- cbind(df, t(as.matrix(counts(sce)[add, , drop = FALSE])))
diag_genes <- intersect(extra_genes, names(df))

Markers not on the panel: Csn2, Lalba, Glycam1.

marker_diag <- df |>
  filter(stage %in% c("Lactation", "Post-weaning")) |>
  group_by(stage, predicted) |>
  summarise(n = n(),
            across(all_of(diag_genes), ~ round(100 * mean(.x > 0), 1)),
            .groups = "drop") |>
  filter(n >= 50) |>
  arrange(stage, desc(n))
write.csv(marker_diag, file.path(OUT_DIR, "marker_diagnostic_by_label.csv"), row.names = FALSE)
kable(marker_diag,
      caption = "% of cells detecting each marker, by original predicted label (labels with n ≥ 50)")
% of cells detecting each marker, by original predicted label (labels with n ≥ 50)
stage predicted n Lipe Fabp4 Fasn Acly Adipoq Plin1 Pecam1 Vwf Egfl7 Itgb1 Plvap Cdh5 Wap
Lactation Secretory alveoli cell_Trf high 128735 9.2 34.6 96.3 91.1 2.7 2.2 1.7 0.2 2.3 32.9 0.2 3.0 100.0
Lactation Endothelial cell_Fabp4 high 30032 43.4 98.5 77.6 61.9 40.1 35.7 26.4 3.0 32.1 33.5 4.1 34.0 93.7
Lactation Secretory alveoli cell_Csn2 high 15565 2.2 21.8 67.5 50.6 0.8 0.7 1.6 0.1 1.8 10.0 0.2 2.2 99.6
Lactation Myoepithelial cell 10614 4.4 47.5 56.2 39.8 2.5 2.4 6.0 1.5 7.9 32.9 1.1 7.6 96.7
Lactation Epithelial cell_Csn1s1 high 9262 10.5 24.4 95.4 91.6 5.9 5.9 0.4 0.2 0.6 27.9 0.1 1.1 99.5
Lactation Endothelial cell_Aqp1 high 8453 4.0 85.7 22.0 12.7 2.5 1.4 23.1 9.7 27.7 18.2 10.8 25.3 80.8
Lactation Fibroblast 6062 8.1 68.2 46.8 32.1 7.3 5.6 10.4 2.4 12.0 42.9 3.3 12.6 95.4
Lactation Macrophage 3060 9.1 56.4 60.2 42.2 7.4 5.5 4.8 0.7 6.1 30.8 0.8 7.1 98.9
Lactation Fibroblast_Mgp high 2162 4.3 53.9 34.6 20.8 5.0 3.7 1.9 0.4 3.1 15.3 0.5 3.2 89.8
Lactation Macrophage_Pf4 high 2118 12.3 49.0 37.2 34.4 13.6 9.2 2.3 0.4 2.5 42.6 0.8 3.4 89.3
Lactation Fibroblast_Pi16 high 1727 4.9 18.5 21.9 34.7 5.1 2.8 0.2 0.0 0.2 42.4 0.1 0.6 90.5
Lactation Luminal cell_Fxyd2 high 1574 6.2 20.0 87.2 78.9 4.9 4.0 0.1 0.3 0.1 4.8 0.0 0.3 95.8
Lactation Fibroblast_Inmt high 1371 8.2 39.9 32.7 18.4 10.4 6.5 2.0 1.1 2.3 24.7 1.0 3.2 73.2
Lactation Secretory alveoli cell_Wfdc3 high 1200 0.3 0.7 0.8 1.0 0.4 0.2 0.1 0.0 0.0 0.1 0.0 0.1 99.0
Lactation Fibroblast_Col3a1 high 1188 42.8 96.6 95.1 75.7 46.0 42.1 3.4 0.4 3.9 37.1 0.7 6.1 97.4
Lactation Secretory alveoli cell_Wap high 1082 0.2 3.0 10.7 9.8 0.0 0.0 0.0 0.0 0.0 0.3 0.0 0.1 97.4
Lactation Luminal cell_Wfdc12 high 800 1.9 17.8 46.0 31.6 2.2 0.8 0.2 0.0 0.4 3.2 0.0 0.4 93.8
Lactation Luminal cell_Csn3 high 786 1.5 25.3 49.6 34.9 0.5 0.4 1.5 0.3 2.0 11.8 0.1 1.3 98.5
Lactation NK cell_Gzma high 772 4.5 41.1 44.2 33.2 3.1 2.7 6.9 0.6 4.8 30.2 1.2 4.8 92.7
Lactation Proliferating T cell 572 0.7 14.5 16.1 8.7 0.9 0.7 0.0 0.2 0.2 1.9 0.0 0.2 78.5
Lactation Macrophage_C1qb high 519 0.2 4.0 6.6 3.5 0.0 0.2 0.0 0.0 0.2 1.0 0.0 0.0 97.3
Lactation Secretory alveoli cell_Hes1 high 482 0.2 1.5 16.2 11.8 0.4 0.0 0.0 0.0 0.0 0.4 0.0 0.0 95.0
Lactation NK cell_Gzmb high 424 1.2 5.4 65.1 55.4 0.0 0.5 0.2 0.0 0.0 0.9 0.0 0.0 97.4
Lactation Luminal cell_Calca high 420 16.0 46.9 90.2 72.6 8.3 8.8 1.4 36.4 3.1 46.0 0.2 1.9 100.0
Lactation Smooth muscle cell_Acta2 high 394 0.3 13.5 34.0 16.5 0.5 0.3 0.5 0.3 0.8 7.9 0.3 0.8 92.1
Lactation Macrophage_Lyz2 high 362 8.0 28.5 21.0 26.8 4.7 1.9 1.1 1.7 1.4 19.9 0.3 0.6 85.4
Lactation Myeloid cell 338 0.6 8.6 13.0 7.4 0.6 0.6 0.9 0.0 0.3 2.7 0.0 0.0 78.7
Lactation Mast cell 335 7.8 46.3 37.6 28.1 10.1 9.9 3.6 0.3 1.2 17.6 0.6 1.5 83.9
Lactation Cd8 positive T cell_Cd8b1 high 319 13.8 71.5 69.9 52.0 11.3 11.0 13.2 2.2 13.5 43.6 3.8 14.1 99.1
Lactation Secretory alveoli cell_Klf2 high 315 0.3 6.3 22.5 16.2 0.3 0.0 0.0 0.0 0.3 1.0 0.0 0.3 85.7
Lactation Secretory alveoli cell_Hspa1a high 291 0.0 7.2 3.4 4.5 0.3 0.0 0.0 0.3 0.0 1.0 0.0 0.0 88.0
Lactation Secretory alveoli cell_Gpx4 high 231 0.0 0.9 2.2 2.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.4 100.0
Lactation Secretory alveoli cell_Lgals1 high 228 0.4 3.1 6.6 6.1 0.4 0.0 0.0 0.0 0.4 0.0 0.0 0.0 93.9
Lactation Dendritic cell 227 1.3 14.1 41.4 30.0 0.0 0.0 0.0 0.0 0.4 0.4 0.4 0.0 95.2
Lactation Proliferating cell 217 8.8 49.8 75.6 57.6 8.8 5.1 0.5 0.0 0.5 3.7 0.0 0.5 81.1
Lactation Secretory alveoli cell_Naxd high 153 0.0 2.0 9.2 2.6 0.7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 96.1
Lactation Secretory alveoli cell_Tubb5 high 151 0.0 2.0 17.2 21.9 0.7 0.0 0.0 0.0 0.0 0.7 0.0 0.0 92.1
Lactation Secretory alveoli cell_Taf1d high 138 0.0 5.1 2.9 2.9 0.7 0.7 0.0 0.0 0.0 2.2 0.0 0.0 87.0
Lactation Secretory alveoli cell_Crip1 high 127 0.8 15.7 5.5 3.9 0.8 0.0 0.0 1.6 0.0 0.8 0.0 0.0 77.2
Lactation Dendritic cell_Fscn1 high 107 1.9 42.1 7.5 9.3 0.0 0.0 3.7 3.7 3.7 21.5 2.8 5.6 85.0
Lactation Secretory alveoli cell_Msrb1 high 103 0.0 1.0 2.9 1.9 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 100.0
Lactation Plasmacytoid dendritic cell_Siglech high 97 4.1 68.0 42.3 22.7 6.2 2.1 7.2 2.1 9.3 30.9 1.0 4.1 95.9
Lactation Secretory alveoli cell_Fcer1g high 95 0.0 3.2 4.2 3.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 90.5
Lactation Luminal cell_Krt19 high 89 0.0 5.6 70.8 59.6 0.0 0.0 0.0 0.0 1.1 11.2 0.0 1.1 100.0
Lactation Secretory alveoli cell_Gm1673 high 86 1.2 17.4 18.6 10.5 0.0 0.0 0.0 0.0 2.3 0.0 0.0 1.2 93.0
Lactation Dendritic cell_Cd209a high 78 0.0 15.4 23.1 9.0 2.6 1.3 0.0 0.0 0.0 5.1 0.0 1.3 91.0
Lactation Secretory alveoli cell_Tyrobp high 72 0.0 51.4 9.7 4.2 0.0 1.4 0.0 2.8 8.3 5.6 0.0 0.0 27.8
Lactation Secretory alveoli cell_Elovl5 high 71 0.0 33.8 8.5 5.6 0.0 0.0 4.2 1.4 8.5 1.4 1.4 5.6 66.2
Lactation Secretory alveoli cell_Spp1 high 65 0.0 15.4 61.5 52.3 1.5 0.0 0.0 0.0 0.0 9.2 0.0 1.5 100.0
Lactation Secretory alveoli cell_Ifi35 high 62 0.0 0.0 21.0 14.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 95.2
Lactation Macrophage_Apoe high 60 0.0 0.0 33.3 18.3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 66.7
Lactation Secretory alveoli cell_Snord118 high 56 0.0 5.4 21.4 16.1 0.0 0.0 0.0 0.0 0.0 1.8 0.0 0.0 91.1
Post-weaning Endothelial cell_Fabp4 high 40860 55.9 99.9 97.3 77.4 58.3 55.5 6.8 0.4 5.3 16.7 0.3 7.3 34.5
Post-weaning Endothelial cell_Aqp1 high 19192 12.3 77.1 50.1 26.2 17.5 12.1 3.4 1.2 2.3 7.0 0.8 2.3 17.7
Post-weaning Fibroblast_Col3a1 high 3042 35.4 98.3 89.2 59.0 42.8 38.8 3.6 0.9 2.9 25.5 1.2 3.3 37.0
Post-weaning Macrophage 1972 16.2 60.4 45.3 34.7 19.5 18.5 1.7 0.6 1.6 27.6 0.6 1.9 37.8
Post-weaning Myoepithelial cell 1778 5.0 32.9 20.5 14.6 3.2 2.9 0.7 2.2 1.3 18.2 0.0 1.1 23.5
Post-weaning Fibroblast_Pi16 high 1586 2.6 19.7 11.9 21.4 4.0 3.3 0.3 0.0 0.1 28.1 0.2 0.1 29.0
Post-weaning Fibroblast 1293 7.2 56.5 30.5 16.5 10.7 6.8 1.4 0.8 0.9 17.3 0.9 0.6 25.9
Post-weaning Secretory alveoli cell_Trf high 1149 9.7 49.3 37.9 29.6 4.6 6.4 0.8 5.0 1.5 25.5 0.1 1.0 39.1
Post-weaning Luminal cell_Krt19 high 119 0.8 18.5 8.4 16.8 0.8 2.5 0.0 5.0 0.8 20.2 0.0 0.0 34.5

Conclusion:

  • Post-weaning “Endothelial cell_Fabp4 high” (the largest class) is adipocyte-like: Fabp4 in nearly all cells and Lipe, Adipoq and Plin1 in 55 to 58%, but endothelial genes in only 5 to 7%. This explains the apparent excess of endothelial cells after weaning.
  • In lactation the same label is mixed, with both adipocyte (36 to 43%) and endothelial (26 to 34%) genes. “Endothelial cell_Aqp1 high” looks like genuine endothelium.
  • “Fibroblast_Col3a1 high” is also adipocyte-like in both stages.
  • Fasn, Acly and Fabp4 are detected in many lactating cell types, so they are not specific enough on their own. Lipe, Adipoq and Plin1 are largely restricted to adipocyte-like labels.
  • Wap is detected in most lactating cells of every type, showing substantial spillover of milk transcripts; low-count calls in lactation should be treated with caution.

8b. Revised cell classes

ADIPO_SPECIFIC <- intersect(c("Lipe", "Adipoq", "Plin1"), names(df))
ADIPO_GENES    <- intersect(c("Lipe", "Adipoq", "Plin1", "Fabp4"), names(df))
ENDO_GENES     <- intersect(c("Pecam1", "Vwf", "Egfl7", "Plvap", "Cdh5"), names(df))
ADIPO_MIN      <- 2   # min number of adipocyte markers detected (at least 1 specific)
ENDO_MIN       <- 2   # min number of endothelial markers detected
EPI_EXCL_MIN   <- 1   # Krt18 + Epcam counts that rule a cell out as adipocyte
EPI_LABELS     <- c("Secretory alveoli cell", "Epithelial cell", "Luminal cell")

df$broad      <- sub("_.*", "", df$predicted)
df$adipo_n    <- rowSums(as.matrix(df[ADIPO_GENES]) > 0)
df$adipo_spec <- rowSums(as.matrix(df[ADIPO_SPECIFIC]) > 0)
df$endo_n     <- rowSums(as.matrix(df[ENDO_GENES]) > 0)
epi_any       <- g("Krt18") + g("Epcam") >= EPI_EXCL_MIN

df$is_adipo <- df$adipo_n >= ADIPO_MIN & df$adipo_spec >= 1 & df$adipo_n > df$endo_n &
  !epi_any & !(df$broad %in% EPI_LABELS) & !df$is_immune
df$is_endo  <- df$endo_n >= ENDO_MIN & df$endo_n > df$adipo_n &
  !df$is_epi & !df$is_immune

df$nb_class <- df$broad
df$nb_class[df$broad %in% EPI_LABELS | df$is_epi] <- "Epithelial"
df$nb_class[grepl("Endothelial", df$nb_class)]    <- "Endothelial label, marker-negative"
df$nb_class[df$is_adipo]                          <- "Adipocyte"
df$nb_class[df$is_endo]                           <- "Endothelial"
df$nb_class[df$is_immune & !df$is_PC]             <- "Immune (Ptprc+)"
df$nb_class[df$is_PC]                             <- "Plasma cell"

class_counts <- as.data.frame.matrix(table(df$nb_class, df$stage))
write.csv(class_counts, file.path(OUT_DIR, "revised_class_counts_by_stage.csv"))
kable(class_counts, caption = "Number of cells in each revised class, by stage")
Number of cells in each revised class, by stage
Virgin Pregnancy Lactation Post-weaning
Adipocyte 29896 23703 18769 43084
B cell 0 626 11 0
Cd8 positive T cell 0 318 162 0
Dendritic cell 0 37 443 0
Endothelial 1380 8806 11222 691
Endothelial label, marker-negative 6394 10632 13663 18143
Epithelial 1239 67402 161937 2712
Fibroblast 2155 6051 9793 3469
Immune (Ptprc+) 614 2949 1797 739
Macrophage 545 2461 4453 1088
Mast cell 39 109 259 0
Myeloid cell 0 8 334 0
Myoepithelial cell 189 9216 9294 823
NK cell 0 206 906 0
Plasma cell 0 18 47 280
Plasmacytoid dendritic cell 0 15 63 0
Proliferating cell 0 1 169 0
Proliferating T cell 0 0 567 0
Smooth muscle cell 0 0 390 0
T cell 0 14 23 0

Conclusion:

  • Adipocytes dominate the fat-pad-rich sections: about 70% of cells in virgin and 61% in post-weaning, against 18% in pregnancy and 8% in lactation, where epithelium expands (51% and 69% of cells).
  • Post-weaning endothelium drops from the apparent 80% to about 1% (691 cells). This is likely an underestimate, as endothelial genes are rarely detected in this section (8a).
  • About a quarter of post-weaning cells (18,143) carry an endothelial label but meet neither the endothelial nor the adipocyte rule. These are left as a separate, unresolved class rather than assigned.
  • B cell, T cell, NK and other immune labels shown separately are only the Ptprc-negative cells carrying those labels; Ptprc+ cells are grouped as “Immune (Ptprc+)”.

8c. Neighbours within 20 µm of PCs

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]]
  near <- near[near != "Plasma cell"]
  d_np <- d[!d$is_PC, ]
  lv   <- sort(unique(d_np$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_np$nb_class, lv))) / nrow(d_np))
}) |>
  bind_rows() |>
  mutate(enrichment = pct_near_PC / pct_section) |>
  filter(pct_near_PC > 0) |>
  arrange(stage, desc(pct_near_PC))
write.csv(pc_neighbours, file.path(OUT_DIR, "pc_neighbours_20um.csv"), row.names = FALSE)

kable(pc_neighbours, digits = 2,
      caption = "Non-PC neighbours within 20 µm of PCs (enrichment > 1 = more common next to PCs)")
Non-PC neighbours within 20 µm of PCs (enrichment > 1 = more common next to PCs)
stage class pct_near_PC pct_section enrichment
Lactation Epithelial 56.48 69.12 0.82
Lactation Adipocyte 13.19 8.01 1.65
Lactation Endothelial label, marker-negative 10.11 5.83 1.73
Lactation Endothelial 7.03 4.79 1.47
Lactation Fibroblast 5.05 4.18 1.21
Lactation Myoepithelial cell 5.05 3.97 1.27
Lactation Immune (Ptprc+) 1.10 0.77 1.43
Lactation Macrophage 1.10 1.90 0.58
Lactation NK cell 0.44 0.39 1.14
Lactation Cd8 positive T cell 0.22 0.07 3.18
Lactation Proliferating T cell 0.22 0.24 0.91
Post-weaning Adipocyte 46.05 60.90 0.76
Post-weaning Endothelial label, marker-negative 30.74 25.64 1.20
Post-weaning Epithelial 10.57 3.83 2.76
Post-weaning Fibroblast 4.70 4.90 0.96
Post-weaning Myoepithelial cell 4.27 1.16 3.67
Post-weaning Immune (Ptprc+) 1.40 1.04 1.34
Post-weaning Macrophage 1.40 1.54 0.91
Post-weaning Endothelial 0.88 0.98 0.90 
pc_neighbours |>
  filter(pct_section >= 0.5 | pct_near_PC >= 2) |>
  ggplot(aes(log2(enrichment), reorder(class, enrichment))) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_point(aes(size = pct_near_PC), colour = "#0072B2") +
  facet_wrap(~stage) +
  labs(x = "log2 enrichment next to PCs (within 20 µm)", y = NULL,
       size = "% of PC\nneighbours")

Classes making up under 0.5% of the section and under 2% of PC neighbours are left out of the plot (they are still in the table).

Conclusion: PCs sit in the stroma during lactation and next to the epithelial structures after weaning.

  • Lactation: PCs are enriched next to adipocytes (1.7×), endothelial cells (1.5×), endothelial-labelled unresolved cells (1.7×), myoepithelium (1.3×) and fibroblasts (1.2×), and depleted next to epithelium (0.8×), which still makes up over half of their neighbours because it dominates the section.
  • Post-weaning: PCs are enriched next to myoepithelium (3.7×) and epithelium (2.8×), and depleted next to adipocytes (0.8×), although adipocytes are still their most common neighbour (46%). Endothelial cells and fibroblasts are at about background (0.9 to 1.0×).
  • These enrichments are descriptive (no statistical test), and classes with few neighbours are unreliable: in lactation, PCs have about 450 neighbours in total, so the immune, macrophage and T cell values rest on 1 to 5 cells. Section 8d tests distances directly for post-weaning.

8d. Post-weaning: distance from each PC to the nearest cell of each type

Why: this somewhat matches the confocal analysis (figure 6B). Each blue dot is one PC. The grey violins show the same distance measured from 5,000 random non-PC cells, so that abundant cell types don’t look close to PCs just because they are common. Bars are medians.

VIOLIN_CLASSES <- c("Epithelial", "Myoepithelial cell", "Adipocyte", "Endothelial",
                    "Fibroblast", "Macrophage", "Immune (Ptprc+)", "Plasma cell")

d      <- df[df$stage == "Post-weaning", ]
xy     <- as.matrix(d[, c("x", "y")])
pc_i   <- which(d$is_PC)
set.seed(SEED)
rand_i <- sample(which(!d$is_PC), min(N_RANDOM, sum(!d$is_PC)))

nearest <- function(target_i, query_i) {
  nn <- RANN::nn2(xy[target_i, , drop = FALSE], xy[query_i, , drop = FALSE], k = 2)$nn.dists
  ifelse(query_i %in% target_i, nn[, 2], nn[, 1])   # skip self
}

dist_df <- lapply(VIOLIN_CLASSES, function(cl) {
  t_i <- which(d$nb_class == cl)
  if (length(t_i) < 2) return(NULL)
  rbind(data.frame(class = cl, from = "Plasma cells", dist = nearest(t_i, pc_i)),
        data.frame(class = cl, from = "Random cells", dist = nearest(t_i, rand_i)))
}) |>
  bind_rows() |>
  mutate(class = factor(class, levels = VIOLIN_CLASSES))

p_dist <- ggplot(dist_df, aes(from, dist, fill = from)) +
  geom_violin(scale = "width", colour = NA, alpha = 0.5) +
  geom_jitter(data = filter(dist_df, from == "Plasma cells"),
              width = 0.15, size = 0.4) +
  stat_summary(fun = median, fun.min = median, fun.max = median,
               geom = "crossbar", width = 0.5, linewidth = 0.3) +
  scale_y_log10() +
  scale_fill_manual(values = c("Plasma cells" = "#0072B2", "Random cells" = "grey60")) +
  facet_wrap(~ class, nrow = 1) +
  labs(x = NULL, y = "Distance to nearest cell (µm, log scale)") +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1))
ggsave(file.path(OUT_DIR, "postweaning_pc_distance_violin.png"), p_dist,
       width = 13, height = 4.5)
p_dist

dist_stats <- dist_df |>
  group_by(class) |>
  summarise(n_PC          = sum(from == "Plasma cells"),
            median_PC     = median(dist[from == "Plasma cells"]),
            median_random = median(dist[from == "Random cells"]),
            p = wilcox.test(dist[from == "Plasma cells"],
                            dist[from == "Random cells"])$p.value,
            .groups = "drop") |>
    mutate(p_adj = p.adjust(p, "BH"),
         across(c(p, p_adj), ~ format.pval(.x, digits = 2, eps = 1e-4)))
write.csv(dist_stats, file.path(OUT_DIR, "postweaning_pc_distance_stats.csv"), row.names = FALSE)
kable(dist_stats, digits = 3,
      caption = "Median distance (µm) to the nearest cell of each type; Wilcoxon test, PCs vs random cells, BH-adjusted")
Median distance (µm) to the nearest cell of each type; Wilcoxon test, PCs vs random cells, BH-adjusted
class n_PC median_PC median_random p p_adj
Epithelial 280 28.544 75.523 < 1e-04 < 1e-04
Myoepithelial cell 280 45.463 89.167 < 1e-04 < 1e-04
Adipocyte 280 8.894 9.800 0.00084 0.00112
Endothelial 280 65.345 72.970 0.00019 0.00030
Fibroblast 280 35.609 43.386 < 1e-04 0.00018
Macrophage 280 70.611 78.703 0.00705 0.00806
Immune (Ptprc+) 280 64.698 69.685 0.07710 0.07710
Plasma cell 280 59.751 118.557 < 1e-04 < 1e-04

Conclusion: PCs are much closer than random cells to epithelial structures and to each other, and only marginally closer to other cell types.

  • Strongest differences: epithelium (median 29 vs 76 µm, about 2.6× closer), myoepithelium (45 vs 89 µm, 2×) and other PCs (60 vs 119 µm, 2×). PC-to-PC distances are bimodal: some PCs sit in tight clusters (another PC within about 10 µm), while the rest are 50 to 200 µm apart.
  • Small differences: fibroblasts (36 vs 43 µm), endothelium (65 vs 73 µm), macrophages (71 vs 79 µm) and adipocytes (8.9 vs 9.8 µm) are statistically significant but differ by only 1 to 8 µm, about one cell width or less. Other immune cells do not differ (65 vs 70 µm, p = 0.08).
  • With 280 PCs and 5,000 random cells, small shifts reach significance, so effect size is the better guide. Endothelial distances are likely inflated because few endothelial cells were called in post-weaning (8b). The random cells are about 60% adipocytes, so “closer than random” partly reflects PCs sitting outside the pure fat pad.

9. Which cells express April?

This uses the supplied predicted_cell_type labels (broad type), except for PCs, which use my calls. Only types with at least 50 cells are shown.

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: April is detected at low levels everywhere: at most 18% of cells of any type have a transcript, and the average is 0.2 counts per cell or less, so any April readout from this panel is noisy.

10. Can PCs respond to April?

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:

PCs can respond to April in every stage.

11. Spatial maps

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)
}

12. Overall conclusions

  1. PC annotation: marker-based calls (Jchain plus a second PC marker, no epithelial signal) give IgA PCs: 0 in virgin, 18 in pregnancy, 47 in lactation, 280 in post-weaning.
  2. Supplied labels: “Endothelial cell_Fabp4 high”, the largest class after weaning, is mostly adipocytes. Endothelial genes are rarely detected in that section.
  3. PC abundance: PCs are most abundant after weaning (about 19 per mm², against about 1 per mm² in pregnancy and lactation).
  4. No co-localisation with April+ epithelium: in any stage, against either baseline, at any radius, or with alternative definitions. A weak post-weaning signal against all epithelium (about 1.3×) disappears when compared with other immune cells.
  5. No April-rich niche: once neighbourhood density is accounted for, PCs are not surrounded by more April than other immune cells.
  6. PCs favour the epithelium after weaning: about 3× enriched next to epithelium and myoepithelium, and about 2.6× closer to epithelium than random cells. During lactation they sit in the stroma.
  7. PCs express BCMA, so they could respond to April, wherever it comes from.

13. Caveats and possible next steps?

Caveats

Possible next steps

Session info

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] gtable_0.3.6        dir.expiry_1.18.0   xfun_0.62          
##  [4] bslib_0.12.0        rhdf5_2.54.1        lattice_0.22-9     
##  [7] rhdf5filters_1.22.0 vctrs_0.7.3         tools_4.5.3        
## [10] parallel_4.5.3      tibble_3.3.1        pkgconfig_2.0.3    
## [13] Matrix_1.7-4        RColorBrewer_1.1-3  S7_0.2.2           
## [16] lifecycle_1.0.5     compiler_4.5.3      farver_2.1.2       
## [19] htmltools_0.5.9     sass_0.4.10         yaml_2.3.12        
## [22] pillar_1.11.1       jquerylib_0.1.4     DelayedArray_0.36.1
## [25] cachem_1.1.0        abind_1.4-8         basilisk_1.22.0    
## [28] tidyselect_1.2.1    digest_0.6.39       purrr_1.2.2        
## [31] labeling_0.4.3      fastmap_1.2.0       grid_4.5.3         
## [34] cli_3.6.6           SparseArray_1.10.10 magrittr_2.0.5     
## [37] S4Arrays_1.10.1     h5mread_1.2.1       withr_3.0.3        
## [40] filelock_1.0.3      scales_1.4.0        rmarkdown_2.32     
## [43] XVector_0.50.0      reticulate_1.47.0   png_0.1-9          
## [46] HDF5Array_1.38.0    evaluate_1.0.5      viridisLite_0.4.3  
## [49] rlang_1.3.0         Rcpp_1.1.2          glue_1.8.1         
## [52] jsonlite_2.0.0      R6_2.6.1            Rhdf5lib_1.32.0