Spatially-resolved single-cell imaging of melanoma brain metastases identifies localized immune patterns predictive of immune checkpoint blockade response.

Code ↔ Paper

18 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 18 matches

[1] § Methods › Survival Analyses ↔ FinalFigures_public.Rmd, lines 4116–4166 · score 0.81 · dummy encoded, categorical variables, continuous variables, missing clinical, clinical variables, imputed
[2] § Results › Phenotypic and Spatial Features of the TME Are Associated With Survival of Patients With MBM ↔ FinalFigures_public.Rmd, lines 3657–3717 · score 0.80 · Myeloid enriched, CD8Tc border, CD8Tc core, Tumor core, CD4Tc, plasma cell
[3] § Results › Highly Multiplexed IMC Analysis of Human MBM Samples ↔ FinalFigures_public.Rmd, lines 119–227 · score 0.79 · extracranial controlled disease, NRAS mutations, histological subtype, wt, anatomical, wild
[4] § Results › TMEs of Intracranial ICI Responders and Non-Responders Are Phenotypically and Spatially Different ↔ FinalFigures_public.Rmd, lines 2120–2201 · score 0.73 · median PFS, median OS, ICI responder, infiltration pattern, postoperative ICI, dichotomizing
[5] § Results › The TMEs of MBM Are Phenotypically Heterogeneous ↔ FinalFigures_public.Rmd, lines 3031–3095 · score 0.72 · CD4Tc, CD8Tc, Neu, exh, Mf, astrocytes
[6] § Results › TMEs of Intracranial ICI Responders and Non-Responders Are Phenotypically and Spatially Different ↔ FinalFigures_public.Rmd, lines 3031–3095 · score 0.70 · NK cells, ICI pretreated, ICI naive, IDO, Arginase, GZB
[7] § Results › Phenotypic and Spatial Features of the TME Are Associated With Survival of Patients With MBM ↔ FinalFigures_public.Rmd, lines 2431–2454 · score 0.70 · plasma cell enriched, CN clusters, tumor core, perivascular, Treg, myeloid
[8] § Results › Highly Multiplexed IMC Analysis of Human MBM Samples ↔ FinalFigures_public.Rmd, lines 119–227 · score 0.68 · BRAF mutation, NRAS mutations, clinical cohort, age, human, ICI
[9] § Methods › Phenotypic and Spatial Analyses ↔ FinalFigures_public.Rmd, lines 1663–1750 · score 0.67 · imcRtools, Spatial interactions, classic, sum, avoidances, scores
[10] § Results › TMEs of Intracranial ICI Responders and Non-Responders Are Phenotypically and Spatially Different ↔ FinalFigures_public.Rmd, lines 2120–2201 · score 0.66 · median OS, ICI responders, ICI naive, ICI response, longer, postoperatively
[11] § Results › The TMEs of MBM Are Phenotypically Heterogeneous ↔ FinalFigures_public.Rmd, lines 1055–1166 · score 0.65 · intraclass correlation coefficient, linear mixed, variance, identity, inter, ROI
[12] § Results › Phenotypic and Spatial Features of the TME Are Associated With Survival of Patients With MBM ↔ FinalFigures_public.Rmd, lines 3252–3293 · score 0.63 · multivariate penalized Cox, Cox proportional hazards, proportional hazards model, Ridge regression, CN, Likelihood
[13] § Methods › Cell Type Annotation ↔ FinalFigures_public.Rmd, lines 699–806 · score 0.63 · MelanA, broad cell, CD3, cytomapper, myeloid, tumor
[14] § Results › TMEs of Intracranial ICI Responders and Non-Responders Are Phenotypically and Spatially Different ↔ FinalFigures_public.Rmd, lines 4116–4166 · score 0.62 · logistic regression models, clinical variables, ICI response, ROC, elastic, responders
[15] § Results › TMEs of Intracranial ICI Responders and Non-Responders Are Phenotypically and Spatially Different ↔ FinalFigures_public.Rmd, lines 2964–3029 · score 0.62 · spatial co localization, ICI pretreated, ICI responders, naive, tumors, cell
[16] § Results › Phenotypic and Spatial Features of the TME Are Associated With Survival of Patients With MBM ↔ FinalFigures_public.Rmd, lines 2964–3029 · score 0.59 · cold samples, median OS, tumor samples, localized, PFS, survival
[17] § Methods › Single-Cell Segmentation, Data Processing and Quality Control ↔ steinbock/measurement/neighbors.py, lines 139–160 · score 0.56 · Euclidean pixel expansion, dmax, Neighboring, Steinbock, masks
[18] § Results › TMEs of Intracranial ICI Responders and Non-Responders Are Phenotypically and Spatially Different ↔ FinalFigures_public.Rmd, lines 2879–2962 · score 0.56 · logFC, ICI responders, ICI naive, FDR, clinical, cell

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R Markdown · 4,281 lines · 195 KB · no license · 17 matches

