Modelling cell-type proportions in repeated-measures designs with sccomp: a TMA case study

Contributor: Pierre-Paul Axisa, CBMR, Copenhagen, Denmark.

Introduction

Testing cell proportional differences from single-cell-resolved count data is central to modern tissue analyses. It becomes a repeated-measures problem as soon as several observations come from the same donor, such as multiple tumour cores in a tissue microarray (TMA) or serial biopsies over a time course. Averaging cores up to the patient and testing one phenotype at a time drops both the compositional constraint and the repeated-measures structure.

sccomp is a Bayesian model for mixed-effect compositional designs that identifies outliers probabilistically, tests for differential variability, and removes unwanted effects. It takes a long tibble with one row per sample and cell type, a count column, and the design columns, and sccomp_estimate(), sccomp_test(), and sccomp_predict() return tibbles ready for dplyr, tidyr, and ggplot2.

Source study and study design

We apply sccomp to the NeoTRIP trial. Wang et al. (2023) used imaging mass cytometry (IMC) to map how multicellular spatial organization in triple-negative breast cancer relates to neoadjuvant immunotherapy response. In the phase III NeoTRIP trial (NCT02620280), early high-risk TNBC patients were randomised 1:1 to carboplatin + nab-paclitaxel (C) or the same regimen plus atezolizumab (C&I) for eight 3-week neoadjuvant cycles before surgery. Research FFPE biopsies:

• Baseline (n = 243)
• On-treatment (first day of cycle 2; n = 207)
• Post-treatment (surgical excision; n = 210)
• Profiled with a 43-marker IMC panel at subcellular resolution
• Deep-learning segmentation and semi-supervised clustering defined 17 epithelial and 20 TME cell phenotypes
• Multiple regions of interest (ROIs) were imaged per biopsy

Each patient contributes several samples. The biopsies are nested in the patient, and the regions of interest are nested in the biopsy. Patient NT193, in the chemotherapy arm with residual disease, was imaged at all three timepoints: 10 ROIs at baseline, 10 on-treatment, and 3 at surgery. In the model, each ROI is one compositional sample and the patient is the repeated-measures group.

mk_spatial_tile <- function(dat, square = TRUE) {
  if (nrow(dat) == 0) {
    dat <- tibble::tibble(
      Location_Center_X = numeric(),
      Location_Center_Y = numeric(),
      draw = character()
    )
  }
  tile <- ggplot(dat, aes(x = Location_Center_X, y = Location_Center_Y)) +
    geom_point(data = filter(dat, draw == "other"), colour = "#39FF14", size = 0.10) +
    geom_point(data = filter(dat, draw == "pdl1"), colour = "#FFFF00", size = 0.16) +
    geom_point(data = filter(dat, draw == "m2"), colour = "#FF00FF", size = 0.16) +
    scale_x_continuous(labels = NULL, guide = "none", expand = c(0, 0)) +
    scale_y_reverse(labels = NULL, guide = "none", expand = c(0, 0)) +
    theme_void(base_size = 8) +
    theme(
      plot.background = element_rect(fill = "white", colour = NA),
      panel.background = element_rect(fill = "black", colour = NA),
      legend.position = "none",
      plot.margin = margin(
        if (square) 1 else 0, if (square) 1 else 0,
        if (square) 1 else 0, if (square) 1 else 0,
        unit = "pt"
      )
    )
  if (square) tile <- tile + theme(aspect.ratio = 1)
  tile
}

roi <- read_csv(
  file.path(results_dir, "sccomp_long_tme_immune_roi.csv"),
  col_select = c(PatientID, BiopsyPhase, ImageID),
  show_col_types = FALSE
) |>
  filter(PatientID == "NT193") |>
  distinct() |>
  mutate(BiopsyPhase = factor(BiopsyPhase, phase_levels)) |>
  arrange(BiopsyPhase, ImageID) |>
  mutate(i = row_number(), .by = BiopsyPhase)

n_phase <- roi |>
  summarise(n = n(), .by = BiopsyPhase) |>
  mutate(x0 = cumsum(n + 2) - (n + 2))
roi <- roi |>
  left_join(n_phase, by = "BiopsyPhase") |>
  mutate(x = i + x0)

cells_nt <- read_fst(file.path(post_dir, "nt193-cells.fst")) |>
  as_tibble()
dx <- 0.46
pad <- 0.62
xmin <- min(roi$x) - pad
xmax <- max(roi$x) + pad
fig_w <- 12
fig_h <- 5.4
target_in <- (2 * dx) / (xmax - xmin) * fig_w
headroom <- 1.9
tile_h <- target_in * headroom / (fig_h - target_in)

phase <- roi |>
  summarise(x = mean(x), .by = BiopsyPhase) |>
  mutate(
    y = tile_h + 0.62,
    label = case_match(
      as.character(BiopsyPhase),
      "Baseline" ~ "Baseline\n10 ROIs",
      "On-treatment" ~ "On-treatment\n10 ROIs",
      "Post-treatment" ~ "Post-treatment\n3 ROIs"
    )
  )
patient <- tibble::tibble(
  x = mean(phase$x), y = phase$y[1] + 0.95, label = "Donor"
)
edges <- bind_rows(
  phase |> mutate(x, y, xend = patient$x, yend = patient$y),
  roi |> mutate(
    x,
    y = tile_h,
    xend = phase$x[match(BiopsyPhase, phase$BiopsyPhase)],
    yend = phase$y[1]
  )
)
phase_cols <- c(
  Baseline = "#1B9E77",
  `On-treatment` = "#D95F02",
  `Post-treatment` = "#7570B3"
)

tree <- ggplot() +
  geom_segment(
    data = edges, aes(x = x, y = y, xend = xend, yend = yend),
    colour = "grey70", linewidth = 0.4
  ) +
  geom_label(
    data = phase,
    aes(x, y, label = label, fill = BiopsyPhase),
    colour = "white", fontface = "bold", size = 3, label.size = 0
  ) +
  geom_label(
    data = patient, aes(x, y, label = label),
    fill = "grey15", colour = "white", fontface = "bold", size = 3.6, label.size = 0
  ) +
  scale_fill_manual(values = phase_cols, guide = "none") +
  coord_cartesian(
    xlim = c(xmin, xmax),
    ylim = c(-0.04, patient$y + 0.4),
    clip = "off"
  ) +
  labs(title = "Patient NT193") +
  theme_void(base_size = 11) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5, size = 12),
    plot.margin = margin(6, 6, 6, 6)
  )

for (i in seq_len(nrow(roi))) {
  tree <- tree + annotation_custom(
    ggplotGrob(mk_spatial_tile(filter(cells_nt, ImageID == roi$ImageID[i]), square = FALSE)),
    xmin = roi$x[i] - dx,
    xmax = roi$x[i] + dx,
    ymin = 0,
    ymax = tile_h
  )
}
tree

Figure 1

Wang et al. linked cell densities, spatial interactions, and proliferative fractions to pathological complete response (pCR), and showed that immunotherapy remodels tissue composition over time, with response-specific dynamics of immune infiltration (notably in their Fig. 5 compositional trajectory analysis).

Quantity	Value
Per-protocol patients	257
ROI-level immune samples	1,740
Immune cell phenotypes	15
Median ROIs per patient	7

sccomp mixed-effects compositional modelling

Wang et al. summarised longitudinal composition as mean phenotype proportions per patient and compared arms, timepoints, and response groups with Wilcoxon tests, Z-scored trend plots (their Fig. 5), and per-phenotype weighted binomial GLMs on tumour-level proportions. Those tests treat each phenotype in isolation and average ROIs to the patient before inference, setting aside compositional closure (proportions sum to one) and the nesting of biopsies and ROIs within patients.

With sccomp, we model the ROI-level cell counts directly on a latent log-ratio scale with mixed effects on the mean and variability. Each ROI is one compositional sample (sample = "ROI_id" in sccomp_estimate()), several ROIs share a biopsy timepoint (BiopsyPhase), and a patient random intercept (1 | PatientID) accounts for baseline differences across individuals and shrinkage across the cohort:

\[\text{composition} \sim \text{Arm} \times \text{BiopsyPhase} \times \text{pCR} + (1 \mid \text{PatientID})\]

sccomp analysis

The blocks below are the minimal tidyverse code to reproduce the analysis locally from the parsed input (run top to bottom), using the immune TME compartment with ROI-level samples and a patient random intercept.

From long-format tidy input

The input long is a tibble with one row per ROI × immune phenotype: ROI_id, a count, the design columns (Arm, BiopsyPhase, pCR), and the grouping PatientID. Several ROIs share one BiopsyPhase. Preview:

library(readr)
library(dplyr)

long <- read_csv(file.path(results_dir, "sccomp_long_tme_immune_roi.csv"),
                 show_col_types = FALSE) |>
  rename(ROI_id = sample) |>
  mutate(
    count = as.integer(count),
    BiopsyPhase = factor(BiopsyPhase, c("Baseline", "On-treatment", "Post-treatment")),
    Arm = factor(Arm, c("C", "C&I")),
    pCR = factor(pCR, c("RD", "pCR"))
  )

long |>
  select(ROI_id, cell_group, count, PatientID, BiopsyPhase, Arm, pCR) |>
  gt::gt_preview() |>
  gt::tab_caption("Long-format sccomp input (1,740 ROIs x 15 immune phenotypes)")

Long-format sccomp input (1,740 ROIs x 15 immune phenotypes)

	ROI_id	cell_group	count	PatientID	BiopsyPhase	Arm	pCR
1	NTImg0011	CD8^+T	19	NT031	Baseline	C&I	RD
2	NTImg0011	CD8^+GZMB^+T	1	NT031	Baseline	C&I	RD
3	NTImg0011	CD8^+TCF1^+T	1	NT031	Baseline	C&I	RD
4	NTImg0011	CD8^+PD1^+T_{Ex}	7	NT031	Baseline	C&I	RD
5	NTImg0011	CD4^+TCF1^+T	4	NT031	Baseline	C&I	RD
6..26099
26100	NTImg1889	PD-L1^+IDO^+APCs	71	NT091	On-treatment	C&I	RD

Fit the model

library(sccomp)

fit <- sccomp_estimate(
  long,
  formula_composition = ~ Arm * BiopsyPhase * pCR + (1 | PatientID),
  formula_variability = ~ 1,
  sample = "ROI_id",
  cell_group = "cell_group",
  abundance = "count",
  inference_method = "hmc",
  cores = parallel::detectCores(),
  verbose = FALSE
) |>
  sccomp_test()

library(sccomp)

fit <- sccomp_estimate(
  long,
  formula_composition = ~ Arm * BiopsyPhase * pCR + (1 | PatientID),
  formula_variability = ~ 1,
  sample = "ROI_id",
  cell_group = "cell_group",
  abundance = "count",
  inference_method = "hmc",
  cores = parallel::detectCores(),
  verbose = FALSE
) |>
  sccomp_test()