---
title: "Spatial single-cell analysis of melanoma brain metastases"
subtitle: "Figure generation code for Voglis et al. 2026, Neuro-Oncology"
output: html_document
date: "2024-01-06"
---
<!--
================================================================================
Analysis and figure-generation code accompanying:
Voglis et al. (2026), Neuro-Oncology.
This document reproduces the main and supplementary figures from imaging mass
cytometry (IMC) data of melanoma brain metastases. The starting point is a
pre-processed SingleCellExperiment (SCE) object together with the corresponding
segmentation masks and images. Each top-level section corresponds to a figure
or analysis block; the workflow proceeds from cohort overview (Figure 1) through
the spatial / survival analyses (Figures 2-4) and the revision analyses.
The code requires a pre-built SCE object that is not redistributed here in raw
form; a publication-ready, de-identified version is exported at the end of this
document for deposition on Zenodo.
================================================================================
-->
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
# --- Single-cell / IMC infrastructure ---------------------------------------
library(dittoSeq)
library(CATALYST)
library(scran)
library(scater)
library(bluster)
library(Rphenoannoy)
library(cytomapper)
library(imcRtools)
library(BiocParallel)
# --- Plotting / layout ------------------------------------------------------
library(patchwork)
library(cowplot)
library(viridis)
library(ComplexHeatmap)
library(scales)
library(RColorBrewer)
library(ggrepel)
library(circlize)
library(ggdendro)
library(ggplotify)
library(ggnewscale)
# --- Data wrangling ---------------------------------------------------------
library(tidyverse)
library(data.table)
# --- Statistics / survival --------------------------------------------------
library(survminer)
library(survival)
library(rstatix)
library(edgeR)
# --- Spatial statistics -----------------------------------------------------
library(spatial)
library(spatstat)
library(sf)
library(concaveman)
```
# Load data
```{r}
# Load project-specific helper functions (e.g. run.coxph, plot.forrest,
# plot.volcano), which are sourced rather than redefined inline.
source("0_helperFunctions.R")
# Set input/output paths
datadir <- "/mnt/projects_output_data/MBM"
output_path <- file.path(datadir, "data")
# Pre-processed SingleCellExperiment object (after spatial preparation) and the
# corresponding segmentation masks.
sce <- readRDS(file.path(output_path, "sce_after_spatialprepare.rds"))
masks <- readRDS(file.path(output_path, "masks.rds"))
```
# Select cell types
```{r}
# Promote selected fine-grained myeloid and T-cell labels into the unified
# `cl_all_spec` annotation used throughout the figures. This keeps the broad
# `cl_all` labels intact while exposing specific functional states
# (e.g. CD38+/IDO+ macrophages and GZB+ activated CD8 T cells).
sce$cl_all_spec[sce$cl_myeloid_combined %in% c("CD38+ Mf", "IDO+ Mf")] <-
sce$cl_myeloid_combined[sce$cl_myeloid_combined %in% c("CD38+ Mf", "IDO+ Mf")]
sce$cl_all_spec[sce$cl_tcells_combined %in% c("GZB+ act.CD8Tc")] <-
sce$cl_tcells_combined[sce$cl_tcells_combined %in% c("GZB+ act.CD8Tc")]
unique(sce$cl_all_spec)
# Build a tumor-only annotation that labels tumor subclones as "Tumor <n>" and
# any remaining tumor cells simply as "Tumor".
sce$cl_all_tumor <- sce$cl_tumor_a
sce$cl_all_tumor[sce$cl_all_tumor != 0 & sce$cl_all_tumor != ""] <-
paste0("Tumor ", sce$cl_tumor_a[sce$cl_all_tumor != 0 & sce$cl_all_tumor != ""])
sce$cl_all_tumor[sce$cl_all_tumor == ""] <- "Tumor"
unique(sce$cl_all_tumor)
```
# Define proliferative cells
```{r}
# Dichotomize cells into high vs. low Ki-67 expression based on the cohort-wide
# median of the arcsinh-transformed Ki-67 signal.
median_ki67 <- median(assays(sce)[["exprs"]]["Ki-67", ])
sce$ki67 <- ifelse(assays(sce)[["exprs"]]["Ki-67", ] > median_ki67,
"high_ki67", "low_ki67")
```
# --------------------------------
# FIGURE 1
## Metadata heatmap
```{r, fig.height = 18}
library(SCpubr)
# Collect the clinical variables to be displayed in the metadata heatmap.
clindata <- metadata(sce)$clinical %>%
select(biosample_id, sex, age, localization, primarius_type, braf, nras,
extracranial_control_dos, rt_bm_preop, rt_bm_postop, ici_preop, ici_postop)
# IMPORTANT: restrict to the final analysis cohort, i.e. samples that survived QC
# and are still present in the SCE object.
clindata <- colData(sce) %>%
as_tibble() %>%
distinct(biosample_id, cond.1, cond.2, cond.3) %>%
left_join(clindata)
# Drop identifier / grouping columns not shown in the heatmap body.
clindata.select <- clindata %>% select(-biosample_id, -age, -cond.1, -cond.2, -cond.3)
# Recode BRAF/NRAS mutation status to human-readable labels.
clindata.select$braf <- as.character(clindata.select$braf)
clindata.select$braf[clindata.select$braf == "mut"] <- "mutated"
clindata.select$braf[clindata.select$braf == "wt"] <- "wild type"
clindata.select$braf[is.na(clindata.select$braf)] <- "not determined"
clindata.select$nras <- as.character(clindata.select$nras)
clindata.select$nras[clindata.select$nras == "mut"] <- "mutated"
clindata.select$nras[clindata.select$nras == "wt"] <- "wild type"
clindata.select$nras[is.na(clindata.select$nras)] <- "not determined"
# Rename columns to publication labels.
clindata_renamed <- clindata.select %>%
rename("Sex" = sex,
"Anatomical localization" = localization,
"Histological subtype" = primarius_type,
"BRAF mutation" = braf,
"NRAS mutation" = nras,
"Extracranial controlled disease" = extracranial_control_dos,
"Radiotherapy preop" = rt_bm_preop,
"Radiotherapy postop" = rt_bm_postop,
"ICI preop" = ici_preop,
"ICI postop" = ici_postop)
# Define colour schemes for the categorical annotations.
col_ann_yesno <- c("yes" = "#E69F00", "no" = "#56B4E9")
col_ann_mut <- c("mutated" = "#D55E00", "wild type" = "#56B4E9", "not determined" = "#D3D3D3")
col_ann_metadata <- list(
"Anatomical localization" = setNames(dittoColors()[seq_along(unique(clindata$localization))],
unique(clindata$localization)),
"Histological subtype" = setNames(dittoColors()[seq_along(unique(clindata$primarius_type))],
unique(clindata$primarius_type)),
"BRAF mutation" = col_ann_mut,
"NRAS mutation" = col_ann_mut,
"Extracranial controlled disease" = col_ann_yesno,
"Radiotherapy preop" = col_ann_yesno,
"Radiotherapy postop" = col_ann_yesno,
"ICI preop" = col_ann_yesno,
"ICI postop" = col_ann_yesno
)
# Build the base metadata heatmap.
(p.clin <- do_MetadataPlot(from_df = T,
df = clindata_renamed,
flip = F,
legend.nrow = 7,
font.size = 22,
heatmap.gap = 1,
cluster = F,
colors.use = col_ann_metadata,
grid.color = "black",
border.color = "black",
legend.title.face = "plain",
legend.position = "bottom"))
# Merge paired tracks (preop/postop, BRAF/NRAS) under shared legends and hide the
# redundant per-track legends.
p.clin[[6]] <- p.clin[[6]] + scale_fill_manual(name = "Extracranial\ncontrolled disease",
breaks = c("yes", "no", NA),
values = col_ann_yesno,
na.value = "#D3D3D3")
p.clin[[7]] <- p.clin[[7]] + scale_fill_manual(name = "Radiotherapy preop/postop",
breaks = c("yes", "no", NA),
values = col_ann_yesno,
na.value = "#D3D3D3")
p.clin[[8]] <- p.clin[[8]] + theme(legend.position = "none")
p.clin[[4]] <- p.clin[[4]] + scale_fill_manual(name = "BRAF/NRAS mutation",
values = col_ann_mut,
breaks = c("wild type", "mutated", "not determined"),
na.value = "#D3D3D3")
p.clin[[5]] <- p.clin[[5]] + theme(legend.position = "none")
p.clin[[9]] <- p.clin[[9]] + scale_fill_manual(name = "ICI preop/postop",
breaks = c("yes", "no", NA),
values = col_ann_yesno,
na.value = "#D3D3D3")
p.clin[[10]] <- p.clin[[10]] + theme(legend.position = "none")
ggsave(plot = p.clin, "figures/fig1_hm_clinical.pdf", width = 25, heigh = 15, dpi = 200)
# Tabular summary of the clinical cohort.
library(gtsummary)
clindata %>% select(-biosample_id) %>% tbl_summary()
```
## Revision: new metadata plot
```{r, fig.width = 18, fig.height = 8}
library(ggh4x)
# Ordered row labels for the redesigned metadata grid (bottom-to-top).
metadata.labels <- c("ICI postop",
"ICI preop",
"RT postop",
"RT preop",
"Extracranial controlled disease",
"NRAS mut",
"BRAF mut",
"Histological subtype",
"Anatomical localization",
"Sex")
# Group sizes (number of samples per ICI-pretreatment status), used to scale the
# facet panel widths so each tile is square.
group_sizes <- clindata.select %>%
count(ici_preop) %>%
arrange(ici_preop)
tile_w <- 0.95 # cm per tile
# Build the grid one variable at a time. Each `geom_tile` adds one metadata row;
# `new_scale_fill()` allows an independent fill scale per row. Paired tracks
# (preop/postop, BRAF/NRAS) share a legend, with the second track's legend hidden.
p <- clindata.select %>%
arrange(ici_preop) %>%
group_by(ici_preop) %>%
mutate(local_id = row_number()) %>%
ungroup() %>%
ggplot(aes(x = local_id)) +
geom_tile(aes(y = 1, fill = ici_postop), color = "black", width = .9, height = 0.8) +
scale_fill_manual(name = "ICI preop/postop",
breaks = c("yes", "no", NA),
values = col_ann_yesno,
na.value = "#D3D3D3") +
new_scale_fill() +
geom_tile(aes(y = 2, fill = ici_preop), color = "black", width = .9, height = 0.8) +
scale_fill_manual(name = "ICI preop/postop",
breaks = c("yes", "no", NA),
values = col_ann_yesno,
na.value = "#D3D3D3",
guide = "none") +
new_scale_fill() +
geom_tile(aes(y = 3, fill = rt_bm_postop), color = "black", width = .9, height = 0.8) +
scale_fill_manual(name = "RT BM preop/postop",
breaks = c("yes", "no", NA),
values = col_ann_yesno,
na.value = "#D3D3D3") +
new_scale_fill() +
geom_tile(aes(y = 4, fill = rt_bm_preop), color = "black", width = .9, height = 0.8) +
scale_fill_manual(name = "RT BM preop/postop",
breaks = c("yes", "no", NA),
values = col_ann_yesno,
na.value = "#D3D3D3",
guide = "none") +
new_scale_fill() +
geom_tile(aes(y = 5, fill = extracranial_control_dos), color = "black", width = .9, height = 0.8) +
scale_fill_manual(name = "Extracranial\ncontrolled disease",
breaks = c("yes", "no", NA),
values = col_ann_yesno,
na.value = "#D3D3D3") +
new_scale_fill() +
geom_tile(aes(y = 6, fill = nras), color = "black", width = .9, height = 0.8) +
scale_fill_manual(name = "BRAF/NRAS mutation",
values = col_ann_mut,
breaks = c("wild type", "mutated", "not determined"),
na.value = "#D3D3D3") +
new_scale_fill() +
geom_tile(aes(y = 7, fill = braf), color = "black", width = .9, height = 0.8) +
scale_fill_manual(name = "BRAF/NRAS mutation",
values = col_ann_mut,
breaks = c("wild type", "mutated", "not determined"),
na.value = "#D3D3D3",
guide = "none") +
new_scale_fill() +
geom_tile(aes(y = 8, fill = primarius_type), color = "black", width = .9, height = 0.8) +
scale_fill_manual(name = "Histological subtype",
values = setNames(dittoColors()[seq_along(unique(clindata$primarius_type))],
unique(clindata$primarius_type)),
na.value = "#D3D3D3") +
new_scale_fill() +
geom_tile(aes(y = 9, fill = localization), color = "black", width = .9, height = 0.8) +
scale_fill_manual(name = "Anatomical localization",
values = setNames(dittoColors()[seq_along(unique(clindata$localization))],
unique(clindata$localization)),
na.value = "#D3D3D3") +
new_scale_fill() +
geom_tile(aes(y = 10, fill = sex), color = "black", width = .9, height = 0.8) +
scale_fill_manual(name = "Sex",
values = c("F" = "#477c9e", "M" = "#9e6845")) +
scale_y_continuous(
breaks = 1:length(metadata.labels),
labels = str_wrap(metadata.labels, width = 15)
) +
scale_x_continuous(
breaks = function(x) seq(ceiling(x[1]), floor(x[2]), by = 1),
expand = expansion(add = 0.5) # half-tile padding on each side
) +
facet_grid(. ~ ici_preop, scales = "free_x",
labeller = as_labeller(c("yes" = "ICI pretreated", "no" = "ICI naive"))) +
theme_minimal() +
theme(
panel.grid = element_blank(),
axis.title = element_blank(),
axis.text.x = element_blank(),
axis.ticks = element_line(color = "black"),
axis.ticks.length = unit(0.1, "cm"),
axis.text.y = element_text(lineheight = 0.7),
legend.box = "horizontal",
legend.box.just = "top",
strip.text = element_text(face = "bold"),
text = element_text(size = 16.5)
)
# Force square tiles by fixing panel sizes to (n tiles x tile width).
p.metadata.final <- p + force_panelsizes(
cols = unit(group_sizes$n * tile_w, "cm"),
rows = unit(length(metadata.labels) * tile_w, "cm")
)
# Export the legend and the (legend-free) grid separately for figure assembly.
legend <- get_legend(p.metadata.final)
plot_grid(legend)
p.metadata.final + theme(legend.position = "none")
ggsave(plot = p.metadata.final + theme(legend.position = "none"),
"figures/Revision_fig1_hm_clinical.pdf", width = 18, heigh = 8, dpi = 300)
ggsave(plot = plot_grid(legend),
"figures/Revision_fig1_hm_clinical_legend.pdf", width = 18, heigh = 8, dpi = 300)
```
## Test: Radial (polar) layout
```{r, fig.width = 12, fig.height = 15}
# Alternative circular rendering of the metadata grid, wrapping the same tiles
# around a polar coordinate system. The empty central hole is created by giving
# the y-axis a negative lower limit; ring tick marks and labels are drawn
# manually in that central region.
ring_ticks <- data.frame(
xstart = 0, # tick start
xend = 0.5, # tick end (gap centre)
y = seq_along(metadata.labels)
)
ring_labels <- data.frame(
local_id = -0.1, # label anchor, just inside the central hole
y = seq_along(metadata.labels),
label = str_wrap(metadata.labels, width = 15)
)
clindata.select %>%
arrange(ici_preop) %>%
group_by(ici_preop) %>%
mutate(local_id = row_number()) %>%
ungroup() %>%
ggplot(aes(x = local_id)) +
geom_tile(aes(y = 1, fill = ici_postop), color = "black", width = .9, height = 0.7) +
scale_fill_manual(name = "ICI preop/postop",
breaks = c("yes", "no", NA),
values = col_ann_yesno,
na.value = "#D3D3D3") +
new_scale_fill() +
geom_tile(aes(y = 2, fill = ici_preop), color = "black", width = .9, height = 0.7) +
scale_fill_manual(name = "ICI preop/postop",
breaks = c("yes", "no", NA),
values = col_ann_yesno,
na.value = "#D3D3D3",
guide = "none") +
new_scale_fill() +
geom_tile(aes(y = 3, fill = rt_bm_postop), color = "black", width = .9, height = 0.7) +
scale_fill_manual(name = "RT BM preop/postop",
breaks = c("yes", "no", NA),
values = col_ann_yesno,
na.value = "#D3D3D3") +
new_scale_fill() +
geom_tile(aes(y = 4, fill = rt_bm_preop), color = "black", width = .9, height = 0.7) +
scale_fill_manual(name = "RT BM preop/postop",
breaks = c("yes", "no", NA),
values = col_ann_yesno,
na.value = "#D3D3D3",
guide = "none") +
new_scale_fill() +
geom_tile(aes(y = 5, fill = extracranial_control_dos), color = "black", width = .9, height = 0.7) +
scale_fill_manual(name = "Extracranial\ncontrolled disease",
breaks = c("yes", "no", NA),
values = col_ann_yesno,
na.value = "#D3D3D3") +
new_scale_fill() +
geom_tile(aes(y = 6, fill = nras), color = "black", width = .9, height = 0.7) +
scale_fill_manual(name = "BRAF/NRAS mutation",
values = col_ann_mut,
breaks = c("wild type", "mutated", "not determined"),
na.value = "#D3D3D3") +
new_scale_fill() +
geom_tile(aes(y = 7, fill = braf), color = "black", width = .9, height = 0.7) +
scale_fill_manual(name = "BRAF/NRAS mutation",
values = col_ann_mut,
breaks = c("wild type", "mutated", "not determined"),
na.value = "#D3D3D3",
guide = "none") +
new_scale_fill() +
geom_tile(aes(y = 8, fill = primarius_type), color = "black", width = .9, height = 0.7) +
scale_fill_manual(name = "Histological subtype",
values = setNames(dittoColors()[seq_along(unique(clindata$primarius_type))],
unique(clindata$primarius_type)),
na.value = "#D3D3D3") +
new_scale_fill() +
geom_tile(aes(y = 9, fill = localization), color = "black", width = .9, height = 0.7) +
scale_fill_manual(name = "Anatomical localization",
values = setNames(dittoColors()[seq_along(unique(clindata$localization))],
unique(clindata$localization)),
na.value = "#D3D3D3") +
new_scale_fill() +
geom_tile(aes(y = 10, fill = sex), color = "black", width = .9, height = 0.7) +
scale_fill_manual(name = "Sex",
values = c("F" = "#477c9e", "M" = "#9e6845")) +
# Negative lower y-limit creates the empty central hole of the ring.
scale_y_continuous(
limits = c(-9, length(metadata.labels) + 0.5),
breaks = 1:length(metadata.labels),
labels = NULL
) +
scale_x_continuous(expand = expansion(add = 0.5)) +
# Manual ring tick marks and labels in the central region.
geom_segment(
data = ring_ticks,
aes(x = xstart, xend = xend, y = y, yend = y),
color = "grey30", linewidth = 0.4, inherit.aes = FALSE
) +
geom_text(
data = ring_labels,
aes(x = local_id, y = y, label = label),
hjust = 1, size = 3.2, lineheight = 0.8,
fontface = "bold", color = "grey30", inherit.aes = FALSE
) +
# facet_wrap works better with coord_polar than facet_grid.
facet_wrap(~ ici_preop,
ncol = 1,
labeller = as_labeller(c("yes" = "ICI pretreated", "no" = "ICI naive"))) +
coord_polar(theta = "x") +
theme_minimal() +
theme(
panel.grid = element_blank(),
axis.title = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.ticks.length = unit(0.1, "cm"),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
legend.box = "horizontal",
legend.box.just = "top",
strip.text = element_text(face = "bold"),
legend.position = "none"
)
```
# --------------------------------
# SUPPL. FIGURE 1
## Import clinical data
Here we load additional clinical data used for the swimmer and survival plots
(follow-up responses and per-regimen ICI treatment characteristics) and restrict
all tables to the final analysis cohort.
```{r}
clin.data <- readRDS(file.path(output_path, "clinical_data.rds"))
load(file.path(output_path, "clinical_data_add.RData"))
# The .RData file provides:
# df = clinical data (as imported above for the SCE object)
# resp = follow-up response data for each available timepoint
# tr_ct.ici = ICI treatment characteristics for each individual regimen
# IMPORTANT: restrict to the final analysis cohort (samples present in the SCE
# object after QC).
df <- colData(sce) %>% as_tibble() %>% distinct(biosample_id) %>% left_join(df)
resp <- colData(sce) %>% as_tibble() %>% distinct(biosample_id) %>% left_join(resp)
tr_ct.ici <- colData(sce) %>% as_tibble() %>% distinct(biosample_id) %>% left_join(tr_ct.ici)
```
## Swimmer plots
### Naive vs. pretreated
```{r}
# Treatment timeline (orange bars) with follow-up response timepoints (coloured
# points) and dates of surgery of other lesions (blue asterisks), for samples
# that received ICI preoperatively.
p1 <- ggplot() +
geom_linerange(data = tr_ct.ici %>% filter(id %in% df$id[df$ici_preop == "yes"]),
aes(ymin = day_start, ymax = day_stop, x = as.factor(id)),
stat = "identity", position = position_dodge(width = 1),
size = 7, colour = "orange") +
geom_point(data = resp %>% filter(id %in% df$id[df$ici_preop == "yes"]),
aes(x = as.factor(id), y = fu_day, color = fu_response), size = 4) +
geom_point(data = df %>% filter(ici_preop == "yes") %>%
select(patnr, id, dos) %>% rowwise() %>%
mutate(other_dos = paste(as.numeric(difftime(df$dos[df$patnr == patnr], dos,
units = "days")), collapse = ";")) %>%
separate_rows(other_dos, sep = ";", convert = T),
aes(x = as.factor(id), y = other_dos), shape = 8, color = "blue", size = 3) +
theme(axis.text.y = element_blank(),
axis.title.y = element_blank()) +
geom_hline(aes(yintercept = 0)) +
coord_flip()
p1.legend <- cowplot::get_legend(p1)
p1 <- p1 + theme(legend.position = "none")
# Patient ID tiles.
p.pat <- ggplot() +
geom_label(data = df %>% filter(ici_preop == "yes") %>% select(id, pid),
aes(x = factor(id), y = "PatientID", label = pid), label.size = NA) +
coord_flip() +
theme_void() + theme(axis.text.x = element_text(angle = 90))
# Sample ID tiles.
p.id <- ggplot() +
geom_label(data = df %>% filter(ici_preop == "yes") %>% select(id, pid),
aes(x = factor(id), y = "SampleID", label = id), label.size = NA) +
coord_flip() +
theme_void() + theme(axis.text.x = element_text(angle = 90))
# Tile annotation: progression after preoperative ICI.
p2 <- ggplot() +
geom_tile(data = df %>% filter(ici_preop == "yes") %>% select(id, ici_preop, pd_after_ici_preop),
aes(x = factor(id), y = 3000, fill = factor(pd_after_ici_preop))) +
scale_fill_manual(name = "pd_after_ici", values = c("no" = "white", "yes" = "blue"),
breaks = c("yes")) +
coord_flip() +
theme_void()
p2.legend <- cowplot::get_legend(p2)
p2 <- p2 + theme(legend.position = "none")
# Tile annotation: progression while on preoperative ICI.
p3 <- ggplot() +
geom_tile(data = df %>% filter(ici_preop == "yes") %>% select(id, ici_preop, pd_while_ici_preop),
aes(x = factor(id), y = 3000, fill = factor(pd_while_ici_preop))) +
scale_fill_manual(name = "pd_while_ici", values = c("no" = "white", "yes" = "green"),
breaks = c("yes")) +
coord_flip() +
theme_void()
p3.legend <- cowplot::get_legend(p3)
p3 <- p3 + theme(legend.position = "none")
# Assemble swimmer plot with side tiles and stacked legends.
ggarrange(p3, p2, p.pat, p.id, p1,
ggarrange(p1.legend, p2.legend, p3.legend, nrow = 3, align = "v"),
ncol = 6, widths = c(1, 1, 3, 3, 100, 15), align = "h")
```
### Naive with response
```{r}
# Note: IDs 6 and 7 can also be treated as effectively ICI-naive (brain
# metastasis diagnosed directly preoperatively, with ICI started thereafter).
# Select ICI-naive cases that received postoperative ICI, plus IDs 6/7.
ids <- df %>% filter((ici_preop == "no" & ici_postop == "yes") | id %in% c(6, 7)) %>% pull(id)
# Treatment timeline with follow-up responses, other surgery dates, and last
# follow-up / death status markers.
p1 <- ggplot() +
geom_linerange(data = tr_ct.ici %>% filter(id %in% ids),
aes(ymin = day_start, ymax = day_stop, x = as.factor(id)),
stat = "identity", position = position_dodge(width = 1),
size = 7, colour = "orange") +
geom_point(data = resp %>% filter(id %in% ids),
aes(x = as.factor(id), y = fu_day, color = fu_response), size = 4) +
geom_point(data = df %>% filter(id %in% ids) %>%
select(patnr, id, dos) %>% rowwise() %>%
mutate(other_dos = paste(as.numeric(difftime(df$dos[df$patnr == patnr], dos,
units = "days")), collapse = ";")) %>%
separate_rows(other_dos, sep = ";", convert = T),
aes(x = as.factor(id), y = other_dos), shape = 8, color = "blue", size = 3) +
geom_point(data = df %>% filter(ici_preop == "no") %>% select(id, days_dos, status),
aes(x = factor(id), y = days_dos, shape = factor(status))) +
theme(axis.text.y = element_blank(),
axis.title.y = element_blank()) +
geom_hline(aes(yintercept = 0)) +
coord_flip()
scale_y_continuous(limits = c(-15, 2000), expand = c(0, 0))
p1
p1.legend <- cowplot::get_legend(p1)
p1 <- p1 + theme(legend.position = "none")
# Tile annotation: postoperative ICI responder status.
p2 <- ggplot() +
geom_tile(data = df %>% filter(id %in% ids) %>% select(id, ici_postop_response),
aes(x = factor(id), y = 1, fill = factor(ici_postop_response))) +
scale_fill_discrete(name = "ICI responder") +
coord_flip() +
theme_void()
p2.legend <- cowplot::get_legend(p2)
p2 <- p2 + theme(legend.position = "none")
# Patient ID tiles.
p.pat <- ggplot() +
geom_label(data = df %>% filter(id %in% ids) %>% select(id, pid),
aes(x = factor(id), y = "PatientID", label = pid), label.size = NA, size = 5) +
coord_flip() +
theme_void() + theme(axis.text.x = element_text(angle = 90))
# Sample ID tiles.
p.id <- ggplot() +
geom_label(data = df %>% filter(id %in% ids) %>% select(id, pid),
aes(x = factor(id), y = "SampleID", label = id), label.size = NA, size = 5) +
coord_flip() +
theme_void() + theme(axis.text.x = element_text(angle = 90))
ggarrange(p2, p.pat, p.id, p1,
ggarrange(p1.legend, p2.legend, nrow = 2, align = "v"),
ncol = 5, widths = c(1, 3, 3, 100, 15), align = "h")
```
## Survival plots
### Overall survival
```{r survivalplot-os}
library(survival)
library(survminer)
# Clinical data from the SCE object.
dff <- metadata(sce)$clinical
# Condition 1: Kaplan-Meier overall survival.
surv.df <- dff %>% filter(!is.na(cond.1)) %>% select(cond.1, days_dos, status)
fit <- survfit(Surv(days_dos, status) ~ cond.1, data = surv.df)
median <- surv_median(fit)
ggsurvplot(fit, risk.table = TRUE, conf.int = F, pval = TRUE,
xscale = "d_m", surv.median.line = "hv",
break.time.by = 365.25, xlab = "Time in months",
risk.table.title = "Patients at risk",
fontsize = 4,
ggtheme = theme_bw())
# Condition 2: Kaplan-Meier overall survival.
surv.df <- dff %>% filter(!is.na(cond.2)) %>% select(cond.2, days_dos, status)
fit <- survfit(Surv(days_dos, status) ~ cond.2, data = surv.df)
median <- surv_median(fit)
ggsurvplot(fit, risk.table = TRUE, conf.int = F, pval = TRUE,
xscale = "d_m", surv.median.line = "hv",
break.time.by = 365.25, xlab = "Time in months",
risk.table.title = "Patients at risk",
fontsize = 4,
ggtheme = theme_bw())
```
### Progression-free survival
```{r survivalplot-pfs}
dff <- metadata(sce)$clinical
# Condition 1: Kaplan-Meier progression-free survival.
surv.df <- dff %>% filter(!is.na(cond.1)) %>% select(cond.1, days_dos_pfs, status_pfs)
fit <- survfit(Surv(days_dos_pfs, status_pfs) ~ cond.1, data = surv.df)
median <- surv_median(fit)
ggsurvplot(fit, risk.table = TRUE, conf.int = F, pval = TRUE,
xscale = "d_m", surv.median.line = "hv",
break.time.by = 365.25, xlab = "Time in months",
risk.table.title = "Patients at risk",
fontsize = 4,
ggtheme = theme_bw())
# Condition 2: Kaplan-Meier progression-free survival.
surv.df <- dff %>% filter(!is.na(cond.2)) %>% select(cond.2, days_dos_pfs, status_pfs)
fit <- survfit(Surv(days_dos_pfs, status_pfs) ~ cond.2, data = surv.df)
median <- surv_median(fit)
ggsurvplot(fit, risk.table = TRUE, conf.int = F, pval = TRUE,
xscale = "d_m", surv.median.line = "hv",
break.time.by = 365.25, xlab = "Time in months",
risk.table.title = "Patients at risk",
fontsize = 4,
ggtheme = theme_bw())
```
# --------------------------------
# FIGURE 2
## Plot images
```{r}
# Load a subset of the multichannel images for visualization.
images <- readRDS(file.path(output_path, "images_subset.rds"))
# Recurring plotting parameters (scale bar, image title, legend styling).
scalebar_param <- list(length = 100, label = "", colour = "white",
position = "bottomright", cex = 1, margin = c(5, 5))
imagetitle_param <- list(position = "topright", colour = "white",
margin = c(5, 5), font = 2, cex = 1)
legend_param <- list(colour_by.title.cex = 0.6, colour_by.labels.cex = 0.6,
colour_by.legend.cex = 0.6, outline_by.title.cex = 0.6,
outline_by.labels.cex = 0.6, outline_by.legend.cex = 0.6,
margin = 0)
# Select representative ROIs and the matching masks/cells.
cur_masks <- masks[names(masks) %in% c("NB16-121_002", "NB19-309_011",
"NB17-357_006", "NB18-818_003")]
tmp_sce <- sce[, sce$image_id %in% names(cur_masks)]
# Collapse T-cell subsets into a single "T cells" class for this overview.
tmp_sce$cl_all[tmp_sce$cl_all %in% c("CD8Tc", "CD4Tc", "DPTc", "Treg")] <- "T cells"
# Normalize images per-image for consistent display.
cur_images <- images[names(cur_masks)]
cur_images <- cytomapper::normalize(cur_images, separateImages = T)
# Outline segmented cells coloured by broad cell type, over selected markers.
plot.cells <- plotCells(
object = tmp_sce,
mask = cur_masks,
img_id = "image_id",
cell_id = "ObjectNumber",
colour_by = c("MelanA", "CD3", "GFAP", "CD20", "CD68", "CD31"),
exprs_values = "exprs",
outline_by = "cl_all",
thick = T,
colour = list(MelanA = c("black", "#1f78b4"),
CD3 = c("black", "#33a02c"),
GFAP = c("black", "yellow"),
CD20 = c("black", "#a6cee3"),
CD68 = c("black", "#6a3d9a"),
CD31 = c("black", "#e31a1c"),
cl_all = c("Tumor" = "#1f78b4",
"T cells" = "#33a02c",
"Vascular" = "#e31a1c",
"Myeloid" = "#6a3d9a",
"Neutrophils" = "#ff7f00",
"Astrocytes" = "yellow",
"B cells" = "#a6cee3",
"Plasma" = "magenta",
"NK cells" = "#b15928",
"unassigned" = "white",
"BnT" = "#00bfc4")),
legend = legend_param,
scale_bar = scalebar_param,
display = "single",
image_title = NULL,
return_plot = T
)
# Pseudo-coloured pixel-level marker images (raw signal, per-channel bcg).
plot.images <- plotPixels(
cur_images,
img_id = "image_id",
cell_id = "ObjectNumber",
colour_by = c("MelanA", "CD3", "GFAP", "CD20", "CD68", "CD31"),
bcg = list(MelanA = c(0, 8, 1),
CD3 = c(0, 4, 1),
GFAP = c(0, 4, 1),
CD20 = c(0, 7, 1),
CD68 = c(0, 4, 1),
CD31 = c(0, 4, 1)),
thick = F,
colour = list(MelanA = c("black", "#1f78b4"),
CD3 = c("black", "#33a02c"),
GFAP = c("black", "yellow"),
CD20 = c("black", "#a6cee3"),
CD68 = c("black", "#6a3d9a"),
CD31 = c("black", "#e31a1c")),
return_plot = T,
legend = legend_param,
scale_bar = scalebar_param,
image_title = NULL,
margin = 2,
display = "single"
)
# Free large image objects from memory.
rm(cur_images, images)
# Assemble: top row = cell outlines, bottom row = pixel images, per ROI + legend.
p.img <-
(
ggdraw(plot.cells$plot$`NB16-121_002`, clip = "on") + theme(plot.margin = margin(2, 5, 0, 2)) |
ggdraw(plot.cells$plot$`NB17-357_006`, clip = "on") + theme(plot.margin = margin(2, 5, 0, 0)) |
ggdraw(plot.cells$plot$`NB18-818_003`, clip = "on") + theme(plot.margin = margin(2, 5, 0, 0)) |
ggdraw(plot.cells$plot$legend, clip = "on")
) /
(
ggdraw(plot.images$plot$`NB16-121_002`, clip = "on") + theme(plot.margin = margin(2, 5, 0, 0)) |
ggdraw(plot.images$plot$`NB17-357_006`, clip = "on") + theme(plot.margin = margin(2, 5, 0, 0)) |
ggdraw(plot.images$plot$`NB18-818_003`, clip = "on") + theme(plot.margin = margin(2, 5, 0, 0)) |
ggdraw(plot.images$plot$legend, clip = "on")
)
ggsave(plot = p.img, "figures/fig2_imgs.pdf", width = 12, heigh = 7, dpi = 200)
```
## Heatmap cell types
```{r}
# Drop nuclear markers; keep markers used for any of the clustering steps.
tmp_sce <- sce[rowData(sce)$marker_groups != "nuclei", ]
selected_markers <- rownames(tmp_sce)[
rowSums(as.data.frame(rowData(tmp_sce)[c("cl_celltype", "cl_immune",
"cl_tcells", "cl_myeloid")])) > 0]
# Colour scheme for the cell-type annotation.
col_ann_hm <- c("Tumor" = "#1f78b4",
"B cells" = "#a6cee3",
"Plasma" = "#fdbf6f",
"Vascular" = "#e31a1c",
"CD8Tc" = "#33a02c",
"DPTc" = "#fb9a99",
"exh.CD8Tc" = "#009E73",
"Myeloid" = "#6a3d9a",
"Treg" = "#DCE319FF",
"NK cells" = "#b15928",
"CD4Tc" = "#b2df8a",
"GZB+ act.CD8Tc" = "#00F6B3",
"BnT" = "#00bfc4",
"unassigned" = "#cab2d6",
"IDO+ Mf" = "#B14380",
"CD38+ Mf" = "#4D4D4D",
"Neutrophils" = "#ff7f00",
"Arginase+ Neu" = "#E69F00",
"Astrocytes" = "#ffff99")
# Compute mean marker expression per cluster ID and derive scaled assays.
mean_tmp.sce <- aggregateAcrossCells(tmp_sce, ids = tmp_sce$cl_all_spec, statistics = "mean")
assay(mean_tmp.sce, "exprs") <- asinh(counts(mean_tmp.sce) / 1)
assay(mean_tmp.sce, "scaled") <- t(scale(t(assay(mean_tmp.sce, "exprs")))) # z-score per marker
assay(mean_tmp.sce, "scaled_ztm") <- t(apply(assay(mean_tmp.sce, "exprs"), 1, scales::rescale)) # zero-to-max
# Proportion of Ki-67 high cells per cluster.
prop.ki67 <- colData(tmp_sce) %>% as_tibble() %>%
group_by(cl_all_spec, ki67) %>% summarise(count = n()) %>%
mutate(prop.ki67 = prop.table(count)) %>% select(-count)
cluster_name <- "cl_all_spec"
# Heatmap body: z-scored mean expression (clusters x markers).
hm_body <- t(assay(mean_tmp.sce, "scaled"))
# Row (cluster) annotation.
row_anno <- colData(mean_tmp.sce) %>% as.data.frame %>% select(all_of(cluster_name), ncells)
row_anno[[cluster_name]] <- as.character(row_anno[[cluster_name]])
# Diverging colour scale for z-scores.
col_1 <- colorRamp2(c(-3, 0, 3), c("#4575B4", "white", "#D73027"))
# Right-side annotation: colour bar + cluster names.
ha_right <- HeatmapAnnotation(
cluster_name = anno_simple(row_anno[[cluster_name]], border = TRUE, col = col_ann_hm),
names = anno_text(row_anno[[cluster_name]]),
annotation_label = c(cluster_name, ""),
annotation_name_rot = 90,
which = "row")
# Draw the heatmap, splitting columns by marker group.
hm <- Heatmap(hm_body,
name = "z-score",
col = col_1,
clustering_method_columns = "complete",
clustering_method_rows = "complete",
column_split = rowData(mean_tmp.sce)["marker_groups"],
cluster_column_slices = F,
show_row_names = F,
right_annotation = c(ha_right))
hm <- draw(hm)
p.hm <- grid.grabExpr(draw(hm))
ggsave(plot = p.hm, "figures/fig2_hm_celltypes.pdf", width = 12, heigh = 7, dpi = 200)
# Keep a copy to later extract the cell-type ordering for other plots.
hm.celltypes <- hm
```
## Plot heterogeneity
```{r}
# Per-cell colData used for the composition barplots.
df_plotting <- colData(sce) %>% as_tibble()
# Cluster samples by their broad cell-type composition to order the barplots, and
# build the matching dendrogram.
df_clust <- df_plotting %>%
group_by(biosample_id, cl_all) %>%
summarise(count = n()) %>% mutate(prop = prop.table(count)) %>%
ungroup()
df_clust_wide <- df_clust %>%
pivot_wider(id_cols = biosample_id, names_from = cl_all, values_from = prop)
df_clust_wide <- df_clust_wide %>% replace(is.na(.), 0) # 0 where a cell type is absent
df_clust_m <- as.matrix(df_clust_wide)
rownames(df_clust_m) <- df_clust_wide$biosample_id
# Hierarchical clustering (Ward.D2 on Euclidean distances) -> sample order.
hc <- hclust(dist(df_clust_m[, -1], method = "euclidian"), "ward.D2")
dend <- as.dendrogram(hc)
dend_data <- dendro_data(dend, type = "rectangle")
(p.dend <- ggplot(dend_data$segments) +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend)) +
geom_text(data = dend_data$labels, aes(x, y, label = label), hjust = 1, angle = 90) +
scale_x_continuous(expand = c(0, 0.5)) +
scale_y_continuous(expand = c(0, 0)) +
theme_void() +
theme(plot.margin = margin(0, 0, 0, 0, unit = "mm")))
# Fixed cell-type stacking order for consistent legends/bars.
celltype_order <- c("Tumor", "Vascular", "Astrocytes", "Myeloid", "Neutrophils",
"Plasma", "B cells", "NK cells", "BnT", "DPTc", "CD4Tc",
"Treg", "CD8Tc", "unassigned")
# Per-ROI cell-type counts with sample/condition annotation.
df_barplots <- df_plotting %>%
group_by(image_id, biosample_id, pid, cond.3, cond.2, cond.4, cl_all) %>%
summarise(count = n()) %>%
ungroup()
# Order samples by the clustering above.
df_barplots$biosample_id <- factor(df_barplots$biosample_id, levels = hc$labels[hc$order])
df_barplots$pid <- as.factor(df_barplots$pid)
# Merge dendrogram x-positions so all panels align.
x <- dend_data$labels[c("x", "label")] %>% as_tibble() %>% rename(biosample_id = label)
df_barplots <- df_barplots %>% left_join(x, by = "biosample_id")
# Composition stacked by sample.
(p1 <-
ggplot(data = df_barplots) +
geom_bar(aes(fill = factor(cl_all, levels = celltype_order), x = x, y = count),
stat = "identity", position = "fill", width = 1.1) +
scale_fill_manual(name = "Cell types",
values = metadata(sce)$col_celltypes$celltype_all,
breaks = celltype_order) +
geom_tile(aes(x = x, y = 0, height = 0.001, width = 1.1), fill = "transparent") +
geom_hline(yintercept = 0) +
scale_y_continuous(breaks = c(seq(0.25, 1, by = 0.25), -0.025),
labels = c(c("0.25", "0.50", "0.75", "1.00"), "")) +
facet_grid(~ x, scales = "free_x", space = "free_x", switch = "x") +
ylab("Proportion") + xlab("Sample") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
panel.spacing = unit(.2, units = "mm"),
strip.background = element_blank(),
strip.text.x = element_blank(),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0, 0, 0, 0, unit = "mm"),
axis.text.y = element_text(face = c("plain", "plain", "plain", "plain",
"bold", "bold", "bold"))))
# Composition stacked by individual ROI.
(p2 <-
ggplot(data = df_barplots, aes(x = image_id)) +
geom_bar(aes(fill = factor(cl_all, levels = celltype_order), y = count),
stat = "identity", position = "fill", width = 1) +
scale_fill_manual(name = "Cell types",
values = metadata(sce)$col_celltypes$celltype_all,
breaks = celltype_order) +
geom_hline(yintercept = 0) +
scale_y_continuous(breaks = c(seq(0.25, 1, by = 0.25)),
labels = c(c("0.25", "0.50", "0.75", "1.00"))) +
facet_grid(~ x, scales = "free_x", space = "free_x", switch = "x") +
ylab("Proportion") + xlab("ROI") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
panel.spacing = unit(.1, units = "lines"),
strip.background = element_blank(),
strip.text.x = element_blank(),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0, 0, 0, 0, unit = "mm"),
axis.text.y = element_text(face = c("plain", "plain", "plain", "plain",
"bold", "bold", "bold", "bold"))))
# Stack dendrogram + sample-level + ROI-level composition.
p.het <- ((p.dend / p1 / p2) + plot_layout(guides = "collect", ncol = 1, heights = c(0.25, 1, 1)))
ggsave(plot = p.het, "figures/fig2_het_stacked.pdf", width = 15, heigh = 7, dpi = 300)
```
### TME correlation per image/sample
```{r}
selected_vars <- unique(df_barplots$cl_all)
# Per-ROI cell-type proportions within each sample/condition.
img.props <- df_barplots %>%
group_by(biosample_id, image_id, cond.2) %>%
mutate(prop = prop.table(count)) %>%
select(biosample_id, image_id, cl_all, prop, cond.2)
# For each cell type, compute the SD of its proportion across ROIs within a sample.
output <- lapply(selected_vars, function(celltype) {
sd <- img.props %>% filter(cl_all == celltype) %>% ungroup() %>%
split(.$biosample_id) %>% map(select, prop)
sd.output <- do.call(rbind, lapply(sd, function(x) sd(x$prop)))
colnames(sd.output) <- celltype
sd.output
})
df <- lapply(output, function(x) data.frame(x, biosample_id = row.names(x))) %>%
reduce(full_join, by = "biosample_id")
# Boxplot of intra-sample (between-ROI) variability per cell type.
(p.sd <- df %>% pivot_longer(cols = -c(biosample_id)) %>%
left_join(df_barplots %>% distinct(biosample_id, cond.3, cond.2)) %>%
mutate(cond.new = if_else(!is.na(cond.2), cond.2, cond.3)) %>%
ggplot(aes(x = name, y = (value), fill = name)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(aes(fill = name), colour = "black", pch = 21, alpha = 0.4, show.legend = FALSE) +
scale_fill_manual(name = "Cell types",
values = metadata(sce)$col_celltypes$celltype_all,
breaks = celltype_order) +
scale_color_manual(values = metadata(sce)$col_celltypes$celltype_all) +
scale_x_discrete(limits = make.names(celltype_order), name = "Cell types") +
scale_y_continuous(limits = c(0, 0.5),
name = "SD of cell type proportions between ROIs per sample") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
legend.position = "none"))
# Sample-level variance of cell-type proportions across samples.
img.props <- df_barplots %>% group_by(biosample_id, cl_all) %>%
summarise(countall = sum(count)) %>% mutate(prop = prop.table(countall))
img.props %>% group_by(cl_all) %>%
summarise(var = var(prop), sd = sd(prop)) %>%
ggplot(aes(x = cl_all, y = var)) + geom_boxplot()
ggsave(plot = p.sd, "figures/fig2_het_var.pdf", width = 6, heigh = 5, dpi = 200)
# Median/mean/SD of between-ROI variances per cell type.
df %>% pivot_longer(cols = -c(biosample_id)) %>%
left_join(df_barplots %>% distinct(biosample_id, cond.3, cond.2)) %>%
mutate(cond.new = if_else(!is.na(cond.2), cond.2, cond.3)) %>%
group_by(name) %>%
summarise(median = median(value, na.rm = TRUE),
mean = mean(value, na.rm = TRUE),
sd = sd(value, na.rm = TRUE)) %>%
arrange(desc(median))
```
### Revision: Inferential statistics (ICC and PERMANOVA)
To quantify intra-patient versus inter-patient variability we performed two
formal analyses: (1) per-cell-type linear mixed-effects models partitioning
variance between patient-level and ROI-level effects, summarised by the
intraclass correlation coefficient (ICC) and a likelihood-ratio test for the
patient effect; and (2) a global PERMANOVA on Bray-Curtis distances testing
whether patient identity structures the overall ROI composition.
```{r}
library(lme4)
library(vegan)
library(performance)
# Frequency table of cell types (cl_all) per ROI (image_id), retaining the
# biosample_id so ROIs can be mapped back to patients.
metadata <- as.data.frame(colData(sce))
roi_counts_table <- metadata %>%
group_by(biosample_id, image_id, cl_all) %>%
tally() %>%
ungroup()
# Wide format and row-wise normalization to per-ROI proportions.
df_proportions <- roi_counts_table %>%
pivot_wider(names_from = cl_all, values_from = n, values_fill = 0) %>%
mutate(across(where(is.numeric), ~ .x / rowSums(across(where(is.numeric)))))
# Make cell-type column names syntactically valid (e.g. "B cells" -> "B.cells").
original_cell_types <- setdiff(colnames(df_proportions), c("biosample_id", "image_id"))
valid_cell_types <- make.names(original_cell_types)
colnames(df_proportions)[colnames(df_proportions) %in% original_cell_types] <- valid_cell_types
# --- 1. Per-cell-type ICC and likelihood-ratio test for the patient effect ---
stats_results <- map_df(valid_cell_types, function(ct) {
# Full model: patient as random intercept.
formula_full <- as.formula(paste0(ct, " ~ (1 | biosample_id)"))
model_full <- lmer(formula_full, data = df_proportions, REML = FALSE)
# Null model: no patient effect.
model_null <- lm(as.formula(paste0(ct, " ~ 1")), data = df_proportions)
# Likelihood-ratio test p-value.
lrt <- anova(model_full, model_null)
p_val <- lrt$`Pr(>Chisq)`[2]
# Adjusted ICC = proportion of variance explained by patient.
icc_val <- icc(model_full)$ICC_adjusted
data.frame(CellType = ct, ICC = icc_val, p_value = p_val)
})
# Mean ICC across lineages.
mean_icc <- mean(icc_results$ICC, na.rm = TRUE)
cat("\n--- Quantitative Summary ---\n")
cat("Mean ICC:", round(mean_icc, 3), "\n")
# --- 2. Global PERMANOVA (Bray-Curtis) testing the patient effect ---
comp_matrix <- df_proportions %>% select(all_of(valid_cell_types))
patient_metadata <- df_proportions$biosample_id
set.seed(42)
permanova <- adonis2(comp_matrix ~ patient_metadata, method = "bray")
global_p <- permanova$`Pr(>F)`[1]
r2_val <- permanova$R2[1]
# --- 3. ICC bar plot, coloured by -log10(p) ---
ggplot(stats_results, aes(x = reorder(CellType, ICC), y = ICC)) +
geom_bar(aes(fill = -log10(p_value)), stat = "identity") +
coord_flip() +
labs(title = "Intraclass Correlation (ICC) by Cell Type",
subtitle = paste0("Global PERMANOVA p = ", global_p, " (R^2 = ", round(r2_val, 3), ")"),
x = "Cell Type", y = "ICC (Proportion of Variance by Patient)") +
theme_minimal() +
ylim(0, 1.1)
# --- 4. Formatted summary table (Table S1) ---
library(gt)
table_data <- stats_results %>%
mutate(
CellType = gsub("\\.", " ", CellType),
p_formatted = ifelse(p_value < 0.001, "< 0.001", sprintf("%.3f", p_value)),
ICC = round(ICC, 3)
) %>%
select(CellType, ICC, `p-value` = p_formatted) %>%
arrange(desc(ICC))
results_table <- table_data %>%
gt() %>%
tab_header(
title = "Quantitative Assessment of Intra-patient Homogeneity",
subtitle = "Intraclass Correlation Coefficients (ICC) and Likelihood Ratio Tests per Cell Type"
) %>%
cols_label(
CellType = "Cell Lineage",
ICC = "ICC (Patient Variance)",
`p-value` = "p-value (Patient Effect)"
) %>%
tab_source_note(
source_note = paste0("Global PERMANOVA (Bray-Curtis): p = ",
round(global_p, 4), " | R^2 = ", round(r2_val, 3))
) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_body(
columns = `p-value`,
rows = `p-value` == "< 0.001" | as.numeric(`p-value`) < 0.05
)
) %>%
opt_stylize(color = "gray", style = 1)
results_table
```
## TSNE + total proportions
```{r, fig.width=12, fig.height=10}
## --- t-SNE on a cell subset across all non-nuclear channels ---
selected_markers <- rownames(sce[rowData(sce)$marker_groups %notin% c("nuclei"), ])
# Subsample cells for a tractable embedding.
set.seed(1234)
cur_cells <- sample(seq_len(ncol(sce)), 200000)
tmp_sce <- sce[rownames(sce) %in% selected_markers, cur_cells]
# Run t-SNE on the batch-corrected (fastMNN) representation.
tmp_sce <- runTSNE(tmp_sce,
exprs_values = "fastMNN",
BPPARAM = MulticoreParam(progressbar = T))
# Collapse all immune subsets into a single "Immune" class for the overview.
tmp_sce$cl_all[tmp_sce$cl_all %notin% c("Astrocytes", "Tumor", "Vascular", "unassigned")] <- "Immune"
col_ann <- metadata(tmp_sce)$col_celltypes$celltype_all
col_ann["Immune"] <- col_ann[["CD8Tc"]]
order <- c("Tumor", "Immune", "Vascular", "Astrocytes", "unassigned")
(p.tsne.cohort <- plotReducedDim(tmp_sce, dimred = "TSNE", colour_by = "cl_all", point_size = .15) +
scale_colour_manual(values = col_ann) +
theme(legend.position = "none"))
## --- Overall cohort composition as a single stacked bar with labels ---
df.barplot.all <- colData(sce) %>% as_tibble()
col_ann <- metadata(sce)$col_celltypes$celltype_all
col_ann["Immune"] <- col_ann[["CD8Tc"]]
order <- c("Tumor", "Immune", "Vascular", "Astrocytes", "unassigned")
plot.data <- df.barplot.all %>%
mutate(cl_all = if_else(cl_all %in% c("Astrocytes", "Tumor", "Vascular", "unassigned"),
cl_all, "Immune")) %>%
group_by(cl_all) %>% summarise(count = n()) %>% mutate(prop = prop.table(count))
(p.prop.cohort <- ggplot(data = plot.data, aes(x = 1, y = prop, fill = factor(cl_all, levels = order))) +
geom_bar(stat = "identity") +
ggrepel::geom_text_repel(aes(label = scales::percent(prop, accuracy = 1, suffix = "%")),
position = position_stack(vjust = 0.5),
direction = "y",
xlim = c(1.5, Inf), ylim = c(-Inf, Inf),
size = 5, hjust = 0,
segment.size = .7, segment.alpha = .5,
segment.linetype = "dotted", box.padding = .4,
segment.curvature = -0.1, segment.ncp = 3, segment.angle = 20) +
coord_cartesian(clip = "off") +
ylab("Proportion") +
scale_fill_manual(name = "Celltypes", values = col_ann, breaks = order) +
theme_minimal() +
theme(legend.position = "right",
plot.margin = unit(c(0, 7, 0, 0), "cm"),
axis.title.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
panel.background = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.x = element_blank()))
p.comb <- p.tsne.cohort | p.prop.cohort + plot_layout(widths = c(1.5, 1))
ggsave(plot = p.comb, "figures/fig2_het_stacked_cohort.png", width = 12, heigh = 7, dpi = 300)
```
### Suppl: Marker expression on t-SNE
```{r, fig.width=20, fig.height=20}
# One t-SNE panel per marker, coloured by z-scored expression.
plot_list <- lapply(sort(rownames(tmp_sce)), function(x) {
p <- plotReducedDim(tmp_sce, dimred = "TSNE", colour_by = x,
by_exprs_values = "scaled", point_size = 0.2)
p + scale_colour_gradient2(low = "#313695", mid = "white", high = "#A50026",
limits = c(-4, 4), oob = scales::squish,
name = "z-score expression") +
ggtitle(x) +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "none",
axis.ticks = element_blank(),
axis.text = element_blank())
})
legend <- get_legend(plot_list[[1]] + theme(legend.position = "right"))
p.wrap <- wrap_plots(c(plot_list, list(legend))) +
plot_layout(axes = "collect", axis_titles = "collect")
ggsave(plot = p.wrap, "figures/suppfig_all_expression.png", width = 20, heigh = 20, dpi = 300)
```
## Suppl: T cell + Myeloid subsets
### t-SNE
#### Lymphoid / T cells
```{r}
## T cells: t-SNE on T-cell clustering markers.
selected_markers <- rownames(sce[rowData(sce)$cl_tcells, ])
# Subset to T-cell populations.
tmp_sce <- sce[, sce$cl_all %in% c("CD4Tc", "CD8Tc", "DPTc", "Treg")]
# Subsample and embed.
set.seed(1234)
cur_cells <- sample(seq_len(ncol(tmp_sce)), 100000)
tmp_sce <- tmp_sce[rownames(tmp_sce) %in% selected_markers, cur_cells]
tmp_sce <- runTSNE(tmp_sce, exprs_values = "fastMNN", BPPARAM = MulticoreParam(progressbar = T))
col_ann <- metadata(tmp_sce)$col_celltypes$celltype_all
# Cell-type embedding.
(p.tsne.lymphoid <- plotReducedDim(tmp_sce, dimred = "TSNE", colour_by = "cl_all_spec",
point_size = .2, text_by = "cl_all_spec") +
scale_colour_manual(values = col_ann_hm) +
ggtitle("Cell types") +
guides(color = guide_legend(override.aes = list(size = 5), title = "Cell types")) +
theme(plot.title = element_text(hjust = 0.5),
axis.ticks = element_blank(),
axis.text = element_blank(),
legend.position = "none")) +
coord_fixed()
legend.celltypes_lymphoid <- get_legend(p.tsne.lymphoid + theme(legend.position = "right"))
# Per-marker expression panels.
plot_list_lymphoid <- lapply(sort(rownames(tmp_sce)), function(x) {
p <- plotReducedDim(tmp_sce, dimred = "TSNE", colour_by = x,
by_exprs_values = "scaled", point_size = 0.2)
p + scale_colour_gradient2(low = "#313695", mid = "white", high = "#A50026",
limits = c(-4, 4), oob = scales::squish,
name = "z-score expression") +
ggtitle(x) +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "none",
axis.ticks = element_blank(),
axis.text = element_blank()) +
coord_fixed()
})
legend_lymphoid <- get_legend(plot_list_lymphoid[[1]] + theme(legend.position = "right"))
```
#### Myeloid
```{r}
## Myeloid: t-SNE on myeloid clustering markers.
selected_markers <- rownames(sce[rowData(sce)$cl_myeloid, ])
# Subset to myeloid + neutrophil populations.
tmp_sce <- sce[, sce$cl_immune_a %in% c("Myeloid", "Neutrophils")]
# Subsample and embed.
set.seed(1234)
cur_cells <- sample(seq_len(ncol(tmp_sce)), 100000)
tmp_sce <- tmp_sce[rownames(tmp_sce) %in% selected_markers, cur_cells]
tmp_sce <- runTSNE(tmp_sce, exprs_values = "fastMNN", BPPARAM = MulticoreParam(progressbar = T))
col_ann <- metadata(tmp_sce)$col_celltypes$celltype_all
# Cell-type embedding.
(p.tsne.myeloid <- plotReducedDim(tmp_sce, dimred = "TSNE", colour_by = "cl_all_spec",
point_size = .15, text_by = "cl_all_spec",
text_colour = "black", text_size = 4) +
scale_colour_manual(values = col_ann_hm) +
ggtitle("Cell types") +
guides(color = guide_legend(override.aes = list(size = 5), title = "Cell types")) +
theme(plot.title = element_text(hjust = 0.5),
axis.ticks = element_blank(),
axis.text = element_blank(),
legend.position = "none")) +
coord_fixed()
legend.celltypes_myeloid <- get_legend(p.tsne.myeloid + theme(legend.position = "right"))
# Per-marker expression panels.
plot_list_myeloid <- lapply(sort(rownames(tmp_sce)), function(x) {
p <- plotReducedDim(tmp_sce, dimred = "TSNE", colour_by = x,
by_exprs_values = "scaled", point_size = 0.2)
p + scale_colour_gradient2(low = "#313695", mid = "white", high = "#A50026",
limits = c(-4, 4), oob = scales::squish,
name = "z-score expression") +
ggtitle(x) +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "none",
axis.ticks = element_blank(),
axis.text = element_blank()) +
coord_fixed()
})
legend_myeloid <- get_legend(plot_list_myeloid[[1]] + theme(legend.position = "right"))
```
### Heatmaps
#### T cells
```{r}
cluster_name <- "cl_all_spec"
selected_markers <- rownames(sce[rowData(sce)$cl_tcells, ])
# Subset to T-cell populations.
tmp_sce <- sce[rownames(sce) %in% selected_markers,
sce$cl_all %in% c("CD4Tc", "CD8Tc", "DPTc", "Treg")]
# Mean expression per cluster and scaled assays.
mean_tmp.sce <- aggregateAcrossCells(tmp_sce, ids = tmp_sce[[cluster_name]], statistics = "mean")
assay(mean_tmp.sce, "exprs") <- asinh(counts(mean_tmp.sce) / 1)
assay(mean_tmp.sce, "scaled") <- t(scale(t(assay(mean_tmp.sce, "exprs"))))
assay(mean_tmp.sce, "scaled_ztm") <- t(apply(assay(mean_tmp.sce, "exprs"), 1, scales::rescale))
# Heatmap body: zero-to-max scaled expression.
hm_body <- t(assay(mean_tmp.sce, "scaled_ztm"))
row_anno <- colData(mean_tmp.sce) %>% as.data.frame %>% select(all_of(cluster_name), ncells)
row_anno[[cluster_name]] <- as.character(row_anno[[cluster_name]])
# Viridis scale for zero-to-max values.
col_1 <- viridis(100)
ha_right <- HeatmapAnnotation(
cluster_name = anno_simple(row_anno[[cluster_name]], border = TRUE, col = col_ann_hm),
names = anno_text(str_wrap(row_anno[[cluster_name]], width = 12), just = "left"),
annotation_label = "",
annotation_name_rot = 90,
which = "row")
hm <- Heatmap(hm_body,
name = "zero-to-max",
col = col_1,
clustering_method_columns = "complete",
clustering_method_rows = "complete",
column_split = fct_recode(rowData(mean_tmp.sce)$marker_groups,
"Cellstate" = "cellstate",
"Lineage" = "lineage"),
cluster_column_slices = F,
show_row_names = F,
right_annotation = c(ha_right))
hm <- draw(hm)
p.hm.tcells <- grid.grabExpr(draw(hm))
```
#### Myeloid
```{r}
cluster_name <- "cl_all_spec"
selected_markers <- rownames(sce[rowData(sce)$cl_myeloid, ])
# Subset to myeloid + neutrophil populations.
tmp_sce <- sce[rownames(sce) %in% selected_markers,
sce$cl_immune_a %in% c("Myeloid", "Neutrophils")]
# Mean expression per cluster and scaled assays.
mean_tmp.sce <- aggregateAcrossCells(tmp_sce, ids = tmp_sce[[cluster_name]], statistics = "mean")
assay(mean_tmp.sce, "exprs") <- asinh(counts(mean_tmp.sce) / 1)
assay(mean_tmp.sce, "scaled") <- t(scale(t(assay(mean_tmp.sce, "exprs"))))
assay(mean_tmp.sce, "scaled_ztm") <- t(apply(assay(mean_tmp.sce, "exprs"), 1, scales::rescale))
# Heatmap body: zero-to-max scaled expression.
hm_body <- t(assay(mean_tmp.sce, "scaled_ztm"))
row_anno <- colData(mean_tmp.sce) %>% as.data.frame %>% select(all_of(cluster_name), ncells)
row_anno[[cluster_name]] <- as.character(row_anno[[cluster_name]])
col_1 <- viridis(100)
ha_right <- HeatmapAnnotation(
cluster_name = anno_simple(row_anno[[cluster_name]], border = TRUE, col = col_ann_hm),
names = anno_text(str_wrap(row_anno[[cluster_name]], width = 12), just = "left"),
annotation_label = "",
annotation_name_rot = 90,
which = "row")
hm <- Heatmap(hm_body,
name = "zero-to-max",
col = col_1,
clustering_method_columns = "complete",
clustering_method_rows = "complete",
column_split = fct_recode(rowData(mean_tmp.sce)$marker_groups,
"Cellstate" = "cellstate",
"Lineage" = "lineage"),
cluster_column_slices = F,
show_row_names = F,
right_annotation = c(ha_right))
hm <- draw(hm)
p.hm.myeloid <- grid.grabExpr(draw(hm))
```
### Barplots
#### T cells
```{r}
# T-cell colData and stacking order.
df.col <- colData(sce) %>% as_tibble() %>%
filter(cl_all %in% c("CD4Tc", "CD8Tc", "DPTc", "Treg"))
celltype_order <- unique(df.col$cl_all_spec)
# Overall T-cell subset composition (single labelled stacked bar).
plot.data <- df.col %>%
group_by(cl_all_spec) %>% summarise(count = n()) %>% mutate(prop = prop.table(count))
(p.prop.tcells <- ggplot(data = plot.data, aes(x = 1, y = prop, fill = cl_all_spec)) +
geom_bar(stat = "identity") +
ggrepel::geom_text_repel(aes(label = scales::percent(prop, accuracy = 1, suffix = "%")),
position = position_stack(vjust = 0.5),
direction = "y",
xlim = c(1.5, Inf), ylim = c(-Inf, Inf),
size = 5, hjust = 0,
segment.size = .7, segment.alpha = .5,
segment.linetype = "dotted", box.padding = .4,
segment.curvature = -0.1, segment.ncp = 3, segment.angle = 20) +
coord_cartesian(clip = "off") +
ylab("Proportion") +
scale_fill_manual(name = "Cell types", values = col_ann_hm) +
theme_minimal() +
theme(legend.position = "right",
plot.margin = unit(c(0, 7, 0, 0), "cm"),
legend.box.margin = margin(0, 0, 0, 25),
axis.title.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
panel.background = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.x = element_blank()))
# Cluster samples by T-cell subset composition and build the dendrogram.
df_clust <- df.col %>%
group_by(biosample_id, cl_all_spec) %>%
summarise(count = n()) %>% mutate(prop = prop.table(count)) %>% ungroup()
df_clust_wide <- df_clust %>%
pivot_wider(id_cols = biosample_id, names_from = cl_all_spec, values_from = prop) %>%
replace(is.na(.), 0)
df_clust_m <- as.matrix(df_clust_wide)
rownames(df_clust_m) <- df_clust_wide$biosample_id
hc <- hclust(dist(df_clust_m[, -1], method = "euclidian"), "ward.D2")
dend <- as.dendrogram(hc)
dend_data <- dendro_data(dend, type = "rectangle")
(p.dend <- ggplot(dend_data$segments) +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend)) +
geom_text(data = dend_data$labels, aes(x, y, label = label), hjust = 1, angle = 90) +
scale_x_continuous(expand = c(0, 0.5)) +
scale_y_continuous(expand = c(0, 0)) +
theme_void() +
theme(plot.margin = margin(0, 0, 0, 0, unit = "mm")))
# Per-ROI counts with sample/condition annotation, ordered by clustering.
df_barplots <- df.col %>%
group_by(image_id, biosample_id, pid, cond.3, cond.2, cond.4, cl_all_spec) %>%
summarise(count = n()) %>% ungroup()
df_barplots$biosample_id <- factor(df_barplots$biosample_id, levels = hc$labels[hc$order])
df_barplots$pid <- as.factor(df_barplots$pid)
x <- dend_data$labels[c("x", "label")] %>% as_tibble() %>% rename(biosample_id = label)
df_barplots <- df_barplots %>% left_join(x, by = "biosample_id")
# Stacked composition by sample, with a condition.2 annotation track.
(p1 <-
ggplot(data = df_barplots) +
geom_bar(aes(fill = factor(cl_all_spec, levels = celltype_order), x = x, y = count),
stat = "identity", position = "fill", width = 1.1) +
scale_fill_manual(name = "Celltypes", values = col_ann_hm, breaks = celltype_order) +
geom_tile(aes(x = x, y = 0, height = 0.001, width = 1.1), fill = "transparent") +
new_scale_fill() +
geom_tile(aes(x = x, y = -0.055, height = 0.025, width = 1.1, fill = cond.2)) +
scale_fill_manual(name = "Condition 2",
values = metadata(sce)$col_clinical$cond.2,
na.value = "transparent") +
scale_y_continuous(breaks = c(seq(0.25, 1, by = 0.25), -0.025),
labels = c(c("0.25", "0.50", "0.75", "1.00"), "")) +
facet_grid(~ x, scales = "free_x", space = "free_x", switch = "x") +
ylab("Proportion") + xlab("Sample") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
panel.spacing = unit(.2, units = "mm"),
strip.background = element_blank(),
strip.text.x = element_blank(),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0, 0, 0, 0, unit = "mm"),
axis.text.y = element_text(face = c("plain", "plain", "plain", "plain",
"bold", "bold", "bold"))))
(p.dend / p1) + plot_layout(guides = "collect", ncol = 1, heights = c(0.25, 1))
```
#### Myeloid
```{r}
# Myeloid colData and stacking order.
df.col <- colData(sce) %>% as_tibble() %>%
filter(cl_immune_a %in% c("Myeloid", "Neutrophils"))
celltype_order <- unique(df.col$cl_all_spec)
# Overall myeloid subset composition (single labelled stacked bar).
plot.data <- df.col %>%
group_by(cl_all_spec) %>% summarise(count = n()) %>% mutate(prop = prop.table(count))
(p.prop.myeloid <- ggplot(data = plot.data, aes(x = 1, y = prop, fill = cl_all_spec)) +
geom_bar(stat = "identity") +
ggrepel::geom_text_repel(aes(label = scales::percent(prop, accuracy = 1, suffix = "%")),
position = position_stack(vjust = 0.5),
direction = "y",
xlim = c(1.5, Inf), ylim = c(-Inf, Inf),
size = 5, hjust = 0,
segment.size = .7, segment.alpha = .5,
segment.linetype = "dotted", box.padding = .4,
segment.curvature = -0.1, segment.ncp = 3, segment.angle = 20) +
coord_cartesian(clip = "off") +
ylab("Proportion") +
scale_fill_manual(name = "Cell types", values = col_ann_hm) +
theme_minimal() +
theme(legend.position = "right",
plot.margin = unit(c(0, 7, 0, 0), "cm"),
legend.box.margin = margin(0, 0, 0, 25),
axis.title.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
panel.background = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major.x = element_blank()))
# Cluster samples by myeloid subset composition and build the dendrogram.
df_clust <- df.col %>%
group_by(biosample_id, cl_all_spec) %>%
summarise(count = n()) %>% mutate(prop = prop.table(count)) %>% ungroup()
df_clust_wide <- df_clust %>%
pivot_wider(id_cols = biosample_id, names_from = cl_all_spec, values_from = prop) %>%
replace(is.na(.), 0)
df_clust_m <- as.matrix(df_clust_wide)
rownames(df_clust_m) <- df_clust_wide$biosample_id
hc <- hclust(dist(df_clust_m[, -1], method = "euclidian"), "ward.D2")
dend <- as.dendrogram(hc)
dend_data <- dendro_data(dend, type = "rectangle")
(p.dend <- ggplot(dend_data$segments) +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend)) +
geom_text(data = dend_data$labels, aes(x, y, label = label), hjust = 1, angle = 90) +
scale_x_continuous(expand = c(0, 0.5)) +
scale_y_continuous(expand = c(0, 0)) +
theme_void() +
theme(plot.margin = margin(0, 0, 0, 0, unit = "mm")))
# Per-ROI counts with sample/condition annotation, ordered by clustering.
df_barplots <- df.col %>%
group_by(image_id, biosample_id, pid, cond.3, cond.2, cond.4, cl_all_spec) %>%
summarise(count = n()) %>% ungroup()
df_barplots$biosample_id <- factor(df_barplots$biosample_id, levels = hc$labels[hc$order])
df_barplots$pid <- as.factor(df_barplots$pid)
x <- dend_data$labels[c("x", "label")] %>% as_tibble() %>% rename(biosample_id = label)
df_barplots <- df_barplots %>% left_join(x, by = "biosample_id")
# Stacked composition by sample, with a condition.2 annotation track.
(p1 <-
ggplot(data = df_barplots) +
geom_bar(aes(fill = factor(cl_all_spec, levels = celltype_order), x = x, y = count),
stat = "identity", position = "fill", width = 1.1) +
scale_fill_manual(name = "Celltypes", values = col_ann_hm, breaks = celltype_order) +
geom_tile(aes(x = x, y = 0, height = 0.001, width = 1.1), fill = "transparent") +
new_scale_fill() +
geom_tile(aes(x = x, y = -0.055, height = 0.025, width = 1.1, fill = cond.2)) +
scale_fill_manual(name = "Condition 2",
values = metadata(sce)$col_clinical$cond.2,
na.value = "transparent") +
scale_y_continuous(breaks = c(seq(0.25, 1, by = 0.25), -0.025),
labels = c(c("0.25", "0.50", "0.75", "1.00"), "")) +
facet_grid(~ x, scales = "free_x", space = "free_x", switch = "x") +
ylab("Proportion") + xlab("Sample") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
panel.spacing = unit(.2, units = "mm"),
strip.background = element_blank(),
strip.text.x = element_blank(),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0, 0, 0, 0, unit = "mm"),
axis.text.y = element_text(face = c("plain", "plain", "plain", "plain",
"bold", "bold", "bold"))))
(p.dend / p1) + plot_layout(guides = "collect", ncol = 1, heights = c(0.25, 1))
```
### Assemble
```{r, fig.width=19, fig.height=20}
# Lymphoid supplementary figure: heatmap + cell-type t-SNE + proportions (top),
# per-marker expression panels (bottom).
(( wrap_plots(p.hm.tcells) | p.tsne.lymphoid + coord_fixed() | p.prop.tcells ) +
plot_layout(widths = c(1.5, 1, .1)) ) /
(wrap_plots(c(plot_list_lymphoid, list(legend_lymphoid))) +
plot_layout(axes = "collect", axis_titles = "collect")) +
plot_annotation(tag_levels = list(c("A", "B", "C", "D"))) +
plot_layout(height = c(.25, 1))
ggsave(filename = "figures/suppfig_lymphoid.png", width = 19, heigh = 20, dpi = 300)
# Myeloid supplementary figure (same layout).
(( wrap_plots(p.hm.myeloid) | p.tsne.myeloid + coord_fixed() | p.prop.myeloid) +
plot_layout(widths = c(1.5, 1, .15)) ) /
(wrap_plots(c(plot_list_myeloid, list(legend_myeloid))) +
plot_layout(axes = "collect", axis_titles = "collect")) +
plot_annotation(tag_levels = list(c("A", "B", "C", "D"))) +
plot_layout(height = c(.25, 1))
ggsave(filename = "figures/suppfig_myeloid.png", width = 19, heigh = 18, dpi = 300)
```
## Pearson correlation + spatial interaction
```{r, fig.width=10}
# Permutation-based neighbourhood interaction test (imcRtools) per sample.
out <- testInteractions(sce[, sce$cl_all_spec != "unassigned"],
group_by = "biosample_id",
label = "cl_all_spec",
colPairName = "neighborhood",
method = "classic",
BPPARAM = MulticoreParam(workers = 1, RNGseed = 221029))
head(out)
rois <- length(unique(out$group_by))
# Summarise significant interaction/avoidance counts per cell-type pair.
df.out <- out %>% as_tibble() %>%
group_by(from_label, to_label) %>%
summarize(sum_sigval = sum(sigval, na.rm = TRUE)) %>%
mutate(prop_sigval = abs(sum_sigval / rois))
# (Standalone z-scored interaction heatmap for the whole cohort.)
p <- df.out %>%
ggplot(aes(from_label, to_label)) +
geom_tile(aes(fill = scale(sum_sigval))) +
scale_fill_gradient2(low = "#313695", mid = "white", high = "#A50026",
limits = c(-2, 2), oob = scales::squish,
name = "z-score expression") +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
ggtitle("Entire cohort")
# Pearson correlation of cell-type proportions across samples.
df.cor <- colData(sce) %>% as_tibble() %>%
filter(cl_all_spec != "unassigned") %>%
group_by(biosample_id, cl_all_spec) %>%
summarise(count = n()) %>% mutate(prop = prop.table(count)) %>%
pivot_wider(id_cols = biosample_id, names_from = cl_all_spec, values_from = prop) %>%
column_to_rownames(var = "biosample_id") %>% replace(is.na(.), 0) %>% cor()
# Reorder a correlation matrix by hierarchical clustering of (1 - r)/2.
reorder_corr.matrix <- function(corr.matrix) {
distance <- as.dist((1 - corr.matrix) / 2)
hc <- hclust(distance)
corr.matrix[hc$order, hc$order]
}
df.cor <- reorder_corr.matrix(df.cor)
# Keep only the lower triangle and reshape to long format.
df.cor[upper.tri(df.cor)] <- NA
df.cor.long <- reshape2::melt(df.cor, rm.na = T)
order <- levels(df.cor.long$Var1)
# Combined plot: tiles = Pearson correlation; points = neighbourhood
# interaction/avoidance (size = % significant samples, fill = z-scored direction).
library(colorspace)
library(ggnewscale)
df.cor.long %>%
left_join(df.out, by = c("Var1" = "from_label", "Var2" = "to_label")) %>%
filter(!is.na(value)) %>%
mutate(sum_sigval_cat = case_when(sum_sigval == 0 ~ "None",
sum_sigval > 0 ~ "Interaction",
sum_sigval < 0 ~ "Avoidance")) %>%
ggplot(aes(factor(Var1, levels = order), factor(Var2, levels = order))) +
geom_tile(aes(fill = value), color = "white") +
scale_fill_gradient2(low = "#313695", mid = "white", high = "#A50026",
limits = c(-1, 1), midpoint = 0, na.value = "transparent",
name = "Pearson Correlation of celltype\nproportions in samples") +
new_scale_fill() +
geom_point(shape = 21,
aes(size = prop_sigval,
fill = scale(sum_sigval, center = T, scale = T),
alpha = prop_sigval == 0),
colour = "black") +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
limits = c(-1, 1), na.value = "transparent",
oob = scales::squish,
name = "Avoidance/Interaction\n(z-score)",
breaks = c(-2, -1, 0, 1, 2),
labels = c("Avoidance", -1, 0, 1, "Interaction")) +
scale_size_continuous(name = "%Samples with significant\ninteraction/avoidance",
breaks = c(0.25, 0.5, 0.75, 1.00)) +
scale_alpha_manual(values = c(1, 0)) +
guides(alpha = "none") +
theme_void() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, size = 12, hjust = 1),
axis.text.y = element_text(size = 12, hjust = 1)) +
coord_fixed()
ggsave(filename = "figures/fig_correlation_interaction.pdf", width = 10, height = 10, dpi = 200)
```
# --------------------------------
# FIGURE 3
Survival analysis of the entire cohort.
## Proportions vs. survival
```{r, fig.width = 7.5, fig.height=6}
# Clinical/survival data restricted to the cohort.
coldata <- colData(sce) %>% as_tibble()
survData <- metadata(sce)$clinical %>%
select(biosample_id, status, status_pfs, days_dg, days_dos, days_dos_pfs)
survData <- coldata %>% distinct(biosample_id) %>% left_join(survData)
sample.ids <- coldata %>% distinct(biosample_id)
# Per-sample cell-type proportions (tumor excluded), used as continuous predictors.
celltype.props <- coldata %>% filter(cl_all_spec != "Tumor") %>%
group_by(biosample_id, cl_all_spec) %>% summarise(count = n()) %>%
mutate(prop = prop.table(count)) %>% ungroup() %>% select(-count) %>%
pivot_wider(names_from = "cl_all_spec", values_from = "prop") %>% replace(is.na(.), 0)
variables <- celltype.props %>% ungroup() %>% select(-biosample_id) %>% colnames()
# --- Univariate Cox: overall survival (OS) ---
time <- "days_dos"
event <- "status"
survData.subset <- survData %>% filter(biosample_id %in% sample.ids$biosample_id) %>%
filter(!is.na(!!rlang::sym(time)) & !is.na(!!rlang::sym(event)))
df <- survData.subset %>% left_join(celltype.props) %>%
mutate(across(all_of(variables), ~ scale(.x)[, 1])) # z-score each predictor
res.os <- run.coxph(df, variables, var.time = time, var.status = event, univariate = T)
res.os$outcome <- "OS"
res.os <- res.os %>% mutate(p.fdr = p.adjust(p.value, method = "fdr"))
# --- Univariate Cox: progression-free survival (PFS) ---
time <- "days_dos_pfs"
event <- "status_pfs"
survData.subset <- survData %>% filter(biosample_id %in% sample.ids$biosample_id) %>%
filter(!is.na(!!rlang::sym(time)) & !is.na(!!rlang::sym(event)))
df <- survData.subset %>% left_join(celltype.props) %>%
mutate(across(all_of(variables), ~ scale(.x)[, 1]))
res.pfs <- run.coxph(df, variables, var.time = time, var.status = event, univariate = T)
res.pfs$outcome <- "PFS"
res.pfs <- res.pfs %>% mutate(p.fdr = p.adjust(p.value, method = "fdr"))
res.prop.all <- rbind(res.os, res.pfs)
# Colour scheme by outcome and cell-type ordering from the Figure 2 heatmap.
col_ann <- setNames(dittoColors()[seq_along(unique(res.prop.all$outcome))],
unique(res.prop.all$outcome))
order <- rev(mean_tmp.sce$cl_all_spec[row_order(hm.celltypes) &
mean_tmp.sce$cl_all_spec %in% res.prop.all$term])
# Panel 0: cell-type colour key.
p0 <- res.prop.all %>% ggplot(aes(x = 0, width = 1, y = term, fill = term, height = 0.8)) +
geom_tile(show.legend = F) +
scale_fill_manual(values = col_ann_hm) +
scale_y_discrete(limits = order) +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.spacing = unit(.2, units = "mm"),
strip.background = element_blank(),
strip.text.x = element_blank(),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
text = element_text(size = 15))
# Panel 1: hazard ratios with confidence intervals.
(p1 <- res.prop.all %>% ggplot(aes(y = term, color = outcome)) +
geom_point(aes(x = (estimate)), shape = 18, size = 5, position = position_dodge(0.7)) +
geom_errorbarh(aes(xmin = (conf.low), xmax = (conf.high)), height = 0.25,
position = position_dodge(0.7)) +
geom_vline(xintercept = (1), color = "black", linetype = "dashed", cex = .5, alpha = 0.5) +
scale_x_continuous(breaks = c(0, 1, 2), labels = c(0, 1, 2), name = "HR") +
scale_y_discrete(limits = order) +
scale_color_manual(values = col_ann) +
theme_minimal() +
theme(legend.position = "none",
axis.title.y = element_blank(),
axis.text.y = element_blank(),
text = element_text(size = 15)))
# Panel 2: -log10(p) bars with FDR-significant markers.
(p2 <- ggplot(data = res.prop.all, aes(y = term, fill = outcome)) +
# Invisible dummy layer to inject the "* FDR < 0.1" legend key.
geom_point(aes(x = -log10(p.value), color = ""), show.legend = TRUE) +
scale_color_manual(name = "* FDR < 0.1", values = "transparent",
guide = guide_legend(override.aes = list(color = NA))) +
geom_bar(aes(x = -log10(p.value)), position = position_dodge(0.7), stat = "identity") +
geom_text(data = res.prop.all %>% filter(p.fdr < 0.1),
aes(x = -log10(p.value), label = "*", group = outcome),
position = position_dodge(0.7), hjust = -0.3, size = 4.5) +
geom_vline(xintercept = -log10(0.05), color = "black", linetype = "dashed",
cex = .5, alpha = 0.5) +
scale_fill_manual(values = col_ann, name = "Outcome") +
scale_y_discrete(limits = (rev(mean_tmp.sce$cl_all_spec[row_order(hm) &
mean_tmp.sce$cl_all_spec %in% res.prop.all$term]))) +
scale_x_continuous(breaks = NULL, name = "-log10(p value)",
sec.axis = dup_axis(breaks = c(-log10(0.05)),
labels = c("p = 0.05"), name = NULL)) +
coord_cartesian(xlim = c(0, max(-log10(res.prop.all$p.value)) + 0.5)) +
theme_minimal() +
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_text(face = "bold"),
axis.ticks.y = element_blank(),
text = element_text(size = 15)))
p.prop.os <- (p0 + ggtitle("Cell type proportions - Cox survival model") + p1 + p2) +
plot_layout(widths = c(0.5, 6, 2))
p.prop.os
ggsave(filename = "figures/suppfig_CoxProportions.pdf", width = 7.5, heigh = 6, dpi = 300)
```
### Table
```{r}
library(gt)
# Format Cox results (HR with CI, raw and FDR-adjusted p-values).
make_tbl <- function(df) {
df %>%
mutate(
`HR (95% CI)` = case_when(
estimate > 1000 | estimate < 0.001 ~
sprintf("%.2e (%.2e-%.2e)", estimate, conf.low, conf.high),
TRUE ~
sprintf("%.2f (%.2f-%.2f)", estimate, conf.low, conf.high)
),
`p-value` = ifelse(p.value < 0.001, "<0.001", sprintf("%.3f", p.value)),
`p-adj (FDR)` = ifelse(p.fdr < 0.001, "<0.001", sprintf("%.3f", p.fdr))
) %>%
select(term, `HR (95% CI)`, `p-value`, `p-adj (FDR)`, p.value, p.fdr)
}
# Join OS and PFS side by side, keeping raw p-values for conditional styling.
tbl_data <- left_join(
res.os %>% make_tbl() %>% rename_with(~ paste0(.x, ".os"), -term),
res.pfs %>% make_tbl() %>% rename_with(~ paste0(.x, ".pfs"), -term),
by = "term"
) %>%
slice(match(rev(order), term))
tbl_data %>%
select(term,
`HR (95% CI).os`, `p-value.os`, `p-adj (FDR).os`,
`HR (95% CI).pfs`, `p-value.pfs`, `p-adj (FDR).pfs`) %>%
gt() %>%
cols_label(
term = "Variable",
`HR (95% CI).os` = "HR (95% CI)",
`p-value.os` = "p-value",
`p-adj (FDR).os` = "Adjusted p-value (FDR)",
`HR (95% CI).pfs` = "HR (95% CI)",
`p-value.pfs` = "p-value",
`p-adj (FDR).pfs` = "Adjusted p-value (FDR)"
) %>%
tab_spanner(label = md("**OS**"), columns = ends_with(".os")) %>%
tab_spanner(label = md("**Local PFS**"), columns = ends_with(".pfs")) %>%
tab_header(title = "Univariate Cox Proportional Hazards Model") %>%
tab_style(style = cell_text(weight = "bold"),
locations = cells_body(columns = `p-value.os`, rows = tbl_data$p.value.os < 0.05)) %>%
tab_style(style = cell_text(weight = "bold"),
locations = cells_body(columns = `p-value.pfs`, rows = tbl_data$p.value.pfs < 0.05)) %>%
tab_style(style = cell_text(weight = "bold"),
locations = cells_body(columns = `p-adj (FDR).os`, rows = tbl_data$p.fdr.os < 0.1)) %>%
tab_style(style = cell_text(weight = "bold"),
locations = cells_body(columns = `p-adj (FDR).pfs`, rows = tbl_data$p.fdr.pfs < 0.1)) %>%
tab_style(style = cell_text(weight = "bold"), locations = cells_column_labels()) %>%
tab_footnote(
footnote = "p-values adjusted using the Benjamini-Hochberg false discovery rate (FDR) method.",
locations = cells_column_labels(columns = ends_with("(FDR)"))
) %>%
tab_options(table.font.size = 12, data_row.padding = px(4))
```
## Lcross function
### CD8Tc/DPTc vs. Tumor
```{r}
# Define the cell types used in the cross-type L-function analysis: combine
# CD8 and double-positive T cells, compared against tumor cells.
sce$celltype_tcells <- sce$celltype_a
sce$celltype_tcells[sce$cl_all == "CD8Tc"] <- "CD8Tc/DPTc"
sce$celltype_tcells[sce$cl_all == "DPTc"] <- "CD8Tc/DPTc"
celltypeA <- "CD8Tc/DPTc"
celltypeB <- "Tumor"
cluster.lvl <- "celltype_tcells"
type.select <- "distal"
# Tag cells with their Lcross role.
sce$Lcross_celltype <- ""
sce[, sce[[cluster.lvl]] == celltypeA]$Lcross_celltype <- celltypeA
sce[, sce[[cluster.lvl]] == celltypeB]$Lcross_celltype <- celltypeB
cur_list <- list()
sce.sub <- sce[, sce[[cluster.lvl]] %in% c(celltypeA, celltypeB)]
# For each image, compute the inhomogeneous cross-L function (tumor -> T cells)
# and integrate the L(r)-r curve over a proximal and a distal radius window.
for (i in unique(sce.sub$image_id)) {
cur_sce <- sce[, sce$image_id == i]
celltypeA_count <- ncol(cur_sce[, cur_sce$Lcross_celltype == celltypeA])
celltypeB_count <- ncol(cur_sce[, cur_sce$Lcross_celltype == celltypeB])
# Require at least 10 cells of each type in the image.
if (celltypeA_count >= 10 & celltypeB_count >= 10) {
cur_sce <- cur_sce[, cur_sce[[cluster.lvl]] %in% c(celltypeA, celltypeB)]
cur_sce$Lcross_celltype <- ""
cur_sce[, cur_sce[[cluster.lvl]] == celltypeA]$Lcross_celltype <- celltypeA
cur_sce[, cur_sce[[cluster.lvl]] == celltypeB]$Lcross_celltype <- celltypeB
# Build a point pattern within a concave hull of the cell coordinates.
cdat <- as.data.frame(colData(cur_sce)[, c("Pos_X", "Pos_Y")])
pt <- cdat %>% st_as_sf(coords = c("Pos_X", "Pos_Y"))
pts <- ppp(x = cdat$Pos_X,
y = cdat$Pos_Y,
marks = as.factor(cur_sce$Lcross_celltype),
window = as.owin(concaveman(points = pt, concavity = 2)))
grid_area <- st_area(concaveman(points = pt, concavity = 2)) / 10^6
# Inhomogeneous cross-L function.
Lfun <- Lcross.inhom(pts, from = celltypeB, to = celltypeA)
radius_close <- 25
radius_far <- 150
# Integrated L(r)-r over proximal (10-25) and distal (25-150) windows, plus
# cell densities per mm^2.
cur_out <- as.data.frame(matrix(nrow = 2, ncol = 4))
colnames(cur_out) <- c("type", "integral",
paste0(celltypeA, "_density"), paste0(celltypeB, "_density"))
cur_out[1, ] <- c("proximal", sum(Lfun$iso[10:radius_close] - Lfun$r[10:radius_close]),
round(celltypeA_count / grid_area), round(celltypeB_count / grid_area))
cur_out[2, ] <- c("distal", sum(Lfun$iso[25:radius_far] - Lfun$r[25:radius_far]),
round(celltypeA_count / grid_area), round(celltypeB_count / grid_area))
cur_list[[i]] <- cur_out
}
}
# Combine per-image results into one data frame.
df <- as.data.frame(do.call(rbind, cur_list))
df$image_id <- str_sub(rownames(df), start = 1, end = -3)
df$integral <- as.numeric(df$integral)
df[[paste0(celltypeA, "_density")]] <- as.numeric(df[[paste0(celltypeA, "_density")]])
df[[paste0(celltypeB, "_density")]] <- as.numeric(df[[paste0(celltypeB, "_density")]])
df.cd8tcells <- df
# Select images with above-median T-cell density and extreme integrals
# (most dispersed vs. most localized).
df <- df.cd8tcells
median.density.target <- median(df[[paste0(celltypeA, "_density")]])
imgids.max <- df %>%
filter(!!rlang::sym(paste0(celltypeA, "_density")) > median.density.target) %>%
filter(type == type.select) %>%
slice_max(order_by = integral, n = 20) %>% pull(image_id)
imgids.min <- df %>%
filter(!!rlang::sym(paste0(celltypeA, "_density")) > median.density.target) %>%
filter(type == type.select) %>%
slice_min(order_by = integral, n = 20) %>% pull(image_id)
# Colour scheme for spatial plots.
col_vec["CD8Tc/DPTc"] <- "red"
col_vec["Tumor"] <- "blue"
# Example spatial plots: dispersed vs. localized T-cell patterns.
(p.disp.tcells <- plotSpatial(sce.sub[, sce.sub$image_id %in% imgids.max], "image_id",
c("Pos_X", "Pos_Y"), node_color_by = cluster.lvl,
node_size_fix = 0.5) +
scale_colour_manual(values = col_vec, name = "Celltypes") +
theme_void() +
theme(strip.text.x = element_text(size = 4), legend.position = "bottom"))
(p.loc.tcells <- plotSpatial(sce.sub[, sce.sub$image_id %in% imgids.min], "image_id",
c("Pos_X", "Pos_Y"), node_color_by = cluster.lvl,
node_size_fix = 0.5) +
scale_colour_manual(values = col_vec, name = "Celltypes") +
theme_void() +
theme(strip.text.x = element_text(size = 4), legend.position = "bottom"))
```
### Assign samples to Lcross category
```{r}
celltypeA <- "CD8Tc/DPTc"
celltypeB <- "Tumor"
cluster.lvl <- "celltype_tcells"
type.select <- "distal"
df <- df.cd8tcells
# Median target-cell (T-cell) density across images.
median.density.target <- median(df[[paste0(celltypeA, "_density")]])
# Median distal integral among images with above-median target-cell density.
median.lcross.value <- median(df$integral[df$type == type.select &
df[paste0(celltypeA, "_density")] >= median.density.target])
# Map images to samples.
sample.ids <- colData(sce) %>% as_tibble() %>% distinct(biosample_id, image_id)
df.samples <- df %>% filter(type == type.select) %>% left_join(sample.ids)
# Per-sample mean density/integral, dichotomized into:
# high = dispersed, low = localized (both above density threshold),
# cold = below density threshold.
df.samples.mean <- df.samples %>% group_by(biosample_id) %>%
summarise(mean.CD8Tc_density = mean(`CD8Tc/DPTc_density`),
mean.integral = mean(integral)) %>%
mutate(lcross.cat.biosample = case_when(
mean.CD8Tc_density >= median.density.target & mean.integral >= median.lcross.value ~ "high",
mean.CD8Tc_density >= median.density.target & mean.integral < median.lcross.value ~ "low",
mean.CD8Tc_density < median.density.target ~ "cold"))
# Same classification at the ROI (image) level, to assess within-sample variability.
df.images <- df.samples %>% group_by(image_id) %>%
mutate(lcross.cat.image = case_when(
`CD8Tc/DPTc_density` >= median.density.target & integral >= median.lcross.value ~ "high",
`CD8Tc/DPTc_density` >= median.density.target & integral < median.lcross.value ~ "low",
`CD8Tc/DPTc_density` < median.density.target ~ "cold"))
# Two samples (NB18-1420, NB17-754) had too few CD8/DP T cells to classify; assign
# them manually as "cold" at both sample and image level.
df.samples.mean <- df.samples.mean %>%
add_row(biosample_id = "NB18-1420", lcross.cat.biosample = "cold") %>%
add_row(biosample_id = "NB17-754", lcross.cat.biosample = "cold")
add.rows <- tibble(image_id = unique(sce$image_id[sce$biosample_id == "NB18-1420"]),
biosample_id = "NB18-1420",
lcross.cat.image = "cold")
add.rows <- add.rows %>% bind_rows(
tibble(image_id = unique(sce$image_id[sce$biosample_id == "NB17-754"]),
biosample_id = "NB17-754",
lcross.cat.image = "cold"))
df.images <- df.images %>% bind_rows(add.rows)
```
#### Suppl plot: spatial assignment ROI vs. Sample
```{r, fig.width= 8, fig.height = 6}
# Compare per-ROI infiltration-pattern assignment against the overall sample
# assignment (how consistent are ROIs within a sample).
df.plot <- df.images %>%
left_join(df.samples.mean) %>%
group_by(biosample_id, lcross.cat.biosample, lcross.cat.image) %>%
summarise(count = n()) %>%
left_join(colData(sce) %>% as_tibble() %>% distinct(biosample_id, pid))
(p.suppl <- df.plot %>%
ggplot(aes(x = reorder(biosample_id, pid))) +
geom_bar(aes(y = count, fill = lcross.cat.image), stat = "identity", position = "fill") +
scale_fill_discrete(name = "Assigned infiltration\npattern per ROI",
labels = c("cold" = "cold", "low" = "localized", "high" = "dispersed")) +
new_scale_fill() +
geom_tile(aes(y = -0.04, height = .04, fill = as.factor(pid))) +
scale_fill_manual(values = metadata(sce)$col_clinical$pid, guide = "none") +
scale_y_continuous(breaks = c(-0.04, 0.00, 0.25, 0.5, 0.75, 1),
labels = c("Patient", "0.00", "0.25", "0.50", "0.75", "1.00")) +
facet_grid(~lcross.cat.biosample, scales = "free_x", space = "free",
labeller = as_labeller(c("cold" = "cold", "low" = "localized", "high" = "dispersed"))) +
ylab("Proportions") + xlab("Sample") +
theme(axis.text.x = element_blank()))
ggsave(plot = p.suppl, filename = "figures/suppfig_spatialAssignement.pdf",
width = 8, height = 4, dpi = 300)
```
#### Suppl plot: spatial assignment vs. ICI responder/non-responder
```{r, fig.height = 10, fig.width = 8}
# Relate the spatial infiltration category to ICI response and pretreatment.
df.cond <- colData(sce) %>% as_tibble() %>%
distinct(biosample_id, cond.2, cond.3) %>% left_join(df.samples.mean)
# Responder vs. non-responder (ICI-naive samples).
p1 <- df.cond %>% filter(!is.na(cond.2) & !is.na(lcross.cat.biosample)) %>%
group_by(cond.2, lcross.cat.biosample) %>%
summarise(count = n()) %>% mutate(prop = prop.table(count)) %>%
ggplot(aes(x = cond.2, y = prop, fill = lcross.cat.biosample)) +
geom_bar(stat = "identity") +
geom_text(aes(label = scales::percent(prop, accuracy = 1, suffix = "%")),
position = position_stack(vjust = 0.5)) +
scale_fill_discrete(name = "Spatial infiltration pattern",
labels = c("cold" = "cold", "low" = "localized", "high" = "dispersed")) +
scale_x_discrete(name = "", labels = c("non-responder" = "Non-responder", "responder" = "Responder")) +
ylab("Proportion") + theme_minimal() + ggtitle("ICI naive with postoperative ICI")
# ICI pretreated vs. naive.
p2 <- df.cond %>% filter(!is.na(cond.3) & !is.na(lcross.cat.biosample)) %>%
group_by(cond.3, lcross.cat.biosample) %>%
summarise(count = n()) %>% mutate(prop = prop.table(count)) %>%
ggplot(aes(x = cond.3, y = prop, fill = lcross.cat.biosample)) +
geom_bar(stat = "identity") +
geom_text(aes(label = scales::percent(prop, accuracy = 1, suffix = "%")),
position = position_stack(vjust = 0.5)) +
scale_fill_discrete(name = "Spatial infiltration pattern",
labels = c("cold" = "cold", "low" = "localized", "high" = "dispersed")) +
scale_x_discrete(name = "", labels = c("ici_naive" = "ICI naive", "ici_pretreated" = "ICI pretreated")) +
ylab("Proportion") + theme_minimal() + ggtitle("ICI pretreated vs. naive")
# Within responders, dichotomize by median OS/PFS and test association with the
# spatial category (Fisher's exact test).
clin <- metadata(sce)$clinical %>%
select(id, biosample_id, days_dos, days_dg, days_dos_pfs, status, status_pfs)
df.clin.cond <- df.cond %>% left_join(clin)
medianOS_cond.2 <- df.clin.cond %>% filter(cond.2 == "responder") %>% summarise(median(days_dos)) %>% pull()
medianPFS_cond.2 <- df.clin.cond %>% filter(cond.2 == "responder") %>% summarise(median(days_dos_pfs)) %>% pull()
df.clin.cond <- df.clin.cond %>% rowwise() %>%
mutate(cond.2.new_OS = case_when(cond.2 == "responder" & days_dos >= medianOS_cond.2 ~ "OS_high",
cond.2 == "responder" & days_dos < medianOS_cond.2 ~ "OS_low"),
cond.2.new_PFS = case_when(cond.2 == "responder" & days_dos_pfs >= medianPFS_cond.2 ~ "PFS_high",
cond.2 == "responder" & days_dos_pfs < medianPFS_cond.2 ~ "PFS_low"))
df.clin.cond %>% filter(cond.2 == "responder") %>% filter(!is.na(lcross.cat.biosample)) %>%
select(cond.2.new_OS, lcross.cat.biosample) %>% ungroup() %>%
summarise(pvalue = fisher.test(cond.2.new_OS, lcross.cat.biosample)$p.value)
df.clin.cond %>% filter(cond.2 == "responder") %>% filter(!is.na(lcross.cat.biosample)) %>%
select(cond.2.new_PFS, lcross.cat.biosample) %>% ungroup() %>%
summarise(pvalue = fisher.test(cond.2.new_PFS, lcross.cat.biosample)$p.value)
# Group sizes.
df.clin.cond %>% filter(cond.2 == "responder") %>% filter(!is.na(lcross.cat.biosample)) %>%
group_by(cond.2.new_OS) %>% summarise(count = n())
df.clin.cond %>% filter(cond.2 == "responder") %>% filter(!is.na(lcross.cat.biosample)) %>%
group_by(cond.2.new_PFS) %>% summarise(count = n())
# Stacked composition of spatial categories by survival cutoff within responders.
p3 <- df.clin.cond %>% filter(cond.2 == "responder") %>% filter(!is.na(lcross.cat.biosample)) %>%
pivot_longer(cols = c("cond.2.new_OS", "cond.2.new_PFS"),
names_to = "cond", names_prefix = "cond.2.new_", values_to = "surv") %>%
rowwise() %>% mutate(surv = str_split(surv, "_")[[1]][2]) %>%
group_by(cond.2, cond, surv, lcross.cat.biosample) %>%
summarise(count = n()) %>% mutate(prop = prop.table(count)) %>%
ggplot(aes(x = surv, y = count, fill = lcross.cat.biosample)) +
geom_bar(stat = "identity") +
geom_text(aes(label = scales::percent(prop, accuracy = 1, suffix = "%")),
position = position_stack(vjust = 0.5)) +
facet_wrap(~ cond) +
scale_fill_discrete(name = "Spatial infiltration pattern",
labels = c("cold" = "cold", "low" = "localized", "high" = "dispersed")) +
scale_x_discrete(name = "Survival cutoff") +
ylab("Proportion") + theme_minimal() + ggtitle("ICI responder")
(p1 | plot_spacer()) / (p2 | plot_spacer()) / p3 +
plot_layout(guides = "collect") + plot_annotation(tag_levels = "A")
ggsave("figures/supplfig_spatialAssignement2.pdf", dpi = 300, width = 8, height = 10)
```
### Example images
```{r}
df <- df.cd8tcells
cluster.lvl <- "celltype_tcells"
# Identify samples per spatial category (dispersed/localized/cold) using the
# mean integral (dispersed = high, localized = low) or density (cold).
biosampleids.dispersed <- df.samples.mean %>% filter(lcross.cat.biosample == "high") %>%
slice_max(order_by = mean.integral, n = 100) %>% pull(biosample_id)
biosampleids.localized <- df.samples.mean %>% filter(lcross.cat.biosample == "low") %>%
slice_min(order_by = mean.integral, n = 100) %>% pull(biosample_id)
biosampleids.cold <- df.samples.mean %>% filter(lcross.cat.biosample == "cold") %>%
slice_min(order_by = mean.CD8Tc_density, n = 100) %>% pull(biosample_id)
# Within those samples, select representative images (above-median target density
# for dispersed/localized; below-median for cold), ranked by integral/density.
imgids.dispersed <- df %>% filter(type == type.select) %>% left_join(sample.ids) %>%
filter(biosample_id %in% biosampleids.dispersed) %>%
filter(!!rlang::sym(paste0(celltypeA, "_density")) > median.density.target) %>%
slice_max(order_by = integral, n = 20) %>%
slice_max(order_by = !!rlang::sym(paste0(celltypeA, "_density")), n = 16) %>% pull(image_id)
imgids.localized <- df %>% filter(type == type.select) %>% left_join(sample.ids) %>%
filter(biosample_id %in% biosampleids.localized) %>%
filter(!!rlang::sym(paste0(celltypeA, "_density")) > median.density.target) %>%
slice_min(order_by = integral, n = 20) %>%
slice_max(order_by = !!rlang::sym(paste0(celltypeA, "_density")), n = 16) %>% pull(image_id)
imgids.cold <- df %>% filter(type == type.select) %>% left_join(sample.ids) %>%
filter(biosample_id %in% biosampleids.cold) %>%
filter(!!rlang::sym(paste0(celltypeA, "_density")) < median.density.target) %>%
slice_max(order_by = !!rlang::sym(paste0(celltypeA, "_density")), n = 16) %>% pull(image_id)
col_vec["CD8Tc/DPTc"] <- "red"
col_vec["Tumor"] <- "blue"
# Spatial example plots for each pattern.
(p.disp.tcells <- plotSpatial(sce.sub[, sce.sub$image_id %in% imgids.dispersed],
"image_id", c("Pos_X", "Pos_Y"),
node_color_by = cluster.lvl, node_size_fix = 0.5, ncol = 4) +
scale_colour_manual(values = col_vec, name = "Celltypes") +
theme_void() +
theme(legend.position = "bottom",
strip.text.x = element_text(size = 4),
strip.background = element_blank(),
strip.text.x.top = element_blank()))
(p.loc.tcells <- plotSpatial(sce.sub[, sce.sub$image_id %in% imgids.localized],
"image_id", c("Pos_X", "Pos_Y"),
node_color_by = cluster.lvl, node_size_fix = 0.5, ncol = 4) +
scale_colour_manual(values = col_vec, name = "Celltypes") +
theme_void() +
theme(legend.position = "bottom",
strip.text.x = element_text(size = 4),
strip.background = element_blank(),
strip.text.x.top = element_blank()))
(p.cold.tcells <- plotSpatial(sce.sub[, sce.sub$image_id %in% imgids.cold],
"image_id", c("Pos_X", "Pos_Y"),
node_color_by = cluster.lvl, node_size_fix = 0.5, ncol = 4) +
scale_colour_manual(values = col_vec, name = "Celltypes") +
theme_void() +
theme(legend.position = "bottom",
strip.text.x = element_text(size = 4),
strip.background = element_blank(),
strip.text.x.top = element_blank()))
```
#### Select survival plot
```{r, fig.height= 4, fig.width = 5}
library(ggsurvfit)
# Combine sample-level spatial category with clinical/survival data.
df.all <- colData(sce) %>% as_tibble() %>% distinct(biosample_id, .keep_all = T) %>%
select(id, biosample_id, contains("cond"))
clin <- metadata(sce)$clinical %>%
select(id, biosample_id, days_dos, days_dg, days_dos_pfs, status, status_pfs)
df.subset_all <- df.all %>% left_join(df.samples.mean) %>% left_join(clin)
# Overall survival by spatial infiltration pattern.
fit <- survfit2(Surv(days_dos / 365.25, status) ~ lcross.cat.biosample, data = df.subset_all)
p.os <- fit %>% ggsurvfit() + add_pvalue() + add_censor_mark() + add_quantile(y_value = 0.5) +
coord_cartesian(xlim = c(0, 8)) +
scale_y_continuous(limits = c(0, 1), labels = scales::percent, expand = c(0.01, 0)) +
scale_color_discrete(labels = c("cold" = "cold", "low" = "localized", "high" = "dispersed")) +
theme_minimal() + theme(text = element_text(size = 15)) +
xlab("Follow-up time, years") + ggtitle("Overall survival")
p.os <- ggsurvfit_build(p.os)
# Median survival per group.
median_table <- fit %>% tidy_survfit() %>% group_by(strata) %>%
summarize(median_survival = time[which(estimate <= 0.5)[1]], .groups = "drop")
print(median_table)
p.os
# Local progression-free survival by spatial infiltration pattern.
fit <- survfit2(Surv(days_dos_pfs / 365.25, status_pfs) ~ lcross.cat.biosample, data = df.subset_all)
p.pfs <- fit %>% ggsurvfit() + add_pvalue() + add_censor_mark() + add_quantile(y_value = 0.5) +
coord_cartesian(xlim = c(0, 8)) +
scale_y_continuous(limits = c(0, 1), labels = scales::percent, expand = c(0.01, 0)) +
scale_color_discrete(labels = c("cold" = "cold", "low" = "localized", "high" = "dispersed")) +
theme_minimal() + theme(text = element_text(size = 15)) +
xlab("Follow-up time, years") + ggtitle("Local progression-free survival")
p.pfs <- ggsurvfit_build(p.pfs)
median_table <- fit %>% tidy_survfit() %>% group_by(strata) %>%
summarize(median_survival = time[which(estimate <= 0.5)[1]], .groups = "drop")
print(median_table)
p.pfs
```
### Plot cell abundances for Lcross high/low
```{r, fig.width = 12, fig.height=10}
coldata <- colData(sce) %>% as_tibble()
cluster.lvl <- "cl_all_spec"
# Within the target T-cell population, compute per-sample subset proportions for
# dispersed ("high") vs. localized ("low") samples (cold excluded).
df <- coldata %>% left_join(df.subset_all %>% select(biosample_id, lcross.cat.biosample)) %>%
filter(biosample_id %in% df.subset_all$biosample_id[!is.na(df.subset_all$lcross.cat.biosample)] &
celltype_tcells == celltypeA) %>%
filter(lcross.cat.biosample != "cold") %>%
group_by(biosample_id, lcross.cat.biosample, cl_all_spec) %>%
summarise(count = n()) %>% mutate(pct = prop.table(count))
# Wilcoxon test per cell type (FDR-adjusted), with formatted p-value labels.
stat.test <- df %>%
group_by(across(cluster.lvl)) %>%
wilcox_test(as.formula(paste("pct", "~", "lcross.cat.biosample"))) %>%
adjust_pvalue(p.col = "p", method = "fdr", output.col = "p.adj") %>%
add_significance(p.col = "p.adj",
cutpoints = c(0, 1e-04, 0.001, 0.01, 0.1, 1),
output.col = "p.adj.signif") %>%
add_xy_position(scales = "free", x = "lcross.cat.biosample") %>%
mutate(x.position = 1.5) %>%
rowwise() %>%
mutate(p.adj.format = p_format(p.adj, add.p = T, leading.zero = F, space = F),
p.format = p_format(p, add.p = T, leading.zero = F, space = F),
p.adj.format = str_replace(p.adj.format, "p=", "p.adj="))
df <- df %>% left_join(stat.test) %>%
rowwise() %>%
mutate(signif = if_else(p.adj.signif != "ns" & p < 0.05, "sign", "ns"),
p.adj.signif.plot = if_else(signif == "ns", "ns", p.adj.signif))
(p.wilcox.lcrossprop <- ggplot(df, aes(x = lcross.cat.biosample, y = pct, fill = cl_all_spec)) +
geom_boxplot(aes(alpha = signif), outlier.shape = NA) +
geom_quasirandom(alpha = 0.35, dodge.width = 0.75) +
scale_alpha_manual(values = c(1, 1), guide = "none") +
geom_text(aes(label = p.adj.signif.plot, x = x.position, y = y.position, vjust = 0.8), size = 3) +
facet_wrap(as.formula(paste("~", cluster.lvl, "+ p.format + p.adj.format")), scales = "free") +
scale_fill_manual(name = "Celltypes", values = col_ann_hm) +
ylab("Proportions") + xlab("") +
scale_x_discrete(guide = guide_axis(n.dodge = 2),
labels = c("high" = "dispersed", "low" = "localized")) +
theme(text = element_text(size = 15)))
```
## Export single panels
```{r, fig.width=5, fig.height=5}
# Figure 3, panel row 1: Cox proportions + example spatial patterns.
(p.prop.os | p.disp.tcells + ggtitle("Dispersed CD8/DP T cells") & theme(legend.position = "none") |
p.loc.tcells + ggtitle("Localized CD8/DP T cells") |
p.cold.tcells + ggtitle("Cold CD8/DP T cells") & theme(legend.position = "none")) +
plot_layout(widths = c(0.5, 1, 1, 1)) +
plot_annotation(tag_levels = list(c("A", "", "", "B", "C", "D")))
ggsave("figures/fig3_1.pdf", dpi = 200, width = 18, height = 6)
# Survival panels.
p.os + plot_annotation(tag_levels = list(c("D")))
ggsave("figures/fig3_2_os.pdf", dpi = 300, width = 5, height = 4)
p.pfs + plot_annotation(tag_levels = list(c("E")))
ggsave("figures/fig3_2_pfs.pdf", dpi = 300, width = 5, height = 4)
p.wilcox.lcrossprop + theme(legend.position = "none") + plot_annotation(tag_levels = list(c("F")))
ggsave("figures/fig3_2_wilcox.pdf", dpi = 300, width = 5, height = 5)
# Figure 3, panel row 3: cellular-neighbourhood boxplot + magnified spatial
# examples. NOTE: this depends on `output.plot` and the `p.min.*` / `p.max.*`
# panels created later in the CN ggmagnify section, so run that section first.
(output.plot | ((p.min.1 | p.min.2) / (p.max.1 | p.max.2) +
plot_layout(guides = "collect") &
theme(legend.position = "bottom", plot.margin = margin(15, 18, 0, 0, "mm")))) +
plot_layout(widths = c(2, 1))
ggsave("figures/fig3_3.eps", dpi = 200, width = 18, height = 7)
p.magnify <- (p.min.1 + theme(legend.position = "none", plot.margin = margin(13, 16, 0, 0, "mm")) |
p.min.2 + theme(legend.position = "none", plot.margin = margin(0, 12, 0, 0, "mm"))) /
(p.max.1 + theme(legend.position = "none", plot.margin = margin(0, 15, 0, 0, "mm")) |
p.max.2 + theme(legend.position = "none", plot.margin = margin(0, 0, 0, 0, "mm")))
p.magnify
ggsave("figures/fig3_3_magnify.pdf", dpi = 600, width = 5, height = 5)
```
# Cellular neighbourhood (CN) analysis
## 1. Aggregate neighbors
```{r}
# Aggregate each cell's neighbourhood into cell-type composition vectors
# (counts of cl_all_spec among neighbours defined by the "neighborhood" graph).
sce <- aggregateNeighbors(sce, colPairName = "neighborhood",
aggregate_by = "metadata", count_by = "cl_all_spec",
name = "aggregatedNeighbors_cl_all_spec")
```
## 2. Optimal cluster number
```{r}
# Estimate a reasonable number of neighbourhood clusters (k) on a 1% subsample,
# using the elbow and silhouette criteria.
cn <- as.data.frame(sce$aggregatedNeighbors_cl_all_spec)
cn_subset <- cn[sample(nrow(cn), 0.01 * nrow(cn)), ]
set.seed(12345)
library(factoextra)
# Elbow method (within-cluster sum of squares).
fviz_nbclust(cn_subset, kmeans, k.max = 25, method = "wss") +
labs(subtitle = "cl_all_spec - Elbow method")
# Silhouette method.
fviz_nbclust(cn_subset, kmeans, k.max = 25, method = "silhouette") +
labs(subtitle = "cl_all_spec - Silhouette method")
```
## 3. Define SCE subset
Exclude unassigned cells for CN detection/analysis.
```{r}
sce.subset <- sce
sce.subset <- sce.subset[, sce.subset$cl_all_spec != "unassigned"]
```
## 4. CN plotting
### 4.1. Define k and subset
```{r}
k <- 7
cond.group <- "lcross.cat.biosample"
CN <- "aggregatedNeighbors_cl_all_spec"
celltypeCluster <- "cl_all_spec"
# Biological interpretation of each CN cluster.
cn.labels <- c("cn_1" = "Perivascular",
"cn_2" = "Tumor core",
"cn_3" = "Tumor/CD8Tc\nborder",
"cn_4" = "Myeloid enriched",
"cn_5" = "CD8Tc core",
"cn_6" = "B cell, Neutrophil,\nMf enriched",
"cn_7" = "CD4Tc/Treg,\nPlasma cell enriched")
# Re-create the subset (exclude unassigned) and run k-means on neighbourhood
# composition vectors.
sce.subset <- sce[, sce$cl_all_spec != "unassigned"]
set.seed(1234)
cn_1 <- kmeans(sce.subset[[CN]], centers = k)
sce.subset$cn_celltypes <- as.factor(paste0("cn_", cn_1$cluster))
```
### 4.2. Heatmap + boxplot
```{r, fig.width=12}
color.palette <- colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(100)
# Cell-type composition per CN (row-normalized), as a clustered heatmap.
for_plot <- prop.table(table(sce.subset$cn_celltypes, sce.subset[[celltypeCluster]]), margin = 1)
hm.matrix <- scale(for_plot, center = T, scale = T)
row.names(hm.matrix) <- as.vector(cn.labels)
hm <- pheatmap(hm.matrix, cluster_rows = T, color = color.palette,
legend = T, annotation_legend = F)
row.order <- row_order(draw(hm))
# CN frequency per sample, with the per-CN cohort median.
df.cn.sum <- colData(sce.subset) %>% as_tibble() %>%
group_by(biosample_id, cn_celltypes, .drop = FALSE) %>%
summarise(count = n()) %>% mutate(prop = prop.table(count)) %>%
group_by(cn_celltypes) %>% mutate(median = median(prop))
# Group samples high/low per CN relative to the median.
df.cn.group <- df.cn.sum %>% distinct(biosample_id, cn_celltypes, .keep_all = T) %>%
rowwise() %>% mutate(cn_group = if_else(prop < median, "low", "high")) %>%
left_join(clin)
# Combine CN frequencies with the spatial infiltration category.
df.comb <- df.cn.sum %>%
left_join(df.samples.mean %>% select(biosample_id, all_of(cond.group))) %>%
filter(!is.na(!!rlang::sym(cond.group)))
df.comb$cn_celltypes <- factor(df.comb$cn_celltypes, levels = rev(rownames(for_plot)[row.order]))
# Boxplot of CN proportions across infiltration categories.
boxplot <- ggboxplot(df.comb, x = "cn_celltypes", y = "prop",
fill = "lcross.cat.biosample", outlier.shape = NA)
stat.test <- df.comb %>% group_by(cn_celltypes) %>%
wilcox_test(as.formula(paste("prop", "~", cond.group))) %>%
adjust_pvalue(p.col = "p", method = "fdr", output.col = "p.adj") %>%
add_significance(p.col = "p.adj", cutpoints = c(0, 1e-04, 0.001, 0.01, 0.1, 1),
output.col = "p.adj.signif") %>%
add_xy_position(x = "cn_celltypes", dodge = 0.8)
p.boxplot <- boxplot +
stat_pvalue_manual(stat.test, label = "p.adj.signif", bracket.nudge.y = -0.1,
tip.length = 0.01, hide.ns = T, coord.flip = T, size = 8) +
coord_flip() +
scale_fill_discrete(labels = c(cold = "cold", high = "dispersed", low = "localized"), name = "") +
ylab("Proportion") + xlab("") +
theme(axis.text.y = element_blank())
(output.plot <- wrap_plots(grid.grabExpr(draw(hm)),
p.boxplot & theme(legend.position = "bottom"), widths = c(7, 2)))
```
### 4.3. Spatial plotting
```{r}
# CN colour scheme.
col_ann <- setNames(dittoColors()[seq_along(unique(sce.subset$cn_celltypes))],
unique(sce.subset$cn_celltypes))
# Example CN spatial maps for dispersed (imgids.max) and localized (imgids.min) ROIs.
plotSpatial(sce.subset[, sce.subset$sample_id %in% imgids.max],
node_color_by = "cn_celltypes", img_id = "image_id", node_size_fix = 0.4) +
scale_color_manual(values = col_ann) +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(),
aspect.ratio = 1, strip.text.x = element_blank())
plotSpatial(sce.subset[, sce.subset$sample_id %in% imgids.min],
node_color_by = "cn_celltypes", img_id = "image_id", node_size_fix = 0.4) +
scale_color_manual(values = col_ann) +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(),
aspect.ratio = 1, strip.text.x = element_blank())
```
### 4.4. ggmagnify
```{r}
# CN to highlight in the magnified examples.
cn.select <- c("cn_3")
df.coldata <- colData(sce.subset) %>% as_tibble()
# Overview of selected example images, highlighting the chosen CN.
image.select <- c("NB16-121_005", "NB18-818_007", "NB21-437_011", "NB21-42_010", "NB21-42_002")
ggplot(df.coldata %>% filter(image_id %in% image.select),
aes(x = Pos_X, y = Pos_Y, alpha = cn_celltypes %in% cn.select)) +
geom_point(aes(fill = cn_celltypes, colour = cn_celltypes), shape = 21, size = 1.) +
scale_colour_manual(values = col_ann) +
scale_fill_manual(values = col_ann) +
scale_alpha_discrete(range = c(0.25, 1), guide = "none") +
facet_wrap(~ image_id) +
theme_void() +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(), aspect.ratio = 1) +
guides(color = guide_legend(override.aes = list(size = 4)))
library(ggmagnify)
# Two "dispersed" examples with a magnified inset over the highlighted CN.
p.max <- ggplot(df.coldata %>% filter(image_id %in% c("NB14-599_004", "NB16-121_011")),
aes(x = Pos_X, y = Pos_Y, alpha = cn_celltypes %in% cn.select)) +
geom_point(aes(color = cn_celltypes), shape = 19, size = 1.8, stroke = 0) +
scale_color_manual(values = col_ann) +
scale_alpha_discrete(range = c(1, 1), guide = "none") +
facet_wrap(~ image_id) +
theme_void() +
coord_cartesian(clip = "off") +
theme(panel.background = element_rect(fill = 'white', color = "transparent"),
panel.spacing.x = unit(5, "lines"), panel.spacing.y = unit(4, "lines"),
axis.text.x = element_blank(), axis.text.y = element_blank(),
aspect.ratio = 1, plot.margin = ggplot2::margin(100, 0, 0, 0)) +
guides(color = guide_legend(override.aes = list(size = 4))) +
geom_magnify(aes(from = cn_celltypes %in% cn.select &
Pos_Y > 150 & Pos_Y < 400 & Pos_X > 150 & Pos_X < 350),
to = c(500, 800, 500, 800)) +
scale_alpha_discrete(range = c(0.2, 1), guide = "none")
# Same examples as individual panels (for flexible figure assembly).
(p.max.1 <- ggplot(df.coldata %>% filter(image_id %in% c("NB14-599_004")),
aes(x = Pos_X, y = Pos_Y, alpha = cn_celltypes %in% cn.select)) +
geom_point(aes(color = cn_celltypes), shape = 19, size = 1.8, stroke = 0, show.legend = F) +
scale_color_manual(values = col_ann, labels = cn.labels, name = "CN") +
scale_alpha_discrete(range = c(1, 1), guide = "none") +
theme_void() + coord_cartesian(clip = "off") +
theme(panel.background = element_rect(fill = 'white', color = "transparent"),
panel.spacing.x = unit(5, "lines"), panel.spacing.y = unit(4, "lines"),
axis.text.x = element_blank(), axis.text.y = element_blank(),
aspect.ratio = 1, plot.margin = ggplot2::margin(100, 0, 0, 0)) +
guides(color = guide_legend(override.aes = list(size = 4))) +
geom_magnify(from = c(250, 400, 100, 250), to = c(500, 800, 500, 800)) +
scale_alpha_discrete(range = c(0.2, 1), guide = "none") +
theme(legend.position = "none"))
(p.max.2 <- ggplot(df.coldata %>% filter(image_id %in% c("NB16-121_011")),
aes(x = Pos_X, y = Pos_Y, alpha = cn_celltypes %in% cn.select)) +
geom_point(aes(color = cn_celltypes), shape = 19, size = 1.8, stroke = 0, show.legend = F) +
scale_color_manual(values = col_ann, labels = cn.labels, name = "CN") +
scale_alpha_discrete(range = c(1, 1), guide = "none") +
theme_void() + coord_cartesian(clip = "off") +
theme(panel.background = element_rect(fill = 'white', color = "transparent"),
panel.spacing.x = unit(5, "lines"), panel.spacing.y = unit(4, "lines"),
axis.text.x = element_blank(), axis.text.y = element_blank(),
aspect.ratio = 1, plot.margin = ggplot2::margin(100, 0, 0, 0)) +
guides(color = guide_legend(override.aes = list(size = 4))) +
geom_magnify(from = c(150, 300, 200, 350), to = c(500, 800, 500, 800)) +
scale_alpha_discrete(range = c(0.2, 1), guide = "none") +
theme(legend.position = "none"))
# Two "localized" examples with magnified insets.
p.min <- ggplot(df.coldata %>% filter(image_id %in% c("NB20-1792_006", "NB14-151_008")),
aes(x = Pos_X, y = Pos_Y, alpha = cn_celltypes %in% cn.select)) +
geom_point(aes(color = cn_celltypes), shape = 19, size = 1.8, stroke = 0) +
scale_color_manual(values = col_ann) +
scale_alpha_discrete(range = c(1, 1), guide = "none") +
facet_wrap(~ image_id) +
theme_void() + coord_cartesian(clip = "off") +
theme(panel.background = element_rect(fill = 'white', color = "transparent"),
panel.spacing.x = unit(5, "lines"), panel.spacing.y = unit(4, "lines"),
axis.text.x = element_blank(), axis.text.y = element_blank(),
aspect.ratio = 1, plot.margin = ggplot2::margin(100, 0, 0, 0)) +
guides(color = guide_legend(override.aes = list(size = 4))) +
geom_magnify(aes(from = cn_celltypes %in% cn.select &
Pos_Y > 150 & Pos_Y < 400 & Pos_X > 150 & Pos_X < 350),
to = c(500, 800, 500, 800)) +
scale_alpha_discrete(range = c(0.2, 1), guide = "none")
(p.min.1 <- ggplot(df.coldata %>% filter(image_id %in% c("NB20-1792_006")),
aes(x = Pos_X, y = Pos_Y, alpha = cn_celltypes %in% cn.select)) +
geom_point(aes(color = cn_celltypes), shape = 19, size = 1.8, stroke = 0, show.legend = F) +
scale_color_manual(values = col_ann, labels = cn.labels, name = "CN") +
scale_alpha_discrete(range = c(1, 1), guide = "none") +
theme_void() + coord_cartesian(clip = "off") +
theme(panel.background = element_rect(fill = 'white', color = "transparent"),
panel.spacing.x = unit(5, "lines"), panel.spacing.y = unit(4, "lines"),
axis.text.x = element_blank(), axis.text.y = element_blank(),
aspect.ratio = 1, plot.margin = ggplot2::margin(100, 0, 0, 0)) +
guides(color = guide_legend(override.aes = list(size = 4))) +
geom_magnify(from = c(250, 400, 250, 400), to = c(500, 800, 500, 800)) +
scale_alpha_discrete(range = c(0.2, 1), guide = "none") +
theme(legend.position = "none"))
(p.min.2 <- ggplot(df.coldata %>% filter(image_id %in% c("NB14-151_008")),
aes(x = Pos_X, y = Pos_Y, alpha = cn_celltypes %in% cn.select)) +
geom_point(aes(color = cn_celltypes), shape = 19, size = 1.8, stroke = 0) +
scale_color_manual(values = col_ann, labels = cn.labels, name = "CN") +
scale_alpha_discrete(range = c(1, 1), guide = "none") +
theme_void() + coord_cartesian(clip = "off") +
theme(panel.background = element_rect(fill = 'white', color = "transparent"),
panel.spacing.x = unit(5, "lines"), panel.spacing.y = unit(4, "lines"),
axis.text.x = element_blank(), axis.text.y = element_blank(),
aspect.ratio = 1, plot.margin = ggplot2::margin(100, 0, 0, 0)) +
guides(color = guide_legend(override.aes = list(size = 4))) +
geom_magnify(from = c(150, 300, 200, 350), to = c(500, 800, 500, 800)) +
scale_alpha_discrete(range = c(0.2, 1), guide = "none"))
# Combined dispersed/localized magnified examples.
p.max / p.min + plot_layout(guides = "collect") & theme(legend.position = "bottom")
```
# --------------------------------
# FIGURE 4
Clinical groups, e.g. ICI responder vs. non-responder.
## Survival for clinical groups
```{r}
library(ggsurvfit)
df.all <- colData(sce) %>% as_tibble() %>% distinct(biosample_id, .keep_all = T) %>%
select(id, biosample_id, contains("cond"))
clin <- metadata(sce)$clinical %>%
select(id, biosample_id, days_dos, days_dg, days_dos_pfs, status, status_pfs)
df.subset_all <- df.all %>% left_join(clin)
# Overall survival by ICI response (cond.2).
fit <- survfit2(Surv(days_dos / 365.25, status) ~ cond.2, data = df.subset_all)
p.os <- fit %>% ggsurvfit() + add_pvalue() + add_censor_mark() + add_quantile(y_value = 0.5) +
coord_cartesian(xlim = c(0, 8)) +
scale_y_continuous(limits = c(0, 1), labels = scales::percent, expand = c(0.01, 0)) +
theme_minimal() + xlab("Follow-up time, years") +
ggtitle("ICI: Overall survival") + theme(text = element_text(size = 15))
p.os <- ggsurvfit_build(p.os)
median_table <- fit %>% tidy_survfit() %>% group_by(strata) %>%
summarize(median_survival = time[which(estimate <= 0.5)[1]], .groups = "drop")
print(median_table)
p.os
# Progression-free survival by ICI response (cond.2).
fit <- survfit2(Surv(days_dos_pfs / 365.25, status_pfs) ~ cond.2, data = df.subset_all)
p.pfs <- fit %>% ggsurvfit() + add_pvalue() + add_censor_mark() + add_quantile(y_value = 0.5) +
coord_cartesian(xlim = c(0, 8)) +
scale_y_continuous(limits = c(0, 1), labels = scales::percent, expand = c(0.01, 0)) +
theme_minimal() + xlab("Follow-up time, years") +
ggtitle("ICI: Progression-free survival") + theme(text = element_text(size = 15))
p.pfs <- ggsurvfit_build(p.pfs)
median_table <- fit %>% tidy_survfit() %>% group_by(strata) %>%
summarize(median_survival = time[which(estimate <= 0.5)[1]], .groups = "drop")
print(median_table)
p.pfs
p.os + p.pfs + plot_layout(ncol = 1, nrow = 4)
```
## Differential abundance
### Wilcoxon
#### cond.2 responder vs. non-responder
```{r, fig.width= 11, fig.height=9}
conditions <- c("cond.2", "cond.3")
cluster.lvl <- "cl_all_spec"
coldata <- colData(sce) %>% as_tibble()
cond <- "cond.2"
# Per-sample cell-type proportions with the ICI-response condition attached.
df <- coldata %>%
filter(!is.na(cond.2)) %>%
group_by(id, across(cluster.lvl)) %>%
dplyr::summarise(count = n()) %>% mutate(pct = prop.table(count)) %>%
filter(!!rlang::sym(cluster.lvl) != 0 & !is.na(!!rlang::sym(cluster.lvl))) %>%
ungroup() %>%
left_join(coldata %>% distinct(id, .keep_all = T) %>% select(id, all_of(cond)), by = "id")
# Wilcoxon test per cell type (FDR-adjusted), with formatted labels.
stat.test <- df %>%
group_by(across(cluster.lvl)) %>%
wilcox_test(as.formula(paste("pct", "~", cond))) %>%
adjust_pvalue(p.col = "p", method = "fdr", output.col = "p.adj") %>%
add_significance(p.col = "p.adj", cutpoints = c(0, 1e-04, 0.001, 0.01, 0.1, 1),
output.col = "p.adj.signif") %>%
add_xy_position(scales = "free", x = cond) %>%
mutate(x.position = 1.5) %>%
rowwise() %>%
mutate(p.adj.format = p_format(p.adj, add.p = T, leading.zero = F, space = F),
p.format = p_format(p, add.p = T, leading.zero = F, space = F),
p.adj.format = str_replace(p.adj.format, "p=", "p.adj="))
df <- df %>% left_join(stat.test) %>%
rowwise() %>%
mutate(signif = if_else(p.adj.signif != "ns" & p < 0.05, "sign", "ns"),
p.adj.signif.plot = if_else(signif == "ns", "", p.adj.signif))
(p.wilcox <- ggplot(df, aes(x = !!rlang::sym(cond), y = pct)) +
geom_boxplot(aes(alpha = signif), outlier.shape = NA) +
geom_quasirandom(aes(alpha = signif), varwidth = T) +
scale_alpha_manual(values = c(0.1, 1)) +
geom_text(aes(label = p.adj.signif.plot, x = x.position, y = y.position, vjust = 0.8), size = 4) +
facet_wrap(as.formula(paste("~", cluster.lvl, "+ p.format + p.adj.format")), scales = "free") +
ylab("Proportions") + xlab("") +
scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
ggtitle(paste(cluster.lvl, ": ", unique(df[cond][!is.na(df[cond])])[1], " vs. ",
unique(df[cond][!is.na(df[cond])])[2])))
```
#### cond.3 ICI naive vs. pretreated
```{r, fig.width= 11, fig.height=9}
cluster.lvl <- "cl_all_spec"
coldata <- colData(sce) %>% as_tibble()
cond <- "cond.3"
# Per-sample cell-type proportions with the ICI-pretreatment condition attached.
df <- coldata %>%
filter(!is.na(cond.3)) %>%
group_by(id, across(cluster.lvl)) %>%
dplyr::summarise(count = n()) %>% mutate(pct = prop.table(count)) %>%
filter(!!rlang::sym(cluster.lvl) != 0 & !is.na(!!rlang::sym(cluster.lvl))) %>%
ungroup() %>%
left_join(coldata %>% distinct(id, .keep_all = T) %>% select(id, all_of(cond)), by = "id")
stat.test <- df %>%
group_by(across(cluster.lvl)) %>%
wilcox_test(as.formula(paste("pct", "~", cond))) %>%
adjust_pvalue(p.col = "p", method = "fdr", output.col = "p.adj") %>%
add_significance(p.col = "p.adj", cutpoints = c(0, 1e-04, 0.001, 0.01, 0.1, 1),
output.col = "p.adj.signif") %>%
add_xy_position(scales = "free", x = cond) %>%
mutate(x.position = 1.5) %>%
rowwise() %>%
mutate(p.adj.format = p_format(p.adj, add.p = T, leading.zero = F, space = F),
p.format = p_format(p, add.p = T, leading.zero = F, space = F),
p.adj.format = str_replace(p.adj.format, "p=", "p.adj="))
df <- df %>% left_join(stat.test) %>%
rowwise() %>%
mutate(signif = if_else(p.adj.signif != "ns" & p < 0.05, "sign", "ns"),
p.adj.signif.plot = if_else(signif == "ns", "", p.adj.signif))
(p.wilcox <- ggplot(df, aes(x = !!rlang::sym(cond), y = pct)) +
geom_boxplot(aes(alpha = signif), outlier.shape = NA) +
geom_quasirandom(aes(alpha = signif), varwidth = T) +
scale_alpha_manual(values = c(0.1, 1)) +
geom_text(aes(label = p.adj.signif.plot, x = x.position, y = y.position, vjust = 0.8), size = 4) +
facet_wrap(as.formula(paste("~", cluster.lvl, "+ p.format + p.adj.format")), scales = "free") +
ylab("Proportions") + xlab("") +
scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
ggtitle(paste(cluster.lvl, ": ", unique(df[cond][!is.na(df[cond])])[1], " vs. ",
unique(df[cond][!is.na(df[cond])])[2])))
```
### edgeR (differential abundance)
Differential-abundance testing of cell-type proportions between clinical groups,
following the OSCA workflow
(https://bioconductor.org/books/3.18/OSCA.multisample/differential-abundance.html).
```{r}
dat <- coldata
dat$id <- as.numeric(dat$id)
# New condition cond.5: ICI non-responder vs. ICI pretreated.
dat$cond.5[dat$cond.2 == "non-responder"] <- "non-responder"
dat$cond.5[dat$cond.3 == "ici_pretreated" & is.na(dat$cond.5)] <- "ici_pretreated"
# Set factor reference levels (first level = reference).
dat$cond.2 <- factor(dat$cond.2, levels = c("non-responder", "responder"))
dat$cond.3 <- factor(dat$cond.3, levels = c("ici_naive", "ici_pretreated"))
dat$cond.4 <- factor(dat$cond.4, levels = c("rt_naive", "rt_bm_preop"))
dat$cond.5 <- factor(dat$cond.5, levels = c("non-responder", "ici_pretreated"))
conditions <- c("cond.2", "cond.3", "cond.4", "cond.5")
# Drop cells with missing cluster labels.
dat <- dat %>% filter(!is.na(!!rlang::sym(cluster.lvl)))
plot.list <- list()
for (cond in conditions) {
dat.l <- dat # fresh copy per iteration
dat.l$group_id <- dat.l[[cond]]
dat.l <- dat.l %>% drop_na(paste(cond))
dat.l$cluster.lvl <- dat.l[[cluster.lvl]]
dat.l$cond <- dat.l[[cond]]
# Cell-type x sample abundance table.
abundances <- table(dat.l$cluster.lvl, dat.l$id)
abundances <- unclass(abundances)
# Attach sample metadata and build the DGEList.
extra.info <- dat.l[match(colnames(abundances), dat.l$id), ]
y.ab <- DGEList(abundances, samples = extra.info)
# Filter low-abundance cell types and fit a quasi-likelihood model.
keep <- filterByExpr(y.ab, group = y.ab$samples$cond)
y.ab <- y.ab[keep, ]
design <- model.matrix(~factor(cond), y.ab$samples)
y.ab <- estimateDisp(y.ab, design, trend = "none")
fit.ab <- glmQLFit(y.ab, design, robust = TRUE, abundance.trend = FALSE)
res <- glmQLFTest(fit.ab, coef = ncol(design))
top <- topTags(res, adjust.method = "BH", n = Inf)
ref.group <- str_split(top$comparison, "1*cond\\)")[[1]][2] # reference group
t.op <- as.data.frame(top)
t.op <- tibble::rownames_to_column(t.op, "Cluster")
# logFC bar plot.
plot.list[["bar"]][[cond]] <-
ggplot(t.op, aes(reorder(Cluster, logFC), logFC)) +
geom_col(aes(fill = FDR < 0.1)) +
coord_flip() +
labs(x = "Cell type", y = "logFC",
title = paste(cluster.lvl, ": ", unique(dat[cond][!is.na(dat[cond])])[1], " vs. ",
unique(dat[cond][!is.na(dat[cond])])[2], " / ref.group: ", ref.group)) +
theme_bw() +
theme(strip.background = element_blank(),
panel.background = element_rect(fill = 'white', colour = 'black'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.y = element_text(size = 10))
# Volcano plot.
plot.list[["scatter"]][[cond]] <-
plot.volcano(t.op, col_ann = NULL) +
labs(title = paste(cluster.lvl, ": ", unique(dat[cond][!is.na(dat[cond])])[1], " vs. ",
unique(dat[cond][!is.na(dat[cond])])[2], " / ref.group: ", ref.group))
}
# Collect results for the "pretreated" comparisons (cond.3, cond.4, cond.5).
data.pretreated <- plot.list$bar$cond.3$data %>%
add_column(condition = "cond.3") %>%
add_column(ref.group = str_split(plot.list$bar$cond.3$labels$title, "ref.group: ")[[1]][2])
data.pretreated <- rbind(data.pretreated,
plot.list$bar$cond.4$data %>% add_column(condition = "cond.4") %>%
add_column(ref.group = str_split(plot.list$bar$cond.4$labels$title, "ref.group: ")[[1]][2]))
data.pretreated <- rbind(data.pretreated,
plot.list$bar$cond.5$data %>% add_column(condition = "cond.5") %>%
add_column(ref.group = str_split(plot.list$bar$cond.5$labels$title, "ref.group: ")[[1]][2]))
# Collect results for the ICI-naive responder comparison (cond.2).
data.naive <- plot.list$bar$cond.2$data %>%
add_column(condition = "cond.2") %>%
add_column(ref.group = str_split(plot.list$bar$cond.2$labels$title, "ref.group: ")[[1]][2])
```
#### Plot
```{r, fig.width=12, fig.height = 3.5}
data.full <- rbind(data.pretreated, data.naive)
# Define the contrasting group label for each comparison.
data.full <- data.full %>%
mutate(alt.group = case_when(
ref.group == "ici_pretreated" & condition == "cond.3" ~ "ici_naive",
ref.group == "ici_pretreated" & condition == "cond.5" ~ "non-responder",
ref.group == "rt_bm_preop" ~ "rt_naive",
ref.group == "responder" ~ "non-responder"))
# Plot only the ICI-naive responder comparison in this figure.
data.subset <- data.full %>% filter(condition %in% c("cond.2"))
labels.conditions <- c("cond.3" = "naive vs. PD-while/-after ICI",
"cond.2" = "ICI-naive: responder vs. non-responder",
"cond.4" = "RT-naive vs. RT-preop",
"cond.5" = "ICI non-responder vs. ICI pretreated")
# Dot plot: point size = p-value, fill = logFC, asterisk = FDR < 0.1.
(p.da.1 <- ggplot(data = data.subset) +
geom_point(data = data.subset %>% filter(PValue >= 0.05),
aes(x = Cluster, y = condition, size = PValue, fill = logFC),
alpha = 0.75, shape = 21, stroke = 2) +
geom_point(data = data.subset %>% filter(PValue < 0.05),
aes(x = Cluster, y = condition, size = PValue, fill = logFC),
alpha = 1, shape = 21, stroke = 2) +
geom_point(data = data.subset %>% filter(FDR < 0.1),
aes(x = Cluster, y = condition), shape = 8) +
scale_size(range = c(10, 1), breaks = c(0.01, 0.05, 0.1)) +
scale_colour_manual(values = c("transparent", "black")) +
scale_fill_gradient2(low = muted("blue"), mid = "white", high = muted("red"),
breaks = c(-3, -2, 0, 2, 4), limits = c(-3, 4),
labels = c("Group 1", 2, "0", 2, "Group 2")) +
scale_y_discrete(labels = labels.conditions) +
geom_point(aes(x = Cluster, y = condition, shape = " "), alpha = 0) + # dummy for legend
scale_shape_manual(name = "* FDR < 0.1", values = 16,
guide = guide_legend(override.aes = list(alpha = 0))) +
xlab("") + ylab("Clinical groups") +
theme_bw() +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(),
axis.text.y = element_blank(), panel.border = element_blank(),
text = element_text(size = 15)))
# Cell-type colour key (x-axis tiles).
p.x.tiles <- ggplot(data = data.subset) +
geom_tile(aes(x = Cluster, y = 0, fill = Cluster)) +
scale_fill_manual(values = col_ann_hm, guide = "none") +
xlab("Cell types") +
theme_void() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
axis.text.y = element_blank(), axis.title.x = element_text(),
panel.border = element_blank(), text = element_text(size = 15))
# Group-label tiles (which clinical group is "Group 1" vs "Group 2").
p.da.2 <- data.subset %>% distinct(ref.group, alt.group, condition) %>%
pivot_longer(cols = c(ref.group, alt.group)) %>%
mutate(value = case_when(value == "ici_pretreated" ~ "ICI pretreated",
value == "ici_pretreated2" ~ "ICI pretreated",
value == "ici_naive" ~ "ICI naive",
value == "rt_bm_preop" ~ "RT preoperative",
value == "rt_naive" ~ "RT naive",
value == "responder" ~ "ICI responder",
value == "non-responder" ~ "ICI non-responder")) %>%
ggplot() +
geom_label(aes(x = name, y = condition, label = value, fill = name),
color = "white", hjust = "left") +
scale_fill_manual(values = c("ref.group" = muted("red"), "alt.group" = muted("blue"))) +
scale_x_discrete(name = "Clinical groups", labels = c("Group 1", "Group 2")) +
theme(axis.text.y = element_blank(), axis.ticks = element_blank(),
axis.title.y = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
legend.position = "none", axis.text.x = element_text(hjust = 0),
text = element_text(size = 15))
p.da.1 + p.da.2 + plot_layout(widths = c(3, 1))
p.da.1.comb <- p.da.1 / p.x.tiles + plot_layout(heights = c(5, .75), ncol = 1, nrow = 2)
p.da.1.comb
p.da.2.comb <- p.da.2 / plot_spacer() + plot_layout(heights = c(5, .1), ncol = 1, nrow = 2)
p.da.2.comb
wrap_plots(p.da.1.comb, p.da.2.comb) + plot_layout(widths = c(2, 1))
```
## SpicyR
```{r}
library(spicyR)
# Clinical data; dichotomize tumor samples by median OS/PFS within each response
# group into "high" vs. "low".
df.clin <- metadata(sce)$clinical
median_os <- df.clin %>% filter(!is.na(cond.2)) %>% group_by(cond.2) %>% summarise(median = median(days_dos))
median_pfs <- df.clin %>% filter(!is.na(cond.2)) %>% group_by(cond.2) %>% summarise(median = median(days_dos_pfs))
survival.cond.2.dich <- df.clin %>% filter(!is.na(cond.2)) %>%
distinct(biosample_id, days_dos, days_dos_pfs, cond.2) %>% rowwise() %>%
mutate(os_dich = if_else(days_dos >= median_os$median[median_os$cond.2 == cond.2], "high", "low"),
pfs_dich = if_else(days_dos_pfs >= median_pfs$median[median_pfs$cond.2 == cond.2], "high", "low")) %>%
select(biosample_id, cond.2, os_dich, pfs_dich)
# Same dichotomization, but using medians computed only over non-cold samples.
median_os_responder_noncold <- df.clin %>% filter(!is.na(cond.2)) %>%
left_join(df.samples.mean %>% select(biosample_id, lcross.cat.biosample)) %>%
filter(lcross.cat.biosample != "cold") %>% group_by(cond.2) %>% summarise(median = median(days_dos))
median_pfs_responder_noncold <- df.clin %>% filter(!is.na(cond.2)) %>%
left_join(df.samples.mean %>% select(biosample_id, lcross.cat.biosample)) %>%
filter(lcross.cat.biosample != "cold") %>% group_by(cond.2) %>% summarise(median = median(days_dos_pfs))
survival.responder_noncold_dich <- df.clin %>% filter(!is.na(cond.2)) %>%
left_join(df.samples.mean %>% select(biosample_id, lcross.cat.biosample)) %>%
filter(lcross.cat.biosample != "cold") %>%
distinct(biosample_id, days_dos, days_dos_pfs, cond.2) %>% rowwise() %>%
mutate(os_dich = if_else(days_dos >= median_os_responder_noncold$median[median_os_responder_noncold$cond.2 == cond.2], "high", "low"),
pfs_dich = if_else(days_dos_pfs >= median_pfs_responder_noncold$median[median_pfs_responder_noncold$cond.2 == cond.2], "high", "low")) %>%
select(biosample_id, cond.2, os_dich, pfs_dich)
# Map the dichotomized survival categories back onto the SCE object.
sce$cond.2.dich.os <- survival.cond.2.dich$os_dich[match(sce$biosample_id, survival.cond.2.dich$biosample_id)]
sce$cond.2.dich.os <- factor(sce$cond.2.dich.os, levels = c("low", "high"))
sce$cond.2.dich.pfs <- survival.cond.2.dich$pfs_dich[match(sce$biosample_id, survival.cond.2.dich$biosample_id)]
sce$cond.2.dich.pfs <- factor(sce$cond.2.dich.pfs, levels = c("low", "high"))
sce$cond.2.dich.os_responder_noncold <- survival.responder_noncold_dich$os_dich[match(sce$biosample_id, survival.responder_noncold_dich$biosample_id)]
sce$cond.2.dich.os_responder_noncold <- factor(sce$cond.2.dich.os_responder_noncold, levels = c("low", "high"))
sce$cond.2.dich.pfs_responder_noncold <- survival.responder_noncold_dich$pfs_dich[match(sce$biosample_id, survival.responder_noncold_dich$biosample_id)]
sce$cond.2.dich.pfs_responder_noncold <- factor(sce$cond.2.dich.pfs_responder_noncold, levels = c("low", "high"))
# spicyR: pairwise spatial co-localization tests.
# cond.2 = ICI responder vs. non-responder.
spicyTest_cond.2 <- spicy(
sce[, (!is.na(sce$cond.2) & sce$cl_all_spec != "unassigned")],
condition = "cond.2",
subject = "biosample_id",
imageID = "image_id",
spatialCoords = c("Pos_X", "Pos_Y"),
cellType = "cl_all_spec",
BPPARAM = MulticoreParam(progressbar = T))
# cond.3 = ICI pretreated vs. naive.
spicyTest_cond.3 <- spicy(
sce[, (!is.na(sce$cond.3) & sce$cl_all_spec != "unassigned")],
condition = "cond.3",
subject = "biosample_id",
imageID = "image_id",
spatialCoords = c("Pos_X", "Pos_Y"),
cellType = "cl_all_spec",
BPPARAM = MulticoreParam(progressbar = T))
```
### For plot
```{r}
library(ggplotify)
# Default spicyR significance/heatmap views.
signifPlot(spicyTest_cond.2, fdr = F, cutoff = 0.05, breaks = c(-2, 2, 1)) +
scale_shape_manual(labels = c("responder", "non-responder"),
values = c(GroupA = "\u25D6", GroupB = "\u25D7")) +
ggtitle("All patients: responder vs. non-responder")
p1 <- as.ggplot(signifPlot(spicyTest_cond.2, type = "heatmap")) +
ggtitle("ICI-naive: responders vs. non-responders")
as.ggplot(signifPlot(spicyTest_cond.2, type = "heatmap")) +
ggtitle("ICI-naive: responders vs. non-responders")
# Rebuild the spicyR heatmap in ggplot for full customization.
celltype_order <- c("Tumor", "DPTc", "CD4Tc", "CD8Tc", "exh.CD8Tc", "GZB+ act.CD8Tc",
"Treg", "NK cells", "BnT", "B cells", "Plasma", "Myeloid",
"CD38+ Mf", "IDO+ Mf", "Arginase+ Neu", "Neutrophils",
"Vascular", "Astrocytes")
# Signed -log10(p): positive = attraction, negative = avoidance (cond.2).
spicyObject <- spicyTest_cond.2
ref.condition <- "responder"
df <- data.frame(coef = spicyObject$coefficient[[paste0("condition", ref.condition)]],
pval = spicyObject$p.value[[paste0("condition", ref.condition)]],
FDR = p.adjust(spicyObject$p.value[[paste0("condition", ref.condition)]], method = "fdr"),
from = spicyObject$comparisons$from,
to = spicyObject$comparisons$to) %>%
as_tibble() %>% mutate(log10.pval = if_else(coef >= 0, -log10(pval), +log10(pval)))
p.spicy.1 <- df %>% ggplot(aes((from), (to))) +
geom_tile(aes(fill = log10.pval), color = "black", lwd = 0.25, linetype = 1) +
scale_fill_gradient2(low = "#4575B4", mid = "white", high = "#D73027",
breaks = c(min(df$log10.pval), -2, -1, 0, 1, 2, max(df$log10.pval)),
labels = c("Avoidance\nin responders", -2, -1, 0, 1, 2, "Attractance\nin responders")) +
scale_x_discrete(limits = celltype_order) +
scale_y_discrete(limits = rev(celltype_order)) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
ggtitle("ICI: Responder vs. non-responder")
p.spicy.1
# Same, for cond.3 (ICI pretreated vs. naive).
spicyObject <- spicyTest_cond.3
ref.condition <- "ici_pretreated"
df <- data.frame(coef = spicyObject$coefficient[[paste0("condition", ref.condition)]],
pval = spicyObject$p.value[[paste0("condition", ref.condition)]],
FDR = p.adjust(spicyObject$p.value[[paste0("condition", ref.condition)]], method = "fdr"),
from = spicyObject$comparisons$from,
to = spicyObject$comparisons$to) %>%
as_tibble() %>% mutate(log10.pval = if_else(coef >= 0, -log10(pval), +log10(pval)))
p.spicy.4 <- df %>% ggplot(aes((from), (to))) +
geom_tile(aes(fill = log10.pval), color = "black", lwd = 0.25, linetype = 1) +
scale_fill_gradient2(low = "#4575B4", mid = "white", high = "#D73027",
breaks = c(min(df$log10.pval), -2, -1, 0, 1, max(df$log10.pval)),
labels = c("Avoidance\nin ici_pretreated", -2, -1, 0, 1, "Attractance\nin ici_pretreated")) +
scale_x_discrete(limits = celltype_order) +
scale_y_discrete(limits = rev(celltype_order)) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
ggtitle("ICI pretreated vs. naive")
p.spicy.4
```
## Survival of ICI-naive stratified by spatial infiltration pattern
```{r}
df.all <- colData(sce) %>% as_tibble() %>% distinct(biosample_id, .keep_all = T) %>%
select(id, biosample_id, contains("cond"))
clin <- metadata(sce)$clinical %>%
select(id, biosample_id, days_dos, days_dg, days_dos_pfs, status, status_pfs)
# ICI-naive samples only.
df.subset_all <- df.all %>%
filter(!is.na(cond.2)) %>%
left_join(df.samples.mean) %>% left_join(clin)
# Overall survival by spatial infiltration pattern.
fit <- survfit2(Surv(days_dos / 365.25, status) ~ lcross.cat.biosample, data = df.subset_all)
p.spat.os <- fit %>% ggsurvfit() + add_pvalue() + add_censor_mark() + add_quantile(y_value = 0.5) +
coord_cartesian(xlim = c(0, 8)) +
scale_y_continuous(limits = c(0, 1), labels = scales::percent, expand = c(0.01, 0)) +
scale_color_discrete(name = "CD8Tc infiltration pattern:",
labels = c("cold" = "cold", "low" = "localized", "high" = "dispersed")) +
theme_minimal() + xlab("Follow-up time, years") +
ggtitle("ICI-naive: Overall survival") + theme(text = element_text(size = 15))
(p.spat.os <- ggsurvfit_build(p.spat.os))
median_table <- fit %>% tidy_survfit() %>% group_by(strata) %>%
summarize(median_survival = time[which(estimate <= 0.5)[1]], .groups = "drop")
print(median_table)
p.spat.os
# Local progression-free survival by spatial infiltration pattern.
fit <- survfit2(Surv(days_dos_pfs / 365.25, status_pfs) ~ lcross.cat.biosample, data = df.subset_all)
p.spat.pfs <- fit %>% ggsurvfit() + add_pvalue() + add_censor_mark() + add_quantile(y_value = 0.5) +
coord_cartesian(xlim = c(0, 8)) +
scale_y_continuous(limits = c(0, 1), labels = scales::percent, expand = c(0.01, 0)) +
scale_color_discrete(name = "CD8Tc infiltration pattern:",
labels = c("cold" = "cold", "low" = "localized", "high" = "dispersed")) +
theme_minimal() + xlab("Follow-up time, years") +
ggtitle("ICI-naive: Local progression-free survival") + theme(text = element_text(size = 15))
(p.spat.pfs <- ggsurvfit_build(p.spat.pfs))
median_table <- fit %>% tidy_survfit() %>% group_by(strata) %>%
summarize(median_survival = time[which(estimate <= 0.5)[1]], .groups = "drop")
print(median_table)
p.spat.pfs
```
## Assemble
```{r, fig.width=20, fig.height=12}
# Main Figure 4 layout.
((((p.os + p.pfs) + plot_layout(ncol = 1, nrow = 2)) |
((p.da.1 + p.da.2 + plot_layout(widths = c(1.5, 0.5))))) +
plot_layout(widths = c(0.5, 3))) /
((p.spicy.1 + theme(text = element_text(size = 14)) |
(p.spat.os / p.spat.pfs) | plot_spacer()) + plot_layout(widths = c(2, 0.5, 1))) +
plot_annotation(tag_levels = list(c("A", "B", "C", "", "D", "E", "F")))
ggsave("figures/fig4.pdf", width = 20, height = 12, dpi = 200)
# Supplementary spicyR figure.
(p.spicy.4) + plot_annotation(tag_levels = list(c("A", "B")))
ggsave("figures/suppfig_spicyR.pdf", width = 10, height = 6, dpi = 200)
# Additional assembled views / single panels.
((p.da.1 + p.da.2 + plot_layout(widths = c(1.5, 1)))) +
(((p.spat.os / p.spat.pfs) | (plot_spacer())) + plot_layout(widths = c(1, 1.5))) +
plot_layout(heights = c(1, 1.5))
ggsave("figures/fig4_add.pdf", width = 15, height = 11, dpi = 300)
((p.os / p.pfs) | (p.spat.os / p.spat.pfs))
ggsave("figures/fig4_survival.pdf", width = 13, height = 8, dpi = 300)
wrap_plots(p.da.1.comb, p.da.2.comb) + plot_layout(widths = c(2, 1))
ggsave("figures/fig4_DA_revised.pdf", width = 12, height = 4.5, dpi = 300)
p.spicy.1 + theme(text = element_text(size = 14))
ggsave("figures/fig4_spicy.pdf", width = 9, height = 5, dpi = 300)
```
# Save workspace
```{r}
save.image(file = file.path(output_path, "FINALFIGURES.RData"))
```
# ------------------------------
# ------------------------------
# REVISION: Survival model
## Load workspace
```{r}
datadir <- "/mnt/projects_output_data/MBM"
output_path <- file.path(datadir, "data")
load(file = file.path(output_path, "FINALFIGURES.RData"))
```
# CLINICAL OUTCOME
## Collect all parameters
```{r}
library(survival)
library(survcomp)
library(survminer)
library(glmnet)
library(pROC)
library(caret)
# --- Spatial features: cellular-neighbourhood proportions ---
# Image-level CN frequencies.
df.cn.sum_image <- colData(sce.subset) %>% as_tibble() %>%
group_by(image_id, cn_celltypes, .drop = FALSE) %>%
summarise(count = n()) %>% mutate(prop = prop.table(count)) %>%
group_by(cn_celltypes) %>% mutate(median = median(prop))
# Combine image-level CN frequencies with the per-image infiltration category.
df.comb_image <- df.cn.sum_image %>%
left_join(df.images %>% select(image_id, lcross.cat.image)) %>%
filter(!is.na(lcross.cat.image))
# Sample-level CN proportions (one row per sample).
df.spatialFeatures_CN <- df.comb %>% ungroup() %>%
select(biosample_id, cn_celltypes, prop, median) %>%
pivot_wider(id_cols = "biosample_id", names_from = "cn_celltypes", values_from = "prop") %>%
replace(is.na(.), 0) %>%
column_to_rownames("biosample_id")
# Sample-level L-cross infiltration category.
df.spatialFeatures_lcross <- df.comb %>% ungroup() %>%
distinct(biosample_id, lcross.cat.biosample) %>%
column_to_rownames("biosample_id")
# --- Cell-type proportions per sample ---
df.celltypeProportions <- colData(sce) %>% as_tibble() %>%
group_by(biosample_id, cl_all) %>%
dplyr::summarise(count = n()) %>% mutate(pct = prop.table(count)) %>%
ungroup() %>%
pivot_wider(id_cols = "biosample_id", names_from = "cl_all", values_from = "pct") %>%
replace(is.na(.), 0) %>%
column_to_rownames("biosample_id")
# --- Clinical features ---
cohort_ids <- colData(sce) %>% as_tibble() %>% distinct(biosample_id, image_id)
cohort_clin.data <- cohort_ids %>% distinct(biosample_id) %>%
left_join(metadata(sce)$clinical) %>%
filter(biosample_id %in% cohort_ids$biosample_id)
df.clinicalData <- cohort_clin.data %>%
column_to_rownames("biosample_id") %>%
select(age, sex, braf, nras, rt_bm_postop, ici_postop, extracranial_control_dos)
# --- Survival outcome ---
survival_outcome <- cohort_clin.data %>%
select(biosample_id, days_dos, status, days_dos_pfs, status_pfs) %>%
column_to_rownames("biosample_id")
# --- Condition labels and cohort alignment ---
cohort.conditions <- colData(sce) %>% as_tibble() %>% distinct(biosample_id, .keep_all = T) %>%
select(biosample_id, contains("cond.")) %>% as.data.frame()
rownames(cohort.conditions) <- cohort.conditions$biosample_id
survival_outcome <- survival_outcome[cohort.conditions$biosample_id, ]
```
# Multivariate penalized Cox proportional hazards model
A ridge-penalized Cox model (glmnet, alpha = 0) was fitted jointly on four feature
blocks (clinical variables, cell-type proportions, cellular-neighbourhood spatial
features, and L-cross spatial features). Factor variables were full-rank dummy
encoded; missing values were imputed (median for continuous, mode for categorical);
zero-variance features were dropped; and all features were z-scored, so hazard
ratios reflect the change in hazard per one standard deviation. Ridge regression
keeps all features with continuously shrunk coefficients (no hard selection). The
penalty lambda was selected by k-fold cross-validated partial-likelihood deviance
(k = 10, reduced to 5 if events < 30). Coefficient uncertainty was quantified by a
non-parametric bootstrap (refitting at the fixed lambda); percentile 95% CIs and
empirical two-sided p-values (FDR-adjusted) were derived from the bootstrap
distribution. Discrimination was summarised by Harrell's C-index, and results are
shown as forest plots on a log scale.
## PFS
```{r}
# Survival endpoint: local progression-free survival.
outcome.variable_days <- "days_dos_pfs"
outcome.variable_status <- "status_pfs"
survival_outcome.subset <- survival_outcome %>%
select(!!rlang::sym(outcome.variable_days), !!rlang::sym(outcome.variable_status)) %>% drop_na()
cohort.subset <- cohort.conditions
df.clinicalData.subset <- df.clinicalData
# Align all feature blocks to a common set of samples.
common_samples <- Reduce(intersect, list(
rownames(df.clinicalData.subset),
rownames(df.celltypeProportions),
rownames(df.spatialFeatures_CN),
rownames(df.spatialFeatures_lcross),
rownames(cohort.subset),
rownames(survival_outcome.subset)
))
clin <- df.clinicalData.subset[common_samples, ]
ctp <- df.celltypeProportions[common_samples, ]
cn <- df.spatialFeatures_CN[common_samples, ]
lcross <- df.spatialFeatures_lcross[common_samples, , drop = F]
days <- survival_outcome.subset[common_samples, outcome.variable_days]
status <- survival_outcome.subset[common_samples, outcome.variable_status]
# Impute clinical NAs (median for numeric, mode for categorical).
clin_imputed <- clin
for (col in colnames(clin_imputed)) {
if (is.numeric(clin_imputed[[col]])) {
clin_imputed[[col]][is.na(clin_imputed[[col]])] <- median(clin_imputed[[col]], na.rm = TRUE)
} else {
mode_val <- names(sort(table(clin_imputed[[col]]), decreasing = TRUE))[1]
clin_imputed[[col]][is.na(clin_imputed[[col]])] <- mode_val
}
}
# Dummy-encode factor blocks and drop the intercept.
clin <- model.matrix(~ . , data = clin_imputed)[, -1]
lcross <- model.matrix(~ . , data = lcross)[, -1]
set.seed(42)
# Assemble and z-score the feature matrix (drop zero-variance columns).
X <- cbind(clin, ctp, cn, lcross)
X <- X[, apply(X, 2, var, na.rm = TRUE) > 0]
X_scaled <- scale(X)
y <- Surv(days, status)
# Annotate each feature with its originating block.
block_map <- bind_rows(
data.frame(feature = colnames(clin), block = "Clinical", stringsAsFactors = FALSE),
data.frame(feature = colnames(ctp), block = "Cell type proportions", stringsAsFactors = FALSE),
data.frame(feature = colnames(cn), block = "Spatial", stringsAsFactors = FALSE),
data.frame(feature = colnames(lcross), block = "Spatial", stringsAsFactors = FALSE)
) %>% filter(feature %in% colnames(X_scaled))
# Ridge Cox with cross-validated lambda (deviance); fall back to AIC if CV fails.
ALPHA <- 0
n_events <- sum(status)
nfolds <- if (n_events < 30) 5 else 10
cat(sprintf("Samples: %d | Events: %d | Features: %d | CV folds: %d\n",
nrow(X_scaled), n_events, ncol(X_scaled), nfolds))
cv_fit <- tryCatch(
cv.glmnet(X_scaled, y, family = "cox", alpha = ALPHA, nfolds = nfolds, type.measure = "deviance"),
error = function(e) {
message("cv.glmnet failed: ", conditionMessage(e), "\nFalling back to glmnet + AIC.")
NULL
}
)
if (is.null(cv_fit)) {
fit_nocv <- glmnet(X_scaled, y, family = "cox", alpha = ALPHA)
aic_vals <- deviance(fit_nocv) + 2 * fit_nocv$df
lambda_chosen <- fit_nocv$lambda[which.min(aic_vals)]
cv_fit <- list(lambda.min = lambda_chosen, lambda.1se = lambda_chosen,
lambda = fit_nocv$lambda, cvm = aic_vals, glmnet.fit = fit_nocv)
class(cv_fit) <- "cv.glmnet.fallback"
cat(sprintf("AIC-chosen lambda: %.5f\n", lambda_chosen))
} else {
lambda_chosen <- cv_fit$lambda.min
cat(sprintf("CV-chosen lambda (lambda.min): %.5f\n", lambda_chosen))
}
# Bootstrap coefficient distribution at the fixed lambda.
N_BOOT <- 1000
CI_LEVEL <- 0.95
alpha_ci <- (1 - CI_LEVEL) / 2
n <- nrow(X_scaled)
boot_coefs <- matrix(NA_real_, nrow = N_BOOT, ncol = ncol(X_scaled),
dimnames = list(NULL, colnames(X_scaled)))
cat(sprintf("\nRunning %d bootstrap iterations (Ridge, alpha=0)...\n", N_BOOT))
pb <- txtProgressBar(min = 0, max = N_BOOT, style = 3)
for (i in seq_len(N_BOOT)) {
idx <- sample(n, n, replace = TRUE)
y_boot <- y[idx]
if (sum(y_boot[, "status"]) < 2) next
if (length(unique(y_boot[y_boot[, "status"] == 1, "time"])) < 2) next
fit_boot <- tryCatch(
glmnet(X_scaled[idx, , drop = FALSE], y_boot, family = "cox", alpha = ALPHA, lambda = lambda_chosen),
error = function(e) NULL
)
if (!is.null(fit_boot)) boot_coefs[i, ] <- as.vector(coef(fit_boot, s = lambda_chosen))
setTxtProgressBar(pb, i)
}
close(pb)
n_success <- sum(complete.cases(boot_coefs))
cat(sprintf("\nSuccessful bootstrap iterations: %d / %d\n", n_success, N_BOOT))
# Point estimates, bootstrap CIs, and empirical p-values.
coef_full <- as.vector(coef(cv_fit, s = lambda_chosen))
names(coef_full) <- colnames(X_scaled)
ci_lower <- apply(boot_coefs, 2, quantile, probs = alpha_ci, na.rm = TRUE)
ci_upper <- apply(boot_coefs, 2, quantile, probs = 1 - alpha_ci, na.rm = TRUE)
emp_p <- sapply(colnames(X_scaled), function(f) {
b <- boot_coefs[, f]; b <- b[!is.na(b)]
if (length(b) == 0) return(NA_real_)
p_raw <- if (coef_full[f] >= 0) mean(b <= 0) else mean(b >= 0)
pmin(2 * p_raw, 1)
})
emp_p <- pmax(emp_p, 1 / n_success) # resolution floor
# C-index on the full data.
lp <- as.vector(predict(cv_fit, newx = X_scaled, s = lambda_chosen, type = "link"))
c_index <- as.numeric(concordance(y ~ lp)$concordance)
cat(sprintf("C-index (full data): %.3f\n", c_index))
# Results table with FDR and CI-crossing flag.
results <- data.frame(
feature = colnames(X_scaled),
HR = exp(coef_full),
HR_lower = exp(ci_lower),
HR_upper = exp(ci_upper),
p = emp_p,
stringsAsFactors = FALSE
) %>%
left_join(block_map, by = "feature") %>%
mutate(block = coalesce(block, "Unknown")) %>%
mutate(p_adj = p.adjust(p, method = "BH"),
ci_cross = HR_lower < 1 & HR_upper > 1,
sig = case_when(p_adj < 0.1 ~ "*", TRUE ~ ""),
label = paste0(feature, sig)) %>%
arrange(block, HR) %>%
mutate(row_id = factor(seq_len(n()), levels = rev(seq_len(n()))))
cat(sprintf("Features with CI crossing HR=1: %d / %d\n", sum(results$ci_cross), nrow(results)))
# Forest plot (grey = CI crosses HR=1).
block_colors <- c("Clinical" = "#1B9E77", "Cell type proportions" = "#D95F02", "Spatial" = "#7570B3")
results <- results %>% mutate(pt_colour = if_else(ci_cross, "grey70", block_colors[block]))
x_lo <- min(results$HR_lower, na.rm = TRUE) * 0.90
x_hi <- max(results$HR_upper, na.rm = TRUE) * 1.10
active_breaks <- exp(pretty(log(c(x_lo, x_hi)), n = 5))
active_breaks <- signif(active_breaks, 2)
active_breaks <- sort(unique(c(1, active_breaks)))
active_breaks <- active_breaks[active_breaks >= x_lo & active_breaks <= x_hi]
size_breaks <- signif(quantile(abs(log(results$HR)), probs = c(0.25, 0.5, 0.75, 1.0), na.rm = TRUE), 2)
size_breaks <- sort(unique(size_breaks[size_breaks > 0]))
subtitle_txt <- sprintf(
paste0("Ridge Cox (\u03b1=0, \u03bb=%.4f) | ALL %d features shown | ",
"C-index = %.3f | %d bootstrap resamples | ",
"Grey = CI crosses HR=1 | * FDR<0.05 ** <0.01 *** <0.001"),
lambda_chosen, nrow(results), c_index, N_BOOT)
p_forest <- ggplot(results, aes(y = row_id)) +
geom_vline(xintercept = 1, linetype = "dashed", colour = "grey40", linewidth = 0.5) +
geom_errorbarh(aes(xmin = HR_lower, xmax = HR_upper, colour = pt_colour), height = 0.25, linewidth = 0.55) +
geom_point(aes(x = HR, colour = pt_colour, size = abs(log(HR))), shape = 18) +
geom_text(aes(x = HR_upper, label = sig), hjust = -0.2, vjust = 0.5, size = 3, colour = "black") +
facet_grid(block ~ ., scales = "free_y", space = "free_y", switch = "y") +
scale_colour_identity() +
scale_size_continuous(name = "|log(HR)|") +
scale_x_log10(limits = c(x_lo, x_hi), breaks = active_breaks, labels = as.character(active_breaks)) +
scale_y_discrete(labels = setNames(results$label, results$row_id)) +
labs(title = "Multivariate Ridge Cox PH \u2013 Forest Plot (all features)",
subtitle = subtitle_txt,
x = "Hazard Ratio (log scale, per-SD)", y = NULL,
caption = paste0("All features shown; no hard selection threshold applied.\n",
sprintf("%d%% bootstrap percentile CIs (n=%d resamples). ", round(CI_LEVEL * 100), N_BOOT),
"Grey = CI includes HR=1. Stars = FDR-adjusted empirical p-value.")) +
theme_bw(base_size = 11) +
theme(strip.placement = "outside",
strip.background = element_rect(fill = "grey93", colour = NA),
strip.text.y.left = element_text(angle = 0, face = "bold", size = 9),
panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(),
panel.spacing = unit(0.35, "lines"), axis.text.y = element_text(size = 7.5),
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(size = 7, colour = "grey35"),
plot.caption = element_text(size = 7, colour = "grey50"),
legend.position = "bottom")
p_forest
results_pfs <- results
```
## OS
Identical workflow to the PFS block above, using overall survival as the endpoint.
```{r}
# Survival endpoint: overall survival.
outcome.variable_days <- "days_dos"
outcome.variable_status <- "status"
survival_outcome.subset <- survival_outcome %>%
select(!!rlang::sym(outcome.variable_days), !!rlang::sym(outcome.variable_status)) %>% drop_na()
cohort.subset <- cohort.conditions
df.clinicalData.subset <- df.clinicalData
# Align all feature blocks to a common set of samples.
common_samples <- Reduce(intersect, list(
rownames(df.clinicalData.subset),
rownames(df.celltypeProportions),
rownames(df.spatialFeatures_CN),
rownames(df.spatialFeatures_lcross),
rownames(cohort.subset),
rownames(survival_outcome.subset)
))
clin <- df.clinicalData.subset[common_samples, ]
ctp <- df.celltypeProportions[common_samples, ]
cn <- df.spatialFeatures_CN[common_samples, ]
lcross <- df.spatialFeatures_lcross[common_samples, , drop = F]
days <- survival_outcome.subset[common_samples, outcome.variable_days]
status <- survival_outcome.subset[common_samples, outcome.variable_status]
# Impute clinical NAs (median for numeric, mode for categorical).
clin_imputed <- clin
for (col in colnames(clin_imputed)) {
if (is.numeric(clin_imputed[[col]])) {
clin_imputed[[col]][is.na(clin_imputed[[col]])] <- median(clin_imputed[[col]], na.rm = TRUE)
} else {
mode_val <- names(sort(table(clin_imputed[[col]]), decreasing = TRUE))[1]
clin_imputed[[col]][is.na(clin_imputed[[col]])] <- mode_val
}
}
clin <- model.matrix(~ . , data = clin_imputed)[, -1]
lcross <- model.matrix(~ . , data = lcross)[, -1]
set.seed(42)
X <- cbind(clin, ctp, cn, lcross)
X <- X[, apply(X, 2, var, na.rm = TRUE) > 0]
X_scaled <- scale(X)
y <- Surv(days, status)
block_map <- bind_rows(
data.frame(feature = colnames(clin), block = "Clinical", stringsAsFactors = FALSE),
data.frame(feature = colnames(ctp), block = "Cell type proportions", stringsAsFactors = FALSE),
data.frame(feature = colnames(cn), block = "Spatial", stringsAsFactors = FALSE),
data.frame(feature = colnames(lcross), block = "Spatial", stringsAsFactors = FALSE)
) %>% filter(feature %in% colnames(X_scaled))
ALPHA <- 0
n_events <- sum(status)
nfolds <- if (n_events < 30) 5 else 10
cat(sprintf("Samples: %d | Events: %d | Features: %d | CV folds: %d\n",
nrow(X_scaled), n_events, ncol(X_scaled), nfolds))
cv_fit <- tryCatch(
cv.glmnet(X_scaled, y, family = "cox", alpha = ALPHA, nfolds = nfolds, type.measure = "deviance"),
error = function(e) {
message("cv.glmnet failed: ", conditionMessage(e), "\nFalling back to glmnet + AIC.")
NULL
}
)
if (is.null(cv_fit)) {
fit_nocv <- glmnet(X_scaled, y, family = "cox", alpha = ALPHA)
aic_vals <- deviance(fit_nocv) + 2 * fit_nocv$df
lambda_chosen <- fit_nocv$lambda[which.min(aic_vals)]
cv_fit <- list(lambda.min = lambda_chosen, lambda.1se = lambda_chosen,
lambda = fit_nocv$lambda, cvm = aic_vals, glmnet.fit = fit_nocv)
class(cv_fit) <- "cv.glmnet.fallback"
cat(sprintf("AIC-chosen lambda: %.5f\n", lambda_chosen))
} else {
lambda_chosen <- cv_fit$lambda.min
cat(sprintf("CV-chosen lambda (lambda.min): %.5f\n", lambda_chosen))
}
N_BOOT <- 1000
CI_LEVEL <- 0.95
alpha_ci <- (1 - CI_LEVEL) / 2
n <- nrow(X_scaled)
boot_coefs <- matrix(NA_real_, nrow = N_BOOT, ncol = ncol(X_scaled),
dimnames = list(NULL, colnames(X_scaled)))
cat(sprintf("\nRunning %d bootstrap iterations (Ridge, alpha=0)...\n", N_BOOT))
pb <- txtProgressBar(min = 0, max = N_BOOT, style = 3)
for (i in seq_len(N_BOOT)) {
idx <- sample(n, n, replace = TRUE)
y_boot <- y[idx]
if (sum(y_boot[, "status"]) < 2) next
if (length(unique(y_boot[y_boot[, "status"] == 1, "time"])) < 2) next
fit_boot <- tryCatch(
glmnet(X_scaled[idx, , drop = FALSE], y_boot, family = "cox", alpha = ALPHA, lambda = lambda_chosen),
error = function(e) NULL
)
if (!is.null(fit_boot)) boot_coefs[i, ] <- as.vector(coef(fit_boot, s = lambda_chosen))
setTxtProgressBar(pb, i)
}
close(pb)
n_success <- sum(complete.cases(boot_coefs))
cat(sprintf("\nSuccessful bootstrap iterations: %d / %d\n", n_success, N_BOOT))
coef_full <- as.vector(coef(cv_fit, s = lambda_chosen))
names(coef_full) <- colnames(X_scaled)
ci_lower <- apply(boot_coefs, 2, quantile, probs = alpha_ci, na.rm = TRUE)
ci_upper <- apply(boot_coefs, 2, quantile, probs = 1 - alpha_ci, na.rm = TRUE)
emp_p <- sapply(colnames(X_scaled), function(f) {
b <- boot_coefs[, f]; b <- b[!is.na(b)]
if (length(b) == 0) return(NA_real_)
p_raw <- if (coef_full[f] >= 0) mean(b <= 0) else mean(b >= 0)
pmin(2 * p_raw, 1)
})
emp_p <- pmax(emp_p, 1 / n_success)
lp <- as.vector(predict(cv_fit, newx = X_scaled, s = lambda_chosen, type = "link"))
c_index <- as.numeric(concordance(y ~ lp)$concordance)
cat(sprintf("C-index (full data): %.3f\n", c_index))
results <- data.frame(
feature = colnames(X_scaled),
HR = exp(coef_full),
HR_lower = exp(ci_lower),
HR_upper = exp(ci_upper),
p = emp_p,
stringsAsFactors = FALSE
) %>%
left_join(block_map, by = "feature") %>%
mutate(block = coalesce(block, "Unknown")) %>%
mutate(p_adj = p.adjust(p, method = "BH"),
ci_cross = HR_lower < 1 & HR_upper > 1,
sig = case_when(p_adj < 0.1 ~ "*", TRUE ~ ""),
label = paste0(feature, sig)) %>%
arrange(block, HR) %>%
mutate(row_id = factor(seq_len(n()), levels = rev(seq_len(n()))))
cat(sprintf("Features with CI crossing HR=1: %d / %d\n", sum(results$ci_cross), nrow(results)))
block_colors <- c("Clinical" = "#1B9E77", "Cell type proportions" = "#D95F02", "Spatial" = "#7570B3")
results <- results %>% mutate(pt_colour = if_else(ci_cross, "grey70", block_colors[block]))
x_lo <- min(results$HR_lower, na.rm = TRUE) * 0.90
x_hi <- max(results$HR_upper, na.rm = TRUE) * 1.10
active_breaks <- exp(pretty(log(c(x_lo, x_hi)), n = 5))
active_breaks <- signif(active_breaks, 2)
active_breaks <- sort(unique(c(1, active_breaks)))
active_breaks <- active_breaks[active_breaks >= x_lo & active_breaks <= x_hi]
size_breaks <- signif(quantile(abs(log(results$HR)), probs = c(0.25, 0.5, 0.75, 1.0), na.rm = TRUE), 2)
size_breaks <- sort(unique(size_breaks[size_breaks > 0]))
subtitle_txt <- sprintf(
paste0("Ridge Cox (\u03b1=0, \u03bb=%.4f) | ALL %d features shown | ",
"C-index = %.3f | %d bootstrap resamples | ",
"Grey = CI crosses HR=1 | * FDR<0.05 ** <0.01 *** <0.001"),
lambda_chosen, nrow(results), c_index, N_BOOT)
p_forest <- ggplot(results, aes(y = row_id)) +
geom_vline(xintercept = 1, linetype = "dashed", colour = "grey40", linewidth = 0.5) +
geom_errorbarh(aes(xmin = HR_lower, xmax = HR_upper, colour = pt_colour), height = 0.25, linewidth = 0.55) +
geom_point(aes(x = HR, colour = pt_colour, size = abs(log(HR))), shape = 18) +
geom_text(aes(x = HR_upper, label = sig), hjust = -0.2, vjust = 0.5, size = 3, colour = "black") +
facet_grid(block ~ ., scales = "free_y", space = "free_y", switch = "y") +
scale_colour_identity() +
scale_size_continuous(name = "|log(HR)|") +
scale_x_log10(limits = c(x_lo, x_hi), breaks = active_breaks, labels = as.character(active_breaks)) +
scale_y_discrete(labels = setNames(results$label, results$row_id)) +
labs(title = "Multivariate Ridge Cox PH \u2013 Forest Plot (all features)",
subtitle = subtitle_txt,
x = "Hazard Ratio (log scale, per-SD)", y = NULL,
caption = paste0("All features shown; no hard selection threshold applied.\n",
sprintf("%d%% bootstrap percentile CIs (n=%d resamples). ", round(CI_LEVEL * 100), N_BOOT),
"Grey = CI includes HR=1. Stars = FDR-adjusted empirical p-value.")) +
theme_bw(base_size = 11) +
theme(strip.placement = "outside",
strip.background = element_rect(fill = "grey93", colour = NA),
strip.text.y.left = element_text(angle = 0, face = "bold", size = 9),
panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(),
panel.spacing = unit(0.35, "lines"), axis.text.y = element_text(size = 7.5),
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(size = 7, colour = "grey35"),
plot.caption = element_text(size = 7, colour = "grey50"),
legend.position = "bottom")
p_forest
results_os <- results
```
### Combined plot
```{r, fig.width = 10, figh.height = 25}
library(ggh4x)
# Combine OS and PFS results and order panels so colours match the block scheme.
results_os$outcome <- "Overall survival"
results_pfs$outcome <- "Local progression-free surrival"
results_all <- rbind(results_os, results_pfs)
results_all$block <- factor(results_all$block, levels = names(block_colors))
# Relabel features with readable names; append "*" for FDR < 0.1.
results_all$label <- results_all$feature
results_all <- results_all %>%
mutate(label = case_when(
label == "sexM" ~ "Sex male",
label == "ici_postopyes" ~ "ICI postop",
label == "nrasmut" ~ "NRAS mutated",
label == "brafmut" ~ "BRAF mutated",
label == "rt_bm_postopyes" ~ "Radiotherapy postop",
label == "age" ~ "Age",
label == "extracranial_control_dosyes" ~ "Extracranial controlled disease",
label == "lcross.cat.biosamplelow" ~ "CD8Tc infiltration: Localized",
label == "lcross.cat.biosamplehigh" ~ "CD8Tc infiltration: Dispersed",
label == "cn_1" ~ "CN: Perivascular",
label == "cn_2" ~ "CN: Tumor core",
label == "cn_3" ~ "CN: Tumor CD8Tc border",
label == "cn_4" ~ "CN: Myeloid enriched",
label == "cn_5" ~ "CN: CD8Tc core",
label == "cn_6" ~ "CN: B cell, Neutrophil, Mf enriched",
label == "cn_7" ~ "CN: CD4Tc/Treg, Plasma cell enriched",
TRUE ~ label)) %>%
mutate(label = if_else(sig == "*", paste0(label, "*"), label))
p.multivariate.coxph <- ggplot(results_all, aes(y = feature)) +
geom_vline(xintercept = 1, linetype = "dashed", colour = "grey40", linewidth = 0.5) +
geom_errorbarh(aes(xmin = HR_lower, xmax = HR_upper, colour = block, alpha = p < 0.05),
height = 0.25, linewidth = 0.55) +
geom_point(aes(x = HR, colour = block, alpha = p < 0.05, size = -log10(p), fill = "")) +
scale_fill_manual(name = "* FDR < 0.1", values = "transparent",
guide = guide_legend(override.aes = list(fill = NA, color = NA),
theme = theme(legend.key = element_blank()))) +
geom_text(data = results_all %>% filter(sig == "*"),
aes(x = HR, label = "*"), vjust = .8, size = 4, colour = "black") +
facet_grid2(block ~ outcome, scales = "free_y", space = "free_y", switch = "y",
strip = strip_themed(background_y = elem_list_rect(fill = block_colors),
text_y = elem_list_text(colour = c("white")))) +
scale_colour_manual(values = block_colors, guide = "none") +
scale_alpha_discrete(name = "p-value < 0.05", range = c(0.4, 1)) +
scale_size_continuous(name = "-log10(p-value)") +
scale_y_discrete(labels = setNames(str_wrap(results_all$label, width = 50), results_all$feature)) +
labs(title = "Multivariate Ridge Cox PH \u2013 Forest Plot (all features)",
x = "Hazard Ratio (per-SD)", y = NULL,
caption = paste0("All features shown; no hard selection threshold applied.\n",
sprintf("%d%% bootstrap percentile CIs (n=%d resamples). ", round(CI_LEVEL * 100), N_BOOT))) +
theme(plot.title = element_text(face = "bold"),
plot.subtitle = element_text(size = 7, colour = "grey35"),
plot.caption = element_text(size = 7, colour = "grey50"),
legend.position = "right",
axis.text.y = element_text(hjust = 0, lineheight = .7))
p.multivariate.coxph
```
# Regularised Cox model comparison (Ridge)
Ridge-regularised Cox models (glmnet, alpha = 0) were fitted for all additive
combinations of the four feature blocks (clinical, cell-type proportions, CN, and
L-cross spatial), yielding the candidate models below. Given the small cohort
(n < 50), leave-one-out cross-validation was used to select lambda and Ridge was
preferred over Lasso for coefficient stability. Discrimination was summarised by
the concordance index (C-index) with bootstrap 95% CIs (1,000 resamples); each
model was compared to the clinical baseline via cindex.comp(), a likelihood-ratio
test on the linear predictors, and a permutation test for the best model.
## PFS
```{r}
outcome.variable_days <- "days_dos_pfs"
outcome.variable_status <- "status_pfs"
survival_outcome.subset <- survival_outcome %>%
select(!!rlang::sym(outcome.variable_days), !!rlang::sym(outcome.variable_status)) %>% drop_na()
# Align all feature blocks to common samples.
common_samples <- Reduce(intersect, list(
rownames(df.clinicalData),
rownames(df.celltypeProportions),
rownames(df.spatialFeatures_CN),
rownames(df.spatialFeatures_lcross),
rownames(survival_outcome.subset)
))
clin <- df.clinicalData[common_samples, ]
ctp <- df.celltypeProportions[common_samples, ]
cn <- df.spatialFeatures_CN[common_samples, ]
lcross <- df.spatialFeatures_lcross[common_samples, , drop = F]
n <- length(common_samples)
# Impute clinical NAs (median numeric, mode categorical) and dummy-encode.
clin_imputed <- clin
for (col in colnames(clin_imputed)) {
if (is.numeric(clin_imputed[[col]])) {
clin_imputed[[col]][is.na(clin_imputed[[col]])] <- median(clin_imputed[[col]], na.rm = TRUE)
} else {
mode_val <- names(sort(table(clin_imputed[[col]]), decreasing = TRUE))[1]
clin_imputed[[col]][is.na(clin_imputed[[col]])] <- mode_val
}
}
clin_mat <- model.matrix(~ . - 1, data = clin_imputed)
lcross <- model.matrix(~ . - 1, data = lcross)
days <- survival_outcome.subset[common_samples, outcome.variable_days]
status <- survival_outcome.subset[common_samples, outcome.variable_status]
surv <- Surv(days, status)
cat(sprintf("Proceeding with %d samples after NA removal.\n", n))
# Candidate feature matrices: each block alone and all additive combinations.
spatial_mat <- cbind(as.matrix(cn), lcross)
mat <- list(
"Clinical" = clin_mat,
"Cell type proportions" = as.matrix(ctp),
"Spatial" = spatial_mat,
"Clinical + Cell type proportions" = cbind(clin_mat, as.matrix(ctp)),
"Clinical + Spatial" = cbind(clin_mat, spatial_mat),
"Cell type proportions + Spatial" = cbind(as.matrix(ctp), spatial_mat),
"Clinical + Cell type proportions + Spatial" = cbind(clin_mat, as.matrix(ctp), spatial_mat)
)
set.seed(42)
# Fit Ridge Cox (LOO-CV), extract C-index (Noether) + bootstrap CI + logLik/AIC.
fit_model <- function(feature_mat, days, status, surv, n, B = 1000) {
cv_fit <- cv.glmnet(feature_mat, surv, family = "cox", alpha = 0, nfolds = n)
lp <- as.vector(predict(cv_fit, newx = feature_mat, s = "lambda.min"))
# Unpenalized Cox on the linear predictor -> logLik / AIC.
std_fit <- coxph(Surv(days, status) ~ lp)
ll <- logLik(std_fit)
aic_val <- AIC(std_fit)
# C-index (Noether method, required for cindex.comp()).
ci_noether <- concordance.index(x = lp, surv.time = days, surv.event = status, method = "noether")
# Bootstrap 95% CI for the C-index.
boot_c <- replicate(B, {
idx <- sample(n, replace = TRUE)
tryCatch(concordance.index(x = lp[idx], surv.time = days[idx], surv.event = status[idx])$c.index,
error = function(e) NA, warning = function(w) NA)
})
n_failed <- sum(is.na(boot_c))
if (n_failed > 0) message(sprintf("Bootstrap: %d/%d resamples failed and were skipped.", n_failed, B))
ci_noether$lower <- quantile(boot_c, 0.025, na.rm = TRUE)
ci_noether$upper <- quantile(boot_c, 0.975, na.rm = TRUE)
list(lp = lp, ci = ci_noether, logLik = ll, aic = aic_val)
}
results <- lapply(mat, fit_model, days = days, status = status, surv = surv, n = n, B = 1000)
# Summary table (C-index, CI, AIC, p) ordered by C-index.
summary_df <- data.frame(
Model = names(results),
C_index = sapply(results, \(r) round(r$ci$c.index, 3)),
CI_lower = sapply(results, \(r) round(r$ci$lower, 3)),
CI_upper = sapply(results, \(r) round(r$ci$upper, 3)),
AIC = sapply(results, \(r) round(r$aic, 2)),
p_value = sapply(results, \(r) signif(r$ci$p.value, 3))
)
summary_df <- summary_df[order(-summary_df$C_index), ]
print(summary_df, row.names = FALSE)
# Pairwise C-index comparison vs. the clinical baseline.
baseline <- results[["Clinical"]]$ci
cat("\n-- C-index comparison vs. Clinical baseline --\n")
for (nm in setdiff(names(results), "Clinical")) {
comp <- cindex.comp(baseline, results[[nm]]$ci)
cat(sprintf(" Clinical vs %-35s p = %.4f\n", nm, comp$p.value))
}
# Likelihood-ratio test vs. clinical (df = 1).
cat("\n-- Likelihood Ratio Test (LRT) vs. Clinical --\n")
ll_baseline <- results[["Clinical"]]$logLik
for (nm in setdiff(names(results), "Clinical")) {
ll_complex <- results[[nm]]$logLik
lrt_stat <- as.numeric(2 * (ll_complex - ll_baseline))
lrt_p <- pchisq(lrt_stat, df = 1, lower.tail = FALSE)
cat(sprintf(" Clinical vs %-20s | dLL: %6.2f | LRT p = %.4f\n",
nm, (ll_complex - ll_baseline), lrt_p))
}
# Permutation test: best combined model vs. clinical.
best_name <- summary_df$Model[1]
if (best_name != "Clinical") {
observed_diff <- results[[best_name]]$ci$c.index - results[["Clinical"]]$ci$c.index
perm_diffs <- replicate(1000, {
idx <- sample(n)
c1 <- concordance.index(results[["Clinical"]]$lp[idx], days, status, method = "noether")$c.index
c2 <- concordance.index(results[[best_name]]$lp[idx], days, status, method = "noether")$c.index
c2 - c1
})
perm_p <- mean(abs(perm_diffs) >= abs(observed_diff))
cat(sprintf("\nPermutation test - Clinical vs %s: dC = %.3f, p = %.4f\n",
best_name, observed_diff, perm_p))
}
# C-index forest plot.
summary_df$Model <- factor(summary_df$Model, levels = rev(summary_df$Model))
ggplot(summary_df, aes(x = C_index, y = Model)) +
geom_vline(xintercept = 0.5, linetype = "dashed", color = "grey60") +
geom_errorbarh(aes(xmin = CI_lower, xmax = CI_upper), height = 0.25, color = "grey40") +
geom_point(size = 3, color = "#E64B35") +
geom_vline(xintercept = results[[best_name]]$ci$c.index) +
geom_vline(xintercept = results[["Clinical"]]$ci$c.index) +
scale_x_continuous(limits = c(0.3, 1), breaks = seq(0.3, 1, 0.1)) +
labs(x = "C-index (bootstrap 95% CI)", y = NULL, title = "Survival prediction: model comparison") +
theme_classic(base_size = 13) +
theme(axis.text.y = element_text(hjust = 1))
# Kaplan-Meier risk groups (median split of linear predictor): clinical vs. best.
risk_group <- function(lp) factor(ifelse(lp > median(lp), "High", "Low"))
df_plot <- data.frame(days = days, status = status,
risk_clin = risk_group(results[["Clinical"]]$lp),
risk_best = risk_group(results[[best_name]]$lp))
p1 <- ggsurvplot(survfit(Surv(days, status) ~ risk_clin, data = df_plot),
data = df_plot, pval = TRUE, risk.table = TRUE,
title = "Clinical only", palette = c("#E64B35", "#4DBBD5"))
p2 <- ggsurvplot(survfit(Surv(days, status) ~ risk_best, data = df_plot),
data = df_plot, pval = TRUE, risk.table = TRUE,
title = best_name, palette = c("#E64B35", "#4DBBD5"))
arrange_ggsurvplots(list(p1, p2), ncol = 2)
# Ridge coefficient heatmap across all models.
library(pheatmap)
get_coefs <- function(feature_mat, surv, n) {
cv_fit <- cv.glmnet(feature_mat, surv, family = "cox", alpha = 0, nfolds = n)
coefs <- as.vector(coef(cv_fit, s = "lambda.min"))
names(coefs) <- rownames(coef(cv_fit, s = "lambda.min"))
coefs
}
set.seed(42)
mat_scaled <- lapply(mat, scale)
coef_list <- lapply(mat, get_coefs, surv = surv, n = n)
# Unify into a variables x models matrix (0 where a feature is absent).
all_vars <- unique(unlist(lapply(coef_list, names)))
coef_mat <- sapply(coef_list, function(coefs) {
out <- setNames(rep(0, length(all_vars)), all_vars)
out[names(coefs)] <- coefs
out
})
pheatmap(coef_mat, cluster_rows = TRUE, cluster_cols = FALSE, scale = "row",
na_col = "grey90", angle_col = 45, main = "Ridge-Cox coefficients by model")
summary_df_pfs <- summary_df
```
## OS
Identical workflow to the PFS block above, using overall survival as the endpoint.
```{r}
outcome.variable_days <- "days_dos"
outcome.variable_status <- "status"
survival_outcome.subset <- survival_outcome %>%
select(!!rlang::sym(outcome.variable_days), !!rlang::sym(outcome.variable_status)) %>% drop_na()
common_samples <- Reduce(intersect, list(
rownames(df.clinicalData),
rownames(df.celltypeProportions),
rownames(df.spatialFeatures_CN),
rownames(df.spatialFeatures_lcross),
rownames(survival_outcome.subset)
))
clin <- df.clinicalData[common_samples, ]
ctp <- df.celltypeProportions[common_samples, ]
cn <- df.spatialFeatures_CN[common_samples, ]
lcross <- df.spatialFeatures_lcross[common_samples, , drop = F]
n <- length(common_samples)
clin_imputed <- clin
for (col in colnames(clin_imputed)) {
if (is.numeric(clin_imputed[[col]])) {
clin_imputed[[col]][is.na(clin_imputed[[col]])] <- median(clin_imputed[[col]], na.rm = TRUE)
} else {
mode_val <- names(sort(table(clin_imputed[[col]]), decreasing = TRUE))[1]
clin_imputed[[col]][is.na(clin_imputed[[col]])] <- mode_val
}
}
clin_mat <- model.matrix(~ . - 1, data = clin_imputed)
lcross <- model.matrix(~ . - 1, data = lcross)
days <- survival_outcome.subset[common_samples, outcome.variable_days]
status <- survival_outcome.subset[common_samples, outcome.variable_status]
surv <- Surv(days, status)
cat(sprintf("Proceeding with %d samples after NA removal.\n", n))
spatial_mat <- cbind(as.matrix(cn), lcross)
mat <- list(
"Clinical" = clin_mat,
"Cell type proportions" = as.matrix(ctp),
"Spatial" = spatial_mat,
"Clinical + Cell type proportions" = cbind(clin_mat, as.matrix(ctp)),
"Clinical + Spatial" = cbind(clin_mat, spatial_mat),
"Cell type proportions + Spatial" = cbind(as.matrix(ctp), spatial_mat),
"Clinical + Cell type proportions + Spatial" = cbind(clin_mat, as.matrix(ctp), spatial_mat)
)
set.seed(42)
fit_model <- function(feature_mat, days, status, surv, n, B = 1000) {
cv_fit <- cv.glmnet(feature_mat, surv, family = "cox", alpha = 0, nfolds = n)
lp <- as.vector(predict(cv_fit, newx = feature_mat, s = "lambda.min"))
std_fit <- coxph(Surv(days, status) ~ lp)
ll <- logLik(std_fit)
aic_val <- AIC(std_fit)
ci_noether <- concordance.index(x = lp, surv.time = days, surv.event = status, method = "noether")
boot_c <- replicate(B, {
idx <- sample(n, replace = TRUE)
tryCatch(concordance.index(x = lp[idx], surv.time = days[idx], surv.event = status[idx])$c.index,
error = function(e) NA, warning = function(w) NA)
})
n_failed <- sum(is.na(boot_c))
if (n_failed > 0) message(sprintf("Bootstrap: %d/%d resamples failed and were skipped.", n_failed, B))
ci_noether$lower <- quantile(boot_c, 0.025, na.rm = TRUE)
ci_noether$upper <- quantile(boot_c, 0.975, na.rm = TRUE)
list(lp = lp, ci = ci_noether, logLik = ll, aic = aic_val)
}
results <- lapply(mat, fit_model, days = days, status = status, surv = surv, n = n, B = 1000)
summary_df <- data.frame(
Model = names(results),
C_index = sapply(results, \(r) round(r$ci$c.index, 3)),
CI_lower = sapply(results, \(r) round(r$ci$lower, 3)),
CI_upper = sapply(results, \(r) round(r$ci$upper, 3)),
AIC = sapply(results, \(r) round(r$aic, 2)),
p_value = sapply(results, \(r) signif(r$ci$p.value, 3))
)
summary_df <- summary_df[order(-summary_df$C_index), ]
print(summary_df, row.names = FALSE)
baseline <- results[["Clinical"]]$ci
cat("\n-- C-index comparison vs. Clinical baseline --\n")
for (nm in setdiff(names(results), "Clinical")) {
comp <- cindex.comp(baseline, results[[nm]]$ci)
cat(sprintf(" Clinical vs %-35s p = %.4f\n", nm, comp$p.value))
}
cat("\n-- Likelihood Ratio Test (LRT) vs. Clinical --\n")
ll_baseline <- results[["Clinical"]]$logLik
for (nm in setdiff(names(results), "Clinical")) {
ll_complex <- results[[nm]]$logLik
lrt_stat <- as.numeric(2 * (ll_complex - ll_baseline))
lrt_p <- pchisq(lrt_stat, df = 1, lower.tail = FALSE)
cat(sprintf(" Clinical vs %-20s | dLL: %6.2f | LRT p = %.4f\n",
nm, (ll_complex - ll_baseline), lrt_p))
}
best_name <- summary_df$Model[1]
if (best_name != "Clinical") {
observed_diff <- results[[best_name]]$ci$c.index - results[["Clinical"]]$ci$c.index
perm_diffs <- replicate(1000, {
idx <- sample(n)
c1 <- concordance.index(results[["Clinical"]]$lp[idx], days, status, method = "noether")$c.index
c2 <- concordance.index(results[[best_name]]$lp[idx], days, status, method = "noether")$c.index
c2 - c1
})
perm_p <- mean(abs(perm_diffs) >= abs(observed_diff))
cat(sprintf("\nPermutation test - Clinical vs %s: dC = %.3f, p = %.4f\n",
best_name, observed_diff, perm_p))
}
summary_df$Model <- factor(summary_df$Model, levels = rev(summary_df$Model))
ggplot(summary_df, aes(x = C_index, y = Model)) +
geom_vline(xintercept = 0.5, linetype = "dashed", color = "grey60") +
geom_errorbarh(aes(xmin = CI_lower, xmax = CI_upper), height = 0.25, color = "grey40") +
geom_point(size = 3, color = "#E64B35") +
geom_vline(xintercept = results[[best_name]]$ci$c.index) +
geom_vline(xintercept = results[["Clinical"]]$ci$c.index) +
scale_x_continuous(limits = c(0.3, 1), breaks = seq(0.3, 1, 0.1)) +
labs(x = "C-index (bootstrap 95% CI)", y = NULL, title = "Survival prediction: model comparison") +
theme_classic(base_size = 13) +
theme(axis.text.y = element_text(hjust = 1))
risk_group <- function(lp) factor(ifelse(lp > median(lp), "High", "Low"))
df_plot <- data.frame(days = days, status = status,
risk_clin = risk_group(results[["Clinical"]]$lp),
risk_best = risk_group(results[[best_name]]$lp))
p1 <- ggsurvplot(survfit(Surv(days, status) ~ risk_clin, data = df_plot),
data = df_plot, pval = TRUE, risk.table = TRUE,
title = "Clinical only", palette = c("#E64B35", "#4DBBD5"))
p2 <- ggsurvplot(survfit(Surv(days, status) ~ risk_best, data = df_plot),
data = df_plot, pval = TRUE, risk.table = TRUE,
title = best_name, palette = c("#E64B35", "#4DBBD5"))
arrange_ggsurvplots(list(p1, p2), ncol = 2)
library(pheatmap)
get_coefs <- function(feature_mat, surv, n) {
cv_fit <- cv.glmnet(feature_mat, surv, family = "cox", alpha = 0, nfolds = n)
coefs <- as.vector(coef(cv_fit, s = "lambda.min"))
names(coefs) <- rownames(coef(cv_fit, s = "lambda.min"))
coefs
}
set.seed(42)
mat_scaled <- lapply(mat, scale)
coef_list <- lapply(mat, get_coefs, surv = surv, n = n)
all_vars <- unique(unlist(lapply(coef_list, names)))
coef_mat <- sapply(coef_list, function(coefs) {
out <- setNames(rep(0, length(all_vars)), all_vars)
out[names(coefs)] <- coefs
out
})
pheatmap(coef_mat, cluster_rows = TRUE, cluster_cols = FALSE, scale = "row",
na_col = "grey90", angle_col = 45, main = "Ridge-Cox coefficients by model")
summary_df_os <- summary_df
```
### Combined plot
```{r, fig.height=8, fig.width = 6.5}
summary_df_os$outcome <- "Overall survival"
summary_df_pfs$outcome <- "Local progression-free survival"
summary_df_all <- rbind(summary_df_os, summary_df_pfs)
# FDR-adjusted p-values and significance flag.
summary_df_all <- summary_df_all %>%
mutate(p_adj = p.adjust(p_value, method = "BH"),
sign = if_else(p_adj < 0.1, "*", ""))
# Helper: colour the feature-block names inside model labels via inline HTML.
colorize_label <- function(label, colors) {
for (category in names(colors)) {
color_val <- colors[category]
replacement <- paste0("<span style='color:", color_val, "'>", category, "</span>")
label <- str_replace_all(label, fixed(category), replacement)
}
label
}
summary_df_all$Model_html <- sapply(summary_df_all$Model, colorize_label, colors = block_colors)
# C-index plots per outcome (ggtext renders the coloured model labels).
p.model.os <- ggplot(summary_df_all %>% filter(outcome == "Overall survival"),
aes(x = C_index, y = reorder(Model_html, C_index))) +
geom_vline(xintercept = 0.5, linetype = "dashed", color = "grey60") +
geom_errorbarh(aes(xmin = CI_lower, xmax = CI_upper), height = 0.25, color = "grey40") +
geom_point(aes(size = -log10(p_value))) +
scale_x_continuous(limits = c(0.3, 1), breaks = seq(0.3, 1, 0.1)) +
labs(x = "C-index (bootstrap 95% CI)", y = NULL, title = "Overall survival") +
theme(axis.text.y = element_markdown(hjust = 1, size = 10, lineheight = 1.2))
p.model.os
p.model.pfs <- ggplot(summary_df_all %>% filter(outcome == "Local progression-free survival"),
aes(x = C_index, y = reorder(Model_html, C_index))) +
geom_vline(xintercept = 0.5, linetype = "dashed", color = "grey60") +
geom_errorbarh(aes(xmin = CI_lower, xmax = CI_upper), height = 0.25, color = "grey40") +
geom_point(aes(size = -log10(p_value))) +
scale_x_continuous(limits = c(0.3, 1), breaks = seq(0.3, 1, 0.1)) +
labs(x = "C-index (bootstrap 95% CI)", y = NULL, title = "Local progression-free survival") +
theme(axis.text.y = element_markdown(hjust = 1, size = 10, lineheight = 1.2))
p.model.os + p.model.pfs + plot_layout(nrow = 2) +
plot_annotation(title = "Survival prediction: Model comparison", tag_levels = "A")
```
# ICI treatment response prediction
To evaluate the predictive value of distinct feature sets for treatment response,
elastic net logistic regression models were trained using the glmnet algorithm as
implemented in the caret R package (v6.0). Feature sets included clinical variables,
cell type proportions, and spatial features (cellular neighborhood composition and
cross-K statistics), each evaluated individually and in combination. Prior to model
training, all features were mean-centered and variance-scaled, and missing clinical
values were imputed using column-wise medians (continuous variables) or modes
(categorical variables). Hyperparameters (alpha and lambda) were selected over a grid
of 10 candidate values. Model performance was assessed using repeated 5-fold
cross-validation (20 repeats), with area under the receiver operating characteristic
curve (AUROC) as the primary metric, computed using the pROC R package. Ninety-five
percent confidence intervals for AUROC values were estimated via DeLong's method, and
pairwise comparisons between feature sets and the clinical baseline model were
performed using DeLong's test with Benjamini-Hochberg correction for multiple testing.
```{r}
# Subset cohort to the ICI-response condition (cond.2).
cohort.subset <- cohort.conditions %>% filter(!is.na(cond.2))
# Remove ici_postop, which would have only one factor level in this subset.
df.clinicalData.subset <- df.clinicalData %>% select(-ici_postop)
# Align all feature blocks to a common set of samples.
common_samples <- Reduce(intersect, list(
rownames(df.clinicalData.subset),
rownames(df.celltypeProportions),
rownames(df.spatialFeatures_CN),
rownames(df.spatialFeatures_lcross),
rownames(cohort.subset)
))
clin <- df.clinicalData.subset[common_samples, ]
ctp <- df.celltypeProportions[common_samples, ]
cn <- df.spatialFeatures_CN[common_samples, ]
lcross <- df.spatialFeatures_lcross[common_samples, , drop = F]
# Binary outcome: responder vs. non-responder.
outcome <- cohort.subset[common_samples, "cond.2"]
outcome <- factor(outcome, levels = c("responder", "non-responder"))
# Impute clinical NAs (median for numeric, mode for categorical) and dummy-encode.
clin_imputed <- clin
for (col in colnames(clin_imputed)) {
if (is.numeric(clin_imputed[[col]])) {
clin_imputed[[col]][is.na(clin_imputed[[col]])] <- median(clin_imputed[[col]], na.rm = TRUE)
} else {
mode_val <- names(sort(table(clin_imputed[[col]]), decreasing = TRUE))[1]
clin_imputed[[col]][is.na(clin_imputed[[col]])] <- mode_val
}
}
clin <- model.matrix(~ . - 1, data = clin_imputed)
lcross <- model.matrix(~ . - 1, data = lcross)
# Individual feature blocks.
X_clin <- clin
X_ctp <- ctp
X_cn <- cn
X_lcross <- lcross
X_spatial <- cbind(cn, lcross)
# Combination blocks.
X_clin_ctp <- cbind(clin, ctp)
X_clin_cn <- cbind(clin, cn)
X_clin_lcross <- cbind(clin, lcross)
X_clin_spatial <- cbind(clin, X_spatial)
X_ctp_spatial <- cbind(ctp, X_spatial)
X_clin_ctp_spatial <- cbind(clin, ctp, X_spatial)
feature_sets <- list(
"Clinical" = X_clin,
"Cell type proportions" = X_ctp,
"Spatial" = X_spatial,
"Clinical + Cell type proportions" = X_clin_ctp,
"Clinical + Spatial" = X_clin_spatial,
"Cell type proportions + Spatial" = X_ctp_spatial,
"Clinical + Cell type proportions + Spatial" = X_clin_ctp_spatial
)
# Cross-validated AUROC via elastic net.
library(MLmetrics)
# Custom CV summary: ROC/Sens/Spec plus PR-AUC.
my_summary <- function(data, lev = NULL, model = NULL) {
roc_stats <- twoClassSummary(data, lev, model)
pr_auc <- PRAUC(data[, lev[1]], ifelse(data$obs == lev[1], 1, 0))
c(roc_stats, PR_AUC = pr_auc)
}
set.seed(42)
cv_ctrl <- trainControl(
method = "repeatedcv",
number = 5,
repeats = 20,
classProbs = TRUE,
summaryFunction = my_summary,
savePredictions = "final"
)
# caret requires syntactically valid class labels.
levels(outcome) <- make.names(levels(outcome))
# Train a glmnet (elastic net) model on one feature block.
fit_model <- function(X, y, ctrl) {
X_scaled <- scale(X) |> as.data.frame() # center + scale
X_scaled[is.na(X_scaled)] <- 0 # impute residual NAs with 0
train(
x = X_scaled,
y = y,
method = "glmnet",
metric = "ROC",
trControl = ctrl,
tuneLength = 10
)
}
results <- map(feature_sets, \(X) fit_model(X, outcome, cv_ctrl))
# Extract AUROC (DeLong CI) and PR-AUC from the held-out CV predictions.
library(PRROC)
extract_metrics <- function(fit, name) {
preds <- fit$pred |>
filter(alpha == fit$bestTune$alpha, lambda == fit$bestTune$lambda)
roc_obj <- roc(preds$obs, preds[, levels(preds$obs)[1]], quiet = TRUE)
ci_roc <- ci.auc(roc_obj)
pos_scores <- preds[preds$obs == levels(preds$obs)[1], levels(preds$obs)[1]]
neg_scores <- preds[preds$obs == levels(preds$obs)[2], levels(preds$obs)[1]]
pr_obj <- pr.curve(scores.class0 = pos_scores, scores.class1 = neg_scores, curve = FALSE)
tibble(model = name,
AUC = as.numeric(auc(roc_obj)),
CI_low = ci_roc[1],
CI_high = ci_roc[3],
PR_AUC = pr_obj$auc.integral)
}
roc_df <- imap_dfr(results, extract_metrics) |> arrange(desc(AUC))
# AUROC bar chart with 95% CI.
roc_df |>
mutate(model = fct_reorder(model, AUC)) |>
ggplot(aes(x = AUC, y = model)) +
geom_col(aes(fill = AUC), width = 0.6, show.legend = FALSE) +
geom_errorbarh(aes(xmin = CI_low, xmax = CI_high), height = 0.25) +
geom_vline(xintercept = 0.5, linetype = "dashed", colour = "red") +
scale_fill_gradient(low = "#91bfdb", high = "#1a6dad") +
scale_x_continuous(limits = c(0, 1), expand = c(0, 0)) +
labs(title = "Predictive performance per feature set",
subtitle = "Elastic net, 5-fold CV x 10 repeats | 95 % CI",
x = "AUROC", y = NULL) +
theme_bw(base_size = 13)
# Overlaid ROC curves for all feature sets.
roc_curves <- imap(results, function(fit, name) {
preds <- fit$pred |>
filter(alpha == fit$bestTune$alpha, lambda == fit$bestTune$lambda)
roc(preds$obs, preds[, levels(outcome)[2]], quiet = TRUE)
})
library(pROC)
roc_plot_df <- imap_dfr(roc_curves, function(r, name) {
data.frame(model = name, specificity = r$specificities, sensitivity = r$sensitivities)
}) |>
mutate(fpr = 1 - specificity,

FinalFigures_public.Rmd at commit 86bf55d, no license · at the source

Overview

Authors: Stefanos Voglis^1,2,3, Daniel Schulz^1,2, Nils Eling^1,2, Natalie De Souza^1,2, Luca Regli³, Marian Christoph Neidert⁴, Marcus Czabanka⁵, Michael Weller⁶, Emilie Le Rhun⁷, Daniela Mihic-Probst⁸, Mitchell Levesque⁹, Bernd Bodenmiller^1,2

ORCID iDs: Stefanos Voglis, Daniel Schulz, Nils Eling, Natalie De Souza, Luca Regli, Marian Christoph Neidert, Michael Weller, Emilie Le Rhun, Mitchell Levesque

Institutions: University of Zurich (Switzerland); ETH Zurich (Switzerland); University Hospital Zurich (Switzerland); Institute of Molecular Health Sciences (Switzerland); Kantonsspital St. Gallen (Switzerland); Goethe University Frankfurt (Germany); University Hospital Frankfurt (Germany)

Journal: Neuro-oncology, volume 28, issue 9, pages 2211-2223

Dates: received 20 January 2026; accepted 26 May 2026; published online 28 May 2026; in print September 2026

Type: Research article · Language: English

License: CC BY-NC

Identifiers: DOI 10.1093/neuonc/noag128 · PMID 42209753 · PMCID PMC13550667 · OpenAlex W7162786992

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: genetics / omics (modality), human (organism), other condition (population), clinical / translational (subfield)

Methods: Statistics, Machine learning, Preprocessing, fMRI & imaging

Keywords: melanoma brain metastases, precision oncology, predictive biomarkers, proteomics, tumor microenvironment

MeSH: Brain Neoplasms*, Immune Checkpoint Inhibitors*, Melanoma*, Single-Cell Analysis*, Tumor Microenvironment*, Adult, Aged, Biomarkers, Tumor, Female, Humans, Male, Middle Aged, Prognosis (* major topic)

Topic: Cancer Immunotherapy and Biomarkers (Oncology, Medicine), according to OpenAlex

Funding: S.V. by the Swiss Academy of Medical Sciences/Gottfried and Julia Bangerter-Rhyner-Foundation (YTCR 58/20); the Swiss Academy of Medical Sciences/Gottfried; the Brihaye EANS Research Fund; Julia Bangerter-Rhyner-Foundation (58/20); “Stiftung Tumorforschung Kopf-Hals.”

Citations: not cited yet (Europe PMC); 44 references in the paper

Abstract

Background: A subset of patients with advanced-stage melanoma with brain metastases responds intracranially to immune-checkpoint inhibitors (ICIs). We reasoned that features of the spatial architecture of the tumor microenvironment correlate with treatment responses.

Methods: We used highly multiplexed single-cell imaging to characterize the TME in 44 samples of MBM patients; both treatment naïve subjects and subjects who had received ICI were included in the cohort. Samples were stained with a 42-plex metal-isotope tagged antibody panel and subsequently imaged using imaging mass cytometry. Downstream analysis focused on phenotypic and spatial characteristics of the tumor microenvironments of ICI responders and non-responders and on correlations of protein expression with overall and local progression-free survival.

Results: Single-cell phenotyping identified more than 1.1 million cells, predominantly tumor (67%) and immune cells (29%). ICI responders exhibited higher baseline proportions of B cells, IDO+ macrophages, and specific T cell subsets (regulatory and exhausted/activated CD8+). Conversely, neutrophil infiltration was negatively associated with intracranial response and survival. Spatial analysis revealed that responders possess a more localized, clustered immune cell pattern compared to the dispersed patterns seen in non-responders. Modeling CD8+ T cell distribution showed that localized infiltration was positively associated with overall and local progression-free survival, while dispersed patterns correlated with poorer outcomes. Furthermore, ICI-pretreated samples showed higher immune-tumor co-localization than treatment-naïve samples.

Conclusions: Our single-cell, spatially resolved characterization of the tumor microenvironment of human melanoma brain metastases demonstrated distinct phenotypic and spatial infiltration patterns of various immune cells that are associated with ICI response and survival outcomes.

Reproduced under the paper's license (CC BY-NC), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 18 matches between paragraphs and lines of code.

BodenmillerGroup/steinbock

License: MIT

State: the link answers, verified on 26 September 2026

Evidence: files inventoried

Commit: 848fc5ed547c5c98f25d69b29b75de9be9824f3c, 27 March 2026

Languages: Python (88), Shell (1)

Size: 153 files, 89 scripts

Software Heritage: not archived

Found in: the text, “Single-Cell Segmentation, Data Processing and Qu”

Holds: README, license file, environment (docker-compose.yml, Dockerfile, pyproject.toml, requirements-deepcell.txt, requirements-napari.txt, requirements.txt, requirements_devel.txt, requirements_docs.txt, requirements_test.txt, setup.cfg), tests, continuous integration, documentation

Not found: CITATION.cff

Tools: NumPy (27 files), pandas (13 files), scikit-image (4 files), SciPy (4 files), Keras (3 files), TensorFlow (3 files), tifffile (3 files), anndata (2 files), NetworkX (2 files), Cellpose (1 file), h5py (1 file), imageio (1 file), napari (1 file), OpenCV (1 file)

Availability: 1 check, the latest on 26 September 2026: the link answers

26 September 2026: the link answers

91 files

conftest.py, Python, 31 lines
entrypoint.sh, Shell, 3 lines
steinbock/__init__.py, Python, 3 lines
steinbock/__main__.py, Python, 17 lines
steinbock/_cli/__init__.py, Python, 3 lines
steinbock/_cli/_cli.py, Python, 35 lines
steinbock/_cli/apps.py, Python, 150 lines
steinbock/_cli/utils.py, Python, 38 lines
steinbock/_cli/visualization.py, Python, 64 lines
steinbock/_env.py, Python, 44 lines
steinbock/_steinbock.py, Python, 7 lines
steinbock/classification/__init__.py, Python, 3 lines
steinbock/classification/_classification.py, Python, 5 lines
steinbock/classification/_cli/__init__.py, Python, 3 lines
steinbock/classification/_cli/_cli.py, Python, 16 lines
steinbock/classification/_cli/ilastik.py, Python, 343 lines
steinbock/classification/ilastik/__init__.py, Python, 37 lines
steinbock/classification/ilastik/_ilastik.py, Python, 573 lines
steinbock/classification/ilastik/data/__init__.py, Python, 1 line
steinbock/export/__init__.py, Python, 1 line
steinbock/export/_cli/__init__.py, Python, 3 lines
steinbock/export/_cli/_cli.py, Python, 159 lines
steinbock/export/_cli/data.py, Python, 284 lines
steinbock/export/_cli/graphs.py, Python, 72 lines
steinbock/export/data.py, Python, 136 lines
steinbock/export/graphs.py, Python, 50 lines
steinbock/io.py, Python, 379 lines
steinbock/measurement/__init__.py, Python, 3 lines
steinbock/measurement/_cli/__init__.py, Python, 3 lines
steinbock/measurement/_cli/_cli.py, Python, 22 lines
steinbock/measurement/_cli/cellprofiler.py, Python, 182 lines
steinbock/measurement/_cli/intensities.py, Python, 98 lines
steinbock/measurement/_cli/neighbors.py, Python, 89 lines
steinbock/measurement/_cli/regionprops.py, Python, 68 lines
steinbock/measurement/_measurement.py, Python, 5 lines
steinbock/measurement/cellprofiler/__init__.py, Python, 3 lines
steinbock/measurement/cellprofiler/_cellprofiler.py, Python, 43 lines
steinbock/measurement/cellprofiler/data/__init__.py, Python, 1 line
steinbock/measurement/intensities.py, Python, 66 lines
steinbock/measurement/neighbors.py, Python, 200 lines, 1 match
steinbock/measurement/regionprops.py, Python, 52 lines
steinbock/preprocessing/__init__.py, Python, 3 lines
steinbock/preprocessing/_cli/__init__.py, Python, 3 lines
steinbock/preprocessing/_cli/_cli.py, Python, 20 lines
steinbock/preprocessing/_cli/external.py, Python, 142 lines
steinbock/preprocessing/_cli/imc.py, Python, 256 lines
steinbock/preprocessing/_preprocessing.py, Python, 5 lines
steinbock/preprocessing/external.py, Python, 87 lines
steinbock/preprocessing/imc.py, Python, 493 lines
steinbock/segmentation/__init__.py, Python, 3 lines
steinbock/segmentation/_cli/__init__.py, Python, 3 lines
steinbock/segmentation/_cli/_cli.py, Python, 22 lines
steinbock/segmentation/_cli/cellpose.py, Python, 253 lines
steinbock/segmentation/_cli/cellprofiler.py, Python, 115 lines
steinbock/segmentation/_cli/deepcell.py, Python, 230 lines
steinbock/segmentation/_segmentation.py, Python, 5 lines
steinbock/segmentation/cellpose.py, Python, 122 lines
steinbock/segmentation/cellprofiler/__init__.py, Python, 3 lines
steinbock/segmentation/cellprofiler/_cellprofiler.py, Python, 43 lines
steinbock/segmentation/cellprofiler/data/__init__.py, Python, 1 line
steinbock/segmentation/deepcell.py, Python, 138 lines
steinbock/utils/__init__.py, Python, 3 lines
steinbock/utils/_cli/__init__.py, Python, 3 lines
steinbock/utils/_cli/_cli.py, Python, 16 lines
steinbock/utils/_cli/expansion.py, Python, 46 lines
steinbock/utils/_cli/matching.py, Python, 47 lines
steinbock/utils/_cli/mosaics.py, Python, 98 lines
steinbock/utils/_utils.py, Python, 5 lines
steinbock/utils/expansion.py, Python, 25 lines
steinbock/utils/matching.py, Python, 50 lines
steinbock/utils/mosaics.py, Python, 125 lines
steinbock/visualization.py, Python, 42 lines
tests/classification/test_ilastik_classification.py, Python, 180 lines
tests/export/test_data_export.py, Python, 48 lines
tests/export/test_graphs_export.py, Python, 43 lines
tests/measurement/test_cellprofiler_measurement.py, Python, 32 lines
tests/measurement/test_intensities_measurement.py, Python, 58 lines
tests/measurement/test_neighbors_measurement.py, Python, 68 lines
tests/measurement/test_regionprops_measurement.py, Python, 47 lines
tests/preprocessing/test_external_preprocessing.py, Python, 9 lines
tests/preprocessing/test_imc_preprocessing.py, Python, 97 lines
tests/segmentation/test_cellpose_segmentation.py, Python, 15 lines
tests/segmentation/test_cellprofiler_segmentation.py, Python, 32 lines
tests/segmentation/test_deepcell_segmentation.py, Python, 35 lines
tests/test_io.py, Python, 61 lines
tests/test_viewer.py, Python, 6 lines
tests/utils/test_expansion_utils.py, Python, 9 lines
tests/utils/test_matching_utils.py, Python, 36 lines
tests/utils/test_mosaics_utils.py, Python, 40 lines
LICENSE, License, 21 lines
README.md, Text, 31 lines

At the source: github.com/BodenmillerGroup/steinbock

StefanosVoglis/Voglis_et_al_2026_Neuro-Oncology

License: none: the authors keep all their rights

State: the link answers, verified on 26 September 2026

Evidence: files inventoried

Commit: 86bf55d7d4faebe861966172590a76694892d500, 11 June 2026

Languages: R (1)

Size: 1 file, 1 script

Software Heritage: not archived

Found in: “Data Availability”

Holds: 1 notebook

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: caret (1 file), circlize (1 file), ComplexHeatmap (1 file), cowplot (1 file), data.table (1 file), easystats (1 file), edgeR (1 file), ggplot2 (1 file), glmnet (1 file), lme4 (1 file), patchwork (1 file), pheatmap (1 file), pROC (1 file), reshape2 (1 file), rstatix (1 file), survival (1 file), tidyverse (1 file)

Availability: 1 check, the latest on 26 September 2026: the link answers

26 September 2026: the link answers

1 file

FinalFigures_public.Rmd, R, 4,281 lines, 17 matches

At the source: github.com/StefanosVoglis/Voglis_et_al_2026_Neuro-Oncology

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
90 scripts, each with its path and the digest of its content;
18 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

zenodo:20539274, at Zenodo; found in DataCite
zenodo:20539275, at Zenodo; found in “Data Availability”

Data Availability

Code used for data analysis and figure/table generation will be made publicly available upon publication via a GitHub repository (https://github.com/StefanosVoglis/Voglis_et_al_2026_Neuro-Oncology). Processed data is available under https://doi.org/10.5281/zenodo.20539275.

Reproduced under the paper's license (CC BY-NC), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 12 authors, 5 keywords, 13 MeSH terms, 5 funders, 42 references.

Cite

This paper

Voglis, S., Schulz, D., Eling, N., De Souza, N., Regli, L., Neidert, M. C., Czabanka, M., Weller, M., Le Rhun, E., Mihic-Probst, D., Levesque, M., & Bodenmiller, B. (2026). Spatially-resolved single-cell imaging of melanoma brain metastases identifies localized immune patterns predictive of immune checkpoint blockade response. Neuro-oncology, 28(9), 2211-2223. https://doi.org/10.1093/neuonc/noag128

BibTeX

@article{voglis2026spatially, author = {Voglis, Stefanos and Schulz, Daniel and Eling, Nils and De Souza, Natalie and Regli, Luca and Neidert, Marian Christoph and Czabanka, Marcus and Weller, Michael and Le Rhun, Emilie and Mihic-Probst, Daniela and Levesque, Mitchell and Bodenmiller, Bernd}, title = {{Spatially-resolved single-cell imaging of melanoma brain metastases identifies localized immune patterns predictive of immune checkpoint blockade response}}, journal = {Neuro-oncology}, year = {2026}, month = sep, volume = {28}, number = {9}, pages = {2211--2223}, publisher = {Oxford University Press}, issn = {1522-8517}, doi = {10.1093/neuonc/noag128}, url = {https://doi.org/10.1093/neuonc/noag128}, pmid = {42209753}, pmcid = {PMC13550667} }

RIS

TY - JOUR AU - Voglis, Stefanos AU - Schulz, Daniel AU - Eling, Nils AU - De Souza, Natalie AU - Regli, Luca AU - Neidert, Marian Christoph AU - Czabanka, Marcus AU - Weller, Michael AU - Le Rhun, Emilie AU - Mihic-Probst, Daniela AU - Levesque, Mitchell AU - Bodenmiller, Bernd TI - Spatially-resolved single-cell imaging of melanoma brain metastases identifies localized immune patterns predictive of immune checkpoint blockade response T2 - Neuro-oncology J2 - Neuro Oncol PY - 2026 DA - 2026/09/01 VL - 28 IS - 9 SP - 2211 EP - 2223 SN - 1522-8517 PB - Oxford University Press DO - 10.1093/neuonc/noag128 UR - https://doi.org/10.1093/neuonc/noag128 LA - en ER -

CSL-JSON

{ "id": "10.1093/neuonc/noag128", "type": "article-journal", "title": "Spatially-resolved single-cell imaging of melanoma brain metastases identifies localized immune patterns predictive of immune checkpoint blockade response", "container-title": "Neuro-oncology", "author": [ { "family": "Voglis", "given": "Stefanos" }, { "family": "Schulz", "given": "Daniel" }, { "family": "Eling", "given": "Nils" }, { "family": "De Souza", "given": "Natalie" }, { "family": "Regli", "given": "Luca" }, { "family": "Neidert", "given": "Marian Christoph" }, { "family": "Czabanka", "given": "Marcus" }, { "family": "Weller", "given": "Michael" }, { "family": "Le Rhun", "given": "Emilie" }, { "family": "Mihic-Probst", "given": "Daniela" }, { "family": "Levesque", "given": "Mitchell" }, { "family": "Bodenmiller", "given": "Bernd" } ], "container-title-short": "Neuro Oncol", "volume": "28", "issue": "9", "page": "2211-2223", "DOI": "10.1093/neuonc/noag128", "PMID": "42209753", "PMCID": "PMC13550667", "ISSN": "1522-8517", "publisher": "Oxford University Press", "URL": "https://doi.org/10.1093/neuonc/noag128", "language": "en", "issued": { "date-parts": [ [ 2026, 9, 1 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: pROC, imageio, edgeR, 21 other tools, genetics / omics, other condition, 2 references
[2] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: glmnet, survival, caret, 18 other tools
[3] doi:10.1093/bioinformatics/btag592 [code]: Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: pROC, survival, edgeR, 15 other tools, genetics / omics, other condition
[4] doi:10.1016/j.isci.2026.115657 [code]: Integration of machine learning to develop a disulfidptosis model for predicting glioma prognosis, immunotherapy response, and drug.
Journal: iScience
In common: glmnet, pROC, survival, 13 other tools, clinical / translational, other condition
[5] doi:10.1016/j.xcrm.2026.102682 [code]: TET CpG sequence-context-specific DNA demethylation shapes progression of IDH-mutant gliomas.
Journal: Cell reports. Medicine
In common: glmnet, pROC, survival, 13 other tools, genetics / omics, other condition
[6] doi:10.3390/ijms27104466 [code]: Uncovering the Key Circuit FOSL2/FOS/EGR3/EGR1, Contributing to the Hyperexcitability of Excitatory Neurons in the Epileptic Temporal Cortex and Hippocampus.
Journal: International journal of molecular sciences
In common: glmnet, pROC, caret, 12 other tools, genetics / omics
[7] doi:10.1002/imt2.70163 [code]: Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: edgeR, rstatix, anndata, 14 other tools, genetics / omics
[8] doi:10.1073/pnas.2609132123 [code]: A human lysosomal storage disorder toolkit for decoding proteome landscapes in cortical-like and dopaminergic-like induced neurons.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: Cellpose, tifffile, rstatix, 13 other tools, genetics / omics
[9] doi:10.1371/journal.pcbi.1014571 [code]: SynAPSeg: A novel dataset and image analysis framework for deep learning-based synapse detection and quantification.
Journal: PLoS computational biology
In common: Cellpose, napari, imageio, 9 other tools, 1 reference
[10] doi:10.1016/j.isci.2026.116439 [code]: Decoding the role of transcriptomic clocks in the human prefrontal cortex.
Journal: iScience
In common: glmnet, Keras, circlize, 12 other tools, genetics / omics

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.