sccomp model
============

Model specifications:
  Family: multi_beta_binomial 
  Composition formula: ~Arm * BiopsyPhase * pCR + (1 | PatientID) 
  Variability formula: ~1 
  Inference method: hmc 

Data: Samples: 1740   Cell groups: 15 

Column prefixes: c_ -> composition parameters  v_ -> variability parameters

Convergence diagnostics:
  For each parameter, n_eff is the effective sample size and R_k_hat is the potential
  scale reduction factor on split chains (at convergence, R_k_hat = 1).

# A tibble: 4,035 × 19
   cell_group parameter  factor c_lower c_effect  c_upper   c_pH0   c_FDR c_rhat
   <chr>      <chr>      <chr>    <dbl>    <dbl>    <dbl>   <dbl>   <dbl>  <dbl>
 1 CD20^+B    (Intercep… <NA>   -0.710   -0.514  -0.314   0       0         1.00
 2 CD20^+B    ArmC&I     Arm    -0.579   -0.292  -0.00208 0.0911  3.66e-2   1.00
 3 CD20^+B    BiopsyPha… Biops… -0.0970   0.134   0.359   0.377   1.30e-1   1.00
 4 CD20^+B    BiopsyPha… Biops…  0.903    1.12    1.33    0       0         1.00
 5 CD20^+B    pCRpCR     pCR    -0.162    0.112   0.390   0.465   2.21e-1   1.00
 6 CD20^+B    ArmC&I:Bi… <NA>   -0.281    0.0545  0.396   0.602   3.17e-1   1.00
 7 CD20^+B    ArmC&I:Bi… <NA>   -0.479   -0.151   0.189   0.384   1.24e-1   1.00
 8 CD20^+B    ArmC&I:pC… <NA>   -0.227    0.164   0.540   0.368   2.55e-1   1.00
 9 CD20^+B    BiopsyPha… <NA>   -0.274    0.0505  0.377   0.625   3.86e-1   1.00
10 CD20^+B    BiopsyPha… <NA>   -0.897   -0.572  -0.259   0.00175 4.17e-4   1.00
# ℹ 4,025 more rows
# ℹ 10 more variables: c_ess_bulk <dbl>, c_ess_tail <dbl>, v_lower <dbl>,
#   v_effect <dbl>, v_upper <dbl>, v_pH0 <dbl>, v_FDR <dbl>, v_rhat <dbl>,
#   v_ess_bulk <dbl>, v_ess_tail <dbl>

sccomp_proportional_fold_change() reports how fitted proportions change between two conditions. The result is descriptive; significance still comes from sccomp_test(). Here, chemoimmunotherapy at surgery, residual disease versus pathological complete response:

fit |>
  sccomp_proportional_fold_change(
    formula_composition = ~ Arm * BiopsyPhase * pCR,
    from = "C&I:Post-treatment:RD",  # chemoimmunotherapy, surgery, residual disease
    to = "C&I:Post-treatment:pCR"    # same arm and visit, pathological complete response
  ) |>
  gt::gt_preview() |>
  gt::tab_caption(
    "Proportional fold change: C&I post-treatment, pCR versus residual disease"
  )

Running standalone generated quantities after 1 MCMC chain, with 1 thread(s) per chain...

Chain 1  Elapsed Time: 0.1 seconds (Generated Quantities) 
Chain 1 finished in 0.0 seconds.

Proportional fold change: C&I post-treatment, pCR versus residual disease

	cell_group	proportion_fold_change	average_uncertainty	statement
1	CD20^+B	-1.300258	0.02649652	1.3-fold decrease (from 0.0645 to 0.0496)
2	CD4^+PD1^+T	-1.595534	0.04749542	1.6-fold decrease (from 0.0439 to 0.0275)
3	CD4^+TCF1^+T	-1.100683	0.01633892	1.1-fold decrease (from 0.0918 to 0.0834)
4	CD56^+NK	1.571215	0.04918382	1.6-fold increase (from 0.0808 to 0.1269)
5	CD79a^+Plasma	1.132715	0.02740352	1.1-fold increase (from 0.0721 to 0.0816)
6..14
15	Treg	-3.406119	0.25870011	3.4-fold decrease (from 0.0275 to 0.0081)

Model checks

Convergence. All 4,035 composition parameters have median \(\hat{R}\) = 1.001 (maximum 1.014). Bulk effective sample size has median 4,837 (minimum 495). Tail effective sample size has median 2,850 (minimum 923).

Variability is an intercept for each phenotype. Those 15 coefficients have maximum \(\hat{R}\) = 1.003, minimum bulk ESS 3,404, and minimum tail ESS 2,451. The composition parameters with the smallest bulk ESS:

fit |>
  filter(is.finite(c_ess_bulk)) |>
  slice_min(c_ess_bulk, n = 6, with_ties = FALSE) |>
  gt::gt_preview() |>
  gt::tab_caption(
    "Composition parameters with the smallest bulk effective sample size"
  )

Composition parameters with the smallest bulk effective sample size

	cell_group	parameter	factor	c_lower	c_effect	c_upper	c_pH0	c_FDR	c_rhat	c_ess_bulk	c_ess_tail	v_lower	v_effect	v_upper	v_pH0	v_FDR	v_rhat	v_ess_bulk	v_ess_tail
1	DCs	ArmC&I	Arm	-0.08375482	0.1601228	0.399447699	0.31431431	0.148577149	1.009882	600.7795	1411.607	NA	NA	NA	NA	NA	NA	NA	NA
2	DCs	(Intercept)	NA	0.77417609	0.9424418	1.116479608	0.00000000	0.000000000	1.011012	697.7882	1552.972	-3.78283	-2.957539	-2.134048	0	0	1.001649	4598.984	2582.758
3	DCs	ArmC&I:pCRpCR	NA	-0.47084651	-0.1203059	0.229162709	0.45070070	0.295827077	1.010088	717.1771	1444.238	NA	NA	NA	NA	NA	NA	NA	NA
4	M2 Mac	ArmC&I	Arm	0.12079694	0.3597968	0.601198700	0.01701702	0.009384384	1.019933	776.9721	1483.882	NA	NA	NA	NA	NA	NA	NA	NA
5	DCs	pCRpCR	pCR	-0.52733358	-0.2649256	-0.008725984	0.10010010	0.057807808	1.007276	858.1370	1509.393	NA	NA	NA	NA	NA	NA	NA	NA
6	M2 Mac	ArmC&I:pCRpCR	NA	-0.74712725	-0.4227944	-0.071407724	0.03678679	0.036786787	1.008612	878.2435	1583.149	NA	NA	NA	NA	NA	NA	NA	NA

Posterior predictive check. remove_unwanted_effects = TRUE removes variation other than biopsy timepoint, so the black boxes are adjusted proportions and the blue boxes are posterior predictive replicates. Fill marks a timepoint contrast with FDR below 0.05.

fit |> 
  sccomp_boxplot(
  factor = "BiopsyPhase",
  significance_statistic = "FDR",
  remove_unwanted_effects = TRUE
  ) +
  labs(
    title = "Posterior predictive check",
    subtitle = "Black: adjusted proportions. Blue: replicates from the fitted model.",
    x = "Biopsy timepoint"
  ) +
  guides(colour = "none", shape = "none") 

Hypothesis testing

We pass contrasts to sccomp_test() as strings.

marginal_CI_post compares responders with residual disease at that visit alone. Any gap the two groups already had at baseline is still in this result.

did_CI_post compares how far each group moved from its own baseline. A head start cancels, so the contrast is the difference in change. We use this contrast because a visit-only comparison still carries the baseline gap. Post-treatment M2 proportions are similar in the two response groups, so that snapshot does not describe the change.

contrasts <- c(
  marginal_CI_post = paste(
    "`pCRpCR`",                                  # pCR versus RD in chemotherapy at baseline
    "`ArmC&I:pCRpCR`",                           # additional pCR versus RD gap in the C&I arm
    "`BiopsyPhasePost-treatment:pCRpCR`",        # additional pCR versus RD gap at post-treatment
    "`ArmC&I:BiopsyPhasePost-treatment:pCRpCR`", # additional pCR versus RD gap for C&I at post-treatment
    sep = " + "
  ),
  did_CI_post = paste(
    "`BiopsyPhasePost-treatment:pCRpCR`",        # change from baseline in the pCR contrast
    "`ArmC&I:BiopsyPhasePost-treatment:pCRpCR`", # further change from baseline in the C&I arm
    sep = " + "
  )
)

res <- sccomp_test(fit, contrasts = contrasts)

res holds those two contrasts, one row per phenotype. The table is the full set of difference-in-differences, strongest first:

hmc_did |>
  filter(contrast_family == "did") |>
  arrange(c_FDR) |>
  gt::gt_preview() |>
  gt::tab_caption("Change from baseline (strongest by FDR)")

Change from baseline (strongest by FDR)

	parameter	cell_group	factor	c_lower	c_effect	c_upper	c_pH0	c_FDR	c_rhat	c_ess_bulk	c_ess_tail	contrast_family	Arm	BiopsyPhase	pCR
1	did__C&I__Post_treatment	CD56^+NK	NA	2.4152900	3.25669670	4.1111366	0.0000000	0.0000000	1.0010947	4247.316	3275.396	did	C&I	Post-treatment	pCR_minus_RD
2	did__C&I__Post_treatment	DCs	NA	0.6045030	1.33237959	2.0188092	0.0000000	0.0000000	0.9996942	5122.662	3061.238	did	C&I	Post-treatment	pCR_minus_RD
3	did__C&I__Post_treatment	Neutrophils	NA	1.6664088	2.63412173	3.5857767	0.0000000	0.0000000	1.0017040	4446.958	3336.609	did	C&I	Post-treatment	pCR_minus_RD
4	did__C&I__Post_treatment	PD-L1^+APCs	NA	-3.9597930	-2.80065090	-1.6193615	0.0000000	0.0000000	1.0002327	3677.553	2867.557	did	C&I	Post-treatment	pCR_minus_RD
5	did__C&I__Post_treatment	Treg	NA	-2.8820452	-1.94723055	-1.0560065	0.0000000	0.0000000	1.0006019	2449.355	2366.402	did	C&I	Post-treatment	pCR_minus_RD
6..59
60	did__C&I__On_treatment	CD8^+TCF1^+T	NA	-0.3600534	0.01963656	0.4067838	0.6599099	0.2803971	1.0012927	5303.619	3036.666	did	C&I	On-treatment	pCR_minus_RD

Population and per-patient predictions

# One row per arm, timepoint, and response group.
pop <- tibble::tribble(
  ~Arm,  ~BiopsyPhase,     ~pCR,
  "C",   "Baseline",       "RD",
  "C",   "Baseline",       "pCR",
  "C",   "On-treatment",   "RD",
  "C",   "On-treatment",   "pCR",
  "C",   "Post-treatment", "RD",
  "C",   "Post-treatment", "pCR",
  "C&I", "Baseline",       "RD",
  "C&I", "Baseline",       "pCR",
  "C&I", "On-treatment",   "RD",
  "C&I", "On-treatment",   "pCR",
  "C&I", "Post-treatment", "RD",
  "C&I", "Post-treatment", "pCR"
) |>
  mutate(
    Arm = factor(Arm, c("C", "C&I")),
    BiopsyPhase = factor(BiopsyPhase, c("Baseline", "On-treatment", "Post-treatment")),
    pCR = factor(pCR, c("RD", "pCR")),
    ROI_id = paste0("pred_", row_number()),
    PatientID = "NT031"
  )

gt::gt_preview(pop) |>
  gt::tab_caption("Population grid passed to sccomp_predict()")

Population grid passed to sccomp_predict()

	Arm	BiopsyPhase	pCR	ROI_id	PatientID
1	C	Baseline	RD	pred_1	NT031
2	C	Baseline	pCR	pred_2	NT031
3	C	On-treatment	RD	pred_3	NT031
4	C	On-treatment	pCR	pred_4	NT031
5	C	Post-treatment	RD	pred_5	NT031
6..11
12	C&I	Post-treatment	pCR	pred_12	NT031

pred_pop <- sccomp_predict(fit, new_data = pop)

Running standalone generated quantities after 1 MCMC chain, with 1 thread(s) per chain...

Chain 1  Elapsed Time: 0.185 seconds (Generated Quantities) 
Chain 1 finished in 0.0 seconds.

pred_pop |>
  gt::gt_preview() |>
  gt::tab_caption("Population predicted proportions")

Population predicted proportions

	Arm	BiopsyPhase	pCR	ROI_id	PatientID	cell_group	proportion_mean	proportion_lower	proportion_upper	unconstrained_mean	unconstrained_lower	unconstrained_upper
1	C	Baseline	RD	pred_1	NT031	CD20^+B	0.025475434	0.011849041	0.04656918	-0.78916098	-1.44397950	-0.18577897
2	C	Baseline	RD	pred_1	NT031	CD4^+PD1^+T	0.050589704	0.032262354	0.07256136	-0.07387369	-0.48652249	0.30169921
3	C	Baseline	RD	pred_1	NT031	CD4^+TCF1^+T	0.081006649	0.051677388	0.11175873	0.39988733	-0.02541777	0.76499292
4	C	Baseline	RD	pred_1	NT031	CD56^+NK	0.032822597	0.017883680	0.05537661	-0.52534173	-1.06513952	0.01087662
5	C	Baseline	RD	pred_1	NT031	CD79a^+Plasma	0.158685350	0.102258774	0.22280513	1.07227269	0.63309480	1.48362000
6..179
180	C&I	Post-treatment	pCR	pred_12	NT031	Treg	0.008723248	0.005836047	0.01223340	-1.53308238	-1.88980277	-1.19372516

# Per-patient grid: one row per patient and timepoint (uses the patient random intercept)
pat <- long |>
  distinct(PatientID, BiopsyPhase, Arm, pCR) |>
  mutate(ROI_id = paste0(PatientID, "_", BiopsyPhase))

pat |>
  gt::gt_preview() |>
  gt::tab_caption("Per-patient prediction grid (624 rows)")

Per-patient prediction grid (624 rows)

	PatientID	BiopsyPhase	Arm	pCR	ROI_id
1	NT031	Baseline	C&I	RD	NT031_Baseline
2	NT277	Baseline	C&I	RD	NT277_Baseline
3	NT277	On-treatment	C&I	RD	NT277_On-treatment
4	NT238	Baseline	C&I	RD	NT238_Baseline
5	NT238	On-treatment	C&I	RD	NT238_On-treatment
6..623
624	NT091	On-treatment	C&I	RD	NT091_On-treatment

pred_pat <- sccomp_predict(fit, new_data = pat)

Running standalone generated quantities after 1 MCMC chain, with 1 thread(s) per chain...

Chain 1  Elapsed Time: 5.094 seconds (Generated Quantities) 
Chain 1 finished in 0.0 seconds.

pred_pat

# A tibble: 9,360 × 12
   PatientID BiopsyPhase Arm   pCR   ROI_id         cell_group   proportion_mean
   <chr>     <fct>       <fct> <fct> <chr>          <chr>                  <dbl>
 1 NT031     Baseline    C&I   RD    NT031_Baseline CD20^+B               0.0175
 2 NT031     Baseline    C&I   RD    NT031_Baseline CD4^+PD1^+T           0.0427
 3 NT031     Baseline    C&I   RD    NT031_Baseline CD4^+TCF1^+T          0.0848
 4 NT031     Baseline    C&I   RD    NT031_Baseline CD56^+NK              0.0353
 5 NT031     Baseline    C&I   RD    NT031_Baseline CD79a^+Plas…          0.139 
 6 NT031     Baseline    C&I   RD    NT031_Baseline CD8^+GZMB^+T          0.0272
 7 NT031     Baseline    C&I   RD    NT031_Baseline CD8^+PD1^+T…          0.0420
 8 NT031     Baseline    C&I   RD    NT031_Baseline CD8^+T                0.0972
 9 NT031     Baseline    C&I   RD    NT031_Baseline CD8^+TCF1^+T          0.0477
10 NT031     Baseline    C&I   RD    NT031_Baseline DCs                   0.123 
# ℹ 9,350 more rows
# ℹ 5 more variables: proportion_lower <dbl>, proportion_upper <dbl>,
#   unconstrained_mean <dbl>, unconstrained_lower <dbl>,
#   unconstrained_upper <dbl>

Results

Predicted immune compositional trajectories

Figure 1 shows baseline-aligned predicted proportions (95% credible intervals) for all 15 immune TME phenotypes. Within each arm × cell-type facet, both response groups are shifted to the mean of their baseline predictions. Panels are split by treatment arm (columns); colour denotes pathological response. Asterisks mark arm–time points with difference-in-differences FDR < 0.05.

# Plot shift. Within each arm and cell type, move pCR and RD to the mean of
# their baseline predictions, on the posterior draws.
draws <- sccomp_predict(
  fit, new_data = pop, number_of_draws = 500, summary_instead_of_draws = FALSE
) |>
  as_tibble() |>
  mutate(across(c(Arm, BiopsyPhase, pCR), as.character))

Running standalone generated quantities after 1 MCMC chain, with 1 thread(s) per chain...

Chain 1  Elapsed Time: 0.183 seconds (Generated Quantities) 
Chain 1 finished in 0.0 seconds.

# Shift pCR and RD to their average baseline within each Arm x cell_group, per draw
shifts <- draws |>
  filter(BiopsyPhase == "Baseline") |>
  select(Arm, cell_group, .draw, pCR, proportion) |>
  pivot_wider(names_from = pCR, values_from = proportion) |>
  mutate(shift_pCR = (pCR + RD) / 2 - pCR,
         shift_RD  = (pCR + RD) / 2 - RD) |>
  select(Arm, cell_group, .draw, shift_pCR, shift_RD)

aligned_summary <- draws |>
  left_join(shifts, by = c("Arm", "cell_group", ".draw")) |>
  mutate(aligned_prop = proportion + if_else(pCR == "pCR", shift_pCR, shift_RD)) |>
  summarise(
    aligned_mean  = mean(aligned_prop),
    aligned_lower = quantile(aligned_prop, 0.025),
    aligned_upper = quantile(aligned_prop, 0.975),
    .by = c(Arm, BiopsyPhase, pCR, cell_group)
  )

cell_order <- c(
  "CD8^+T", "CD8^+GZMB^+T", "CD8^+TCF1^+T", "CD8^+PD1^+T_{Ex}",
  "CD4^+TCF1^+T", "CD4^+PD1^+T", "Treg", "CD20^+B", "CD79a^+Plasma",
  "CD56^+NK", "M2 Mac", "DCs", "Neutrophils", "PD-L1^+APCs", "PD-L1^+IDO^+APCs"
)
cols_response <- c(RD = "#377EB8", pCR = "#E41A1C")

pred_main <- aligned_summary |>
  mutate(
    BiopsyPhase = factor(BiopsyPhase, phase_levels),
    Arm = factor(Arm, arm_levels),
    pCR = factor(pCR, c("RD", "pCR")),
    cell_group = factor(cell_group, levels = cell_order),
    # A few intervals sit just below zero after the baseline shift.
    aligned_mean = pmax(aligned_mean, 0),
    aligned_lower = pmax(aligned_lower, 0)
  )

# Asterisk just above the ribbon where the DiD FDR is below 0.05.
did_sig_fig1 <- hmc_did |>
  filter(contrast_family == "did", is.finite(c_FDR), c_FDR < 0.05) |>
  distinct(cell_group, Arm, BiopsyPhase) |>
  mutate(cell_group = factor(cell_group, levels = cell_order)) |>
  inner_join(
    pred_main |>
      summarise(ymax = max(aligned_upper), .by = c(cell_group, Arm, BiopsyPhase)),
    by = c("cell_group", "Arm", "BiopsyPhase")
  )

ggplot(pred_main, aes(x = BiopsyPhase, y = aligned_mean, colour = pCR, group = pCR)) +
  geom_ribbon(
    aes(ymin = aligned_lower, ymax = aligned_upper, fill = pCR),
    alpha = 0.10, colour = NA
  ) +
  geom_line(linewidth = 0.75) +
  geom_point(size = 1.2) +
  geom_text(
    data = did_sig_fig1,
    aes(x = BiopsyPhase, y = ymax, label = "*"),
    inherit.aes = FALSE,
    size = 4.5, colour = "black", fontface = "bold", vjust = -0.2
  ) +
  facet_grid(cell_group ~ Arm, scales = "free_y") +
  scale_y_continuous(labels = scales::label_number(accuracy = 0.01)) +
  coord_cartesian(clip = "off") +
  scale_colour_manual(values = cols_response, name = "Response") +
  scale_fill_manual(values = cols_response, name = "Response") +
  labs(
    title = "Baseline-aligned immune compositional trajectories",
    x = NULL,
    y = "Baseline-aligned predicted proportion"
  ) +
  theme_pub() +
  theme(strip.text.y = element_text(size = 7))

Figure 2

Change from baseline

For each cell type and treatment arm, we test whether responders and patients with residual disease changed by different amounts from baseline. Figure 2 shows that difference at the on-treatment and post-treatment visits. Bordered tiles have FDR < 0.05.

did_hm <- hmc_did |>
  filter(contrast_family == "did") |>
  mutate(
    cell_group = factor(cell_group, levels = rev(cell_order)),
    BiopsyPhase = factor(BiopsyPhase, levels = c("On-treatment", "Post-treatment")),
    Arm = factor(Arm, arm_levels)
  ) |>
  filter(!is.na(cell_group)) |>
  mutate(
    sig = is.finite(c_FDR) & c_FDR < 0.05,
    facet_lab = factor(
      paste(Arm, BiopsyPhase, sep = " · "),
      levels = c(
        "C · On-treatment", "C · Post-treatment",
        "C&I · On-treatment", "C&I · Post-treatment"
      )
    )
  )

lim <- max(abs(did_hm$c_effect), na.rm = TRUE)

ggplot(did_hm, aes(x = facet_lab, y = cell_group, fill = c_effect)) +
  geom_tile(colour = "grey90", linewidth = 0.3) +
  geom_tile(
    data = filter(did_hm, sig),
    colour = "black", linewidth = 0.7, fill = NA
  ) +
  geom_text(
    aes(label = if_else(sig, sprintf("%.2f", c_effect), "")),
    size = 2.6
  ) +
  scale_fill_gradient2(
    low = cols_response[["RD"]], mid = "white", high = cols_response[["pCR"]],
    midpoint = 0, limits = c(-lim, lim), name = "center_log_ratio"
  ) +
  labs(
    title = "Response-associated change from baseline (difference-in-differences)",
    x = NULL, y = NULL
  ) +
  theme_pub() +
  theme(
    axis.text.x = element_text(angle = 30, hjust = 1),
    panel.grid = element_blank()
  )

Figure 3

Finding: Under the difference-in-differences contrast, 14 immune phenotypes change from baseline at post-treatment, and 9 at the on-treatment visit.

Finding: In the chemoimmunotherapy arm at post-treatment, M2 macrophages increase more from baseline in patients who reach pathological complete response (FDR 0.007), while PD-L1-positive antigen-presenting cells decrease (FDR < 0.0002).

Finding: The M2 macrophage change is not detected at the on-treatment visit (FDR 0.21).

M2 macrophage changes by response

Figure 3 compares pCR versus residual disease at each visit with the change from baseline. Points above the dotted line are FDR < 0.05.

m2_volcano <- bind_rows(
  hmc_marginal_pcr |>
    filter(cell_group == "M2 Mac") |>
    mutate(Arm, BiopsyPhase, c_effect, c_FDR, family = "Marginal pCR vs RD"),
  hmc_did |>
    filter(contrast_family == "did", cell_group == "M2 Mac") |>
    mutate(Arm, BiopsyPhase, c_effect, c_FDR, family = "DiD (change from T0)")
) |>
  mutate(label = paste(Arm, BiopsyPhase, sep = " · "))

ggplot(m2_volcano, aes(x = c_effect, y = -log10(pmax(c_FDR, 1e-12)))) +
  geom_vline(xintercept = 0, linetype = 2, colour = "grey70") +
  geom_hline(yintercept = -log10(0.05), linetype = 2, colour = "grey70") +
  geom_point(aes(colour = family, shape = BiopsyPhase), size = 2.8) +
  geom_text(aes(label = label), size = 2.3, vjust = -0.9, colour = "grey25") +
  scale_colour_manual(values = c(
    "Marginal pCR vs RD" = "#999999",
    "DiD (change from T0)" = "#d73027"
  )) +
  labs(
    title = "M2 Mac: visit contrast and change from baseline",
    x = "center_log_ratio",
    y = expression(-log[10](FDR)),
    colour = NULL, shape = NULL
  ) +
  theme_pub()

Figure 4

Between-patient heterogeneity in M2 Mac

The fixed effect is the population estimate. Patient random intercepts describe individual baselines and trajectories. Figure 4 shows the per-patient predicted M2 Mac trajectories.

pred_pat <- hmc_pred_patient |>
  filter(cell_group == "M2 Mac") |>
  mutate(cell_group = factor(cell_group))

pcr_cols <- c("RD" = "#4575b4", "pCR" = "#d73027")

ggplot(pred_pat, aes(x = BiopsyPhase, y = proportion_mean, group = PatientID, colour = pCR)) +
  geom_line(linewidth = 0.4, alpha = 0.5) +
  facet_grid(cell_group ~ Arm + pCR, scales = "free_y") +
  scale_colour_manual(values = pcr_cols, name = "Response") +
  labs(
    title = "Between-patient variation in predicted M2 Mac proportion",
    x = NULL,
    y = "Predicted proportion"
  ) +
  theme_pub() +
  theme(legend.position = "bottom", strip.text = element_text(size = 8))

Figure 5

Visualise results in spatial context

Figure 5 shows the two chemoimmunotherapy post-treatment phenotypes from the findings above on the tissue, for 20 representative patients per response group (those closest to each group’s median post-treatment M2 proportion). All post-treatment regions of interest are shown for each patient (pathological complete response on the left, residual disease on the right; magenta = M2 macrophages, yellow = PD-L1-positive antigen-presenting cells, green = other immune cells).

cores <- read_csv(file.path(post_dir, "spatial-cores.csv"), show_col_types = FALSE)
cells <- read_fst(file.path(post_dir, "spatial-cells.fst")) |> as_tibble()

mk_roi <- function(image_id) {
  if (is.na(image_id)) return(patchwork::plot_spacer())
  mk_spatial_tile(filter(cells, ImageID == image_id))
}

mk_row <- function(pid) {
  image_ids <- cores |>
    filter(PatientID == pid) |>
    arrange(roi) |>
    pull(ImageID)
  ids <- c(image_ids, rep(NA_character_, 3L - length(image_ids)))
  patchwork::wrap_plots(
    ggplot() +
      annotate("text", x = 0.5, y = 0.5, label = pid, fontface = "bold", size = 2.8) +
      scale_x_continuous(limits = c(0, 1), expand = c(0, 0)) +
      scale_y_continuous(limits = c(0, 1), expand = c(0, 0)) +
      theme_void(),
    patchwork::wrap_plots(lapply(ids, mk_roi), nrow = 1, ncol = 3),
    ncol = 1,
    heights = c(0.16, 1)
  )
}

mk_block <- function(response, title) {
  one_column <- function(which_side) {
    pids <- cores |>
      filter(pCR == response, side == which_side) |>
      arrange(patient) |>
      distinct(PatientID) |>
      pull(PatientID)
    patchwork::wrap_plots(lapply(pids, mk_row), ncol = 1)
  }
  patchwork::wrap_plots(one_column("left"), one_column("right"), ncol = 2) +
    patchwork::plot_annotation(
      title = title,
      theme = theme(
        plot.title = element_text(face = "bold", size = 12, hjust = 0.5, margin = margin(b = 4))
      )
    )
}

divider <- ggplot() +
  annotate("segment", x = 0.5, xend = 0.5, y = 0, yend = 1, linewidth = 0.7, colour = "grey35") +
  scale_x_continuous(limits = c(0, 1), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, 1), expand = c(0, 0)) +
  theme_void()

patchwork::wrap_plots(
  mk_block("pCR", "pCR (responder)"),
  divider,
  mk_block("RD", "RD (non-responder)"),
  ncol = 3,
  widths = c(1, 0.03, 1)
) +
  patchwork::plot_annotation(
    title = "Opposite-direction immune remodelling: C&I post-treatment",
    theme = theme(plot.title = element_text(face = "bold", size = 13))
  )

Figure 6

Post-treatment M2 levels alone look similar between the two groups (21.9% in pCR vs 18.2% in RD among the displayed patients; 22.8% vs 20.9% across the cohort). That is why the visit-only contrast is a poor summary, and why the difference-in-differences contrast was used.

---

sessionInfo()

R version 4.6.1 (2026-06-24)
Platform: aarch64-apple-darwin23
Running under: macOS Tahoe 26.6.2

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8

time zone: Australia/Sydney
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] fstcore_0.10.0  tidyr_1.3.2     dplyr_1.2.1     readr_2.2.0    
[5] patchwork_1.3.2 here_1.0.2      fst_0.9.8       ggplot2_4.0.3  

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1            farver_2.1.2               
 [3] S7_0.2.2                    fastmap_1.2.0              
 [5] SingleCellExperiment_1.34.0 tensorA_0.36.2.1           
 [7] digest_0.6.39               lifecycle_1.0.5            
 [9] sccomp_2.5.2                processx_3.9.0             
[11] magrittr_2.0.5              posterior_1.7.0            
[13] compiler_4.6.1              sass_0.4.10                
[15] rlang_1.3.0                 tools_4.6.1                
[17] gt_1.4.0                    utf8_1.2.6                 
[19] yaml_2.3.12                 data.table_1.18.4          
[21] knitr_1.51                  S4Arrays_1.12.0            
[23] labeling_0.4.3              htmlwidgets_1.6.4          
[25] bit_4.6.0                   DelayedArray_0.38.2        
[27] xml2_1.6.0                  RColorBrewer_1.1-3         
[29] cmdstanr_0.9.0              abind_1.4-8                
[31] withr_3.0.3                 purrr_1.2.2                
[33] BiocGenerics_0.58.1         grid_4.6.1                 
[35] stats4_4.6.1                scales_1.4.0               
[37] SummarizedExperiment_1.42.0 cli_3.6.6                  
[39] rmarkdown_2.31              crayon_1.5.3               
[41] generics_0.1.4              otel_0.2.0                 
[43] tzdb_0.5.0                  stringr_1.6.0              
[45] parallel_4.6.1              XVector_0.52.0             
[47] matrixStats_1.5.0           vctrs_0.7.3                
[49] Matrix_1.7-5                jsonlite_2.0.0             
[51] callr_3.8.0                 IRanges_2.46.0             
[53] hms_1.1.4                   S4Vectors_0.50.1           
[55] bit64_4.8.2                 ggrepel_0.9.8              
[57] glue_1.8.1                  ggside_0.4.1               
[59] distributional_0.8.1        stringi_1.8.7              
[61] gtable_0.3.6                GenomicRanges_1.64.0       
[63] instantiate_0.2.3           tibble_3.3.1               
[65] pillar_1.11.1               htmltools_0.5.9            
[67] Seqinfo_1.2.0               R6_2.6.1                   
[69] rprojroot_2.1.1             vroom_1.7.1                
[71] evaluate_1.0.5              Biobase_2.72.0             
[73] lattice_0.22-9              backports_1.5.1            
[75] Rcpp_1.1.2                  SparseArray_1.12.2         
[77] checkmate_2.3.4             xfun_0.60                  
[79] MatrixGenerics_1.24.0       fs_2.1.0                   
[81] forcats_1.0.1               pkgconfig_2.0.3            

© 2025 tidyomics. Content is published under Creative Commons CC-BY-4.0 License for the text and BSD 3-Clause License for any code. | R-Bloggers