ADC Visium and Xu Atlas GE1-GE10 Analysis

Cancer epithelial program states before and after antibody-drug conjugate treatment

knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE,
  fig.width = 9,
  fig.height = 6,
  dpi = 150
)

options(stringsAsFactors = FALSE)
set.seed(123)

knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(root.dir = normalizePath("/work/InternalMedicine/s239947/ADC"))
getwd()
knitr::opts_chunk$get("root.dir")

Purpose and analysis design

This notebook compares the 10 Xu cancer epithelial gene elements (GE1-GE10) across ADC-treated breast cancer biopsies. It is designed for the currently known cohort of 17 biopsies (6 pre-ADC and 11 post-ADC), including three matched sample pairs:

• ADC005 and ADC008
• ADC009 and ADC010
• ADC017 and ADC018

The direction of each pair is taken from clinical metadata; it is not inferred from the sample number.

The principal questions are:

1. Which Xu GE states are represented in the ADC Visium tumor regions?
2. Are any GE scores or GE-state proportions higher post-ADC than pre-ADC?
3. Which GE profiles appear stable within the three matched pre/post pairs?
4. Do post-ADC samples contain residual states that are uncommon in the Xu primary-breast reference?
5. Once response metadata are available, which pretreatment or residual GE states are associated with response?

Important interpretation

Visium observations are spots, not individual cells. Therefore, this notebook clusters cancer-epithelial-enriched spots by their GE profiles. It does not claim that each Visium spot is a single cancer epithelial cell. If a cancer-cell fraction from RCTD, cell2location, Spotiphy, or another deconvolution method is available, use it as the primary tumor-spot filter below. Actual single-cell validation should be performed in matched Xenium or scRNA-seq data.

The integrated object contains zeros for every GE in non-cancer-epithelial spots. Those spots must be removed before scaling, clustering, or pre/post comparison. Otherwise, the shared zero vector will create an artificial low-GE cluster.

Statistical unit

Spots are used for spatial visualization and state discovery. The biopsy/sample is the statistical unit for pre/post and response analyses. This avoids pseudoreplication from treating thousands of spots as thousands of patients. The three matched pairs are summarized as within-patient changes; the other samples support a secondary, exploratory cross-sectional comparison.

Xu method reproduced here

Xu et al. calculated UCell scores for each of the 10 GEs, standardized each GE across cancer epithelial cells, and assigned each cell to the GE with the highest standardized score. This notebook reproduces that dominant-GE labeling and also performs secondary unsupervised clustering of the full 10-dimensional profile.

Primary references:

• Xu et al., Cell Reports Medicine 2024
• Chan Lab BreastCancer_Integrated code
• Figure 3 analysis code

1. Packages

This notebook does not install packages automatically. Install any missing packages once, then rerun this section.

required_packages <- c(
  "Seurat",
  "dplyr",
  "tidyr",
  "purrr",
  "tibble",
  "readr",
  "stringr",
  "forcats",
  "ggplot2",
  "patchwork",
  "pheatmap",
  "cluster",
  "scales"
)

missing_packages <- setdiff(required_packages, rownames(installed.packages()))

if (length(missing_packages) > 0) {
  stop(
    "Install these packages before continuing: ",
    paste(missing_packages, collapse = ", ")
  )
}

invisible(lapply(required_packages, library, character.only = TRUE))

2. Configuration

Edit this section first. The Xu object is optional for the ADC-only analysis, but required for reference-anchored comparison and out-of-reference screening.

The GE column vectors may be left as NULL for automatic detection. If automatic detection reports ambiguity, provide the exact 10 columns in GE1 to GE10 order, for example paste0("raw_GE", 1:10).

CFG <- list(
  # Input files ---------------------------------------------------------------
  adc_rds = "/project/InternalMedicine/Chan_lab/shared/shared_spatial/integrated_ADC_objects/20260804_integrated_ADC_Visium_object.rds",
  xu_rds = "/work/InternalMedicine/s239947/BMBC/objects/xu_mapping/xu_epithelial_reference_GE_scored_v2.rds",
  clinical_csv = "/work/InternalMedicine/s239947/ADC/adc_clinical_metadata.csv",

  # Output folder -------------------------------------------------------------
  output_dir = "ADC_GE_outputs",

  # ADC Visium metadata -------------------------------------------------------
  sample_col = "sample",
  patient_col = "patient_id",
  timepoint_col = "timepoint",
  adc_drug_col = "adc_drug",
  metastatic_site_col = NULL,
  response_group_col = "response_group",
  batch_col = NULL,

  # Exact columns in GE1-GE10 order, or NULL for automatic detection.
  adc_ge_cols = paste0("raw_GE", 1:10),

  # Preferred tumor-spot filter. Use a deconvolved cancer epithelial fraction
  # when available. Otherwise specify an epithelial annotation column. If both
  # are NULL, the fallback is any nonzero GE score.
  tumor_fraction_col = NULL,
  min_tumor_fraction = 0.20,
  epithelial_col = NULL,
  cancer_epithelial_values = c(
    "Cancer epithelial",
    "Cancer Epithelial",
    "Malignant epithelial",
    "Tumor",
    "Epithelial"
  ),

  # Xu reference metadata ----------------------------------------------------
  # Set xu_object_is_cancer_epithelial = TRUE only if the object is already a
  # cancer epithelial subset. The deposited Xu metadata does not necessarily
  # contain a final cancer/CNV call.
  xu_object_is_cancer_epithelial = TRUE,
  xu_epithelial_col = NULL,
  xu_cancer_values = c("Cancer epithelial", "Cancer Epithelial", "Tumor"),
  xu_sample_col = "Sample",
  xu_ge_cols = stats::setNames(
  paste0("GE", 1:10, "_UCell"),
  paste0("GE", 1:10)),

  # Clustering and display ----------------------------------------------------
  max_spots_per_sample_for_clustering = 500,
  max_xu_cells_per_ge_for_display = 500,
  k_ge_clusters = NULL,       # NULL selects k from the silhouette screen.
  k_range = 2:8,
  silhouette_max_n = 2000,
  out_of_reference_quantile = 0.95,

  # A descriptive stability screen, not a formal equivalence test.
  conservation_margin_z = 0.25,

  # Saving the full Seurat object can be large; leave FALSE until needed.
  save_augmented_rds = FALSE
)

GE_NAMES <- paste0("GE", 1:10)
GE_RAW_COLS <- paste0("raw_", GE_NAMES)
GE_Z_ADC_COLS <- paste0(GE_NAMES, "_z_adc")
GE_Z_REF_COLS <- paste0(GE_NAMES, "_z_ref")

dir.create(CFG$output_dir, recursive = TRUE, showWarnings = FALSE)

3. Helper functions

is_provided <- function(x) {
  length(x) == 1 && !is.na(x) && nzchar(x)
}

extract_seurat_metadata <- function(object) {
  if (!inherits(object, "Seurat")) {
    stop("Expected a Seurat object, but received: ", paste(class(object), collapse = ", "))
  }

  object[[]] |>
    tibble::rownames_to_column("spot_id")
}

resolve_ge_columns <- function(metadata, explicit = NULL, object_name = "object") {
  if (!is.null(explicit)) {
    if (length(explicit) != 10) {
      stop(object_name, ": explicit GE column vector must contain exactly 10 columns.")
    }

    if (is.null(names(explicit))) {
      names(explicit) <- GE_NAMES
    }

    explicit <- explicit[GE_NAMES]
    missing <- setdiff(unname(explicit), colnames(metadata))

    if (length(missing) > 0) {
      stop(object_name, ": missing configured GE columns: ", paste(missing, collapse = ", "))
    }

    return(explicit)
  }

  resolved <- setNames(rep(NA_character_, 10), GE_NAMES)

  for (i in seq_along(GE_NAMES)) {
    patterns <- c(
      sprintf("^raw[._-]?(GE|X)[._-]?%d([._-]?UCell)?$", i),
      sprintf("^GE[._-]?%d[._-]?(UCell|score)$", i),
      sprintf("^GE%d$", i)
    )

    hits <- unique(unlist(lapply(
      patterns,
      function(pattern) grep(pattern, colnames(metadata), value = TRUE, ignore.case = TRUE)
    )))

    hits <- hits[!grepl("z|scaled|label|cluster", hits, ignore.case = TRUE)]

    if (length(hits) != 1) {
      stop(
        object_name, ": could not uniquely resolve ", GE_NAMES[i], ". Found: ",
        ifelse(length(hits) == 0, "none", paste(hits, collapse = ", ")),
        ". Set the exact GE columns in CFG."
      )
    }

    resolved[i] <- hits
  }

  resolved
}

add_canonical_ge_scores <- function(metadata, ge_map, object_name = "object") {
  for (ge in names(ge_map)) {
    source_col <- ge_map[[ge]]
    target_col <- paste0(ge, "_raw")
    source_values <- metadata[[source_col]]

    if (is.factor(source_values)) {
      source_values <- as.character(source_values)
    }

    metadata[[target_col]] <- suppressWarnings(as.numeric(source_values))
  }

  if (any(!is.finite(as.matrix(metadata[, GE_RAW_COLS])) &
          !is.na(as.matrix(metadata[, GE_RAW_COLS])))) {
    stop(object_name, ": one or more GE values are infinite.")
  }

  metadata
}

derive_tumor_flag <- function(metadata, cfg) {
  ge_matrix <- as.matrix(metadata[, GE_RAW_COLS, drop = FALSE])
  score_flag <- rowSums(abs(ge_matrix), na.rm = TRUE) > 0

  if (is_provided(cfg$tumor_fraction_col)) {
    if (!cfg$tumor_fraction_col %in% colnames(metadata)) {
      stop("Configured tumor_fraction_col is absent: ", cfg$tumor_fraction_col)
    }

    primary_flag <- as.numeric(metadata[[cfg$tumor_fraction_col]]) >= cfg$min_tumor_fraction
    criterion <- paste0(
      cfg$tumor_fraction_col,
      " >= ",
      cfg$min_tumor_fraction
    )
  } else if (is_provided(cfg$epithelial_col)) {
    if (!cfg$epithelial_col %in% colnames(metadata)) {
      stop("Configured epithelial_col is absent: ", cfg$epithelial_col)
    }

    primary_flag <- as.character(metadata[[cfg$epithelial_col]]) %in%
      cfg$cancer_epithelial_values
    criterion <- paste0("annotation in ", cfg$epithelial_col)
  } else {
    primary_flag <- score_flag
    criterion <- "fallback: at least one GE score is nonzero"
  }

  primary_flag[is.na(primary_flag)] <- FALSE
  final_flag <- primary_flag & score_flag

  list(
    flag = final_flag,
    primary_flag = primary_flag,
    score_flag = score_flag,
    criterion = criterion
  )
}

sample_balanced_moments <- function(score_matrix, sample_id) {
  sample_id <- as.character(sample_id)
  sample_sizes <- table(sample_id)
  weights <- 1 / as.numeric(sample_sizes[sample_id])

  weighted_mean <- vapply(seq_len(ncol(score_matrix)), function(j) {
    keep <- is.finite(score_matrix[, j])
    sum(weights[keep] * score_matrix[keep, j]) / sum(weights[keep])
  }, numeric(1))

  weighted_sd <- vapply(seq_len(ncol(score_matrix)), function(j) {
    keep <- is.finite(score_matrix[, j])
    centered <- score_matrix[keep, j] - weighted_mean[j]
    sqrt(sum(weights[keep] * centered^2) / sum(weights[keep]))
  }, numeric(1))

  names(weighted_mean) <- colnames(score_matrix)
  names(weighted_sd) <- colnames(score_matrix)

  if (any(!is.finite(weighted_sd)) || any(weighted_sd == 0)) {
    stop("At least one GE has zero or undefined variance after tumor filtering.")
  }

  list(mean = weighted_mean, sd = weighted_sd)
}

ordinary_moments <- function(score_matrix) {
  mu <- colMeans(score_matrix, na.rm = TRUE)
  sigma <- apply(score_matrix, 2, stats::sd, na.rm = TRUE)

  if (any(!is.finite(sigma)) || any(sigma == 0)) {
    stop("At least one Xu GE has zero or undefined variance.")
  }

  list(mean = mu, sd = sigma)
}

apply_z_scaling <- function(metadata, raw_cols, moments, suffix) {
  score_matrix <- as.matrix(metadata[, raw_cols, drop = FALSE])
  z_matrix <- sweep(score_matrix, 2, moments$mean, FUN = "-")
  z_matrix <- sweep(z_matrix, 2, moments$sd, FUN = "/")
  colnames(z_matrix) <- paste0(GE_NAMES, suffix)

  metadata |>
    dplyr::select(-dplyr::any_of(colnames(z_matrix))) |>
    dplyr::bind_cols(as.data.frame(z_matrix))
}

add_dominant_ge <- function(metadata, z_cols, label_col, margin_col) {
  z_matrix <- as.matrix(metadata[, z_cols, drop = FALSE])
  usable <- rowSums(is.finite(z_matrix)) == ncol(z_matrix)

  label <- rep(NA_character_, nrow(metadata))
  margin <- rep(NA_real_, nrow(metadata))
  maximum <- rep(NA_real_, nrow(metadata))

  if (any(usable)) {
    z_use <- z_matrix[usable, , drop = FALSE]
    winner <- max.col(z_use, ties.method = "first")
    ordered <- t(apply(z_use, 1, sort, decreasing = TRUE))

    label[usable] <- GE_NAMES[winner]
    margin[usable] <- ordered[, 1] - ordered[, 2]
    maximum[usable] <- ordered[, 1]
  }

  metadata[[label_col]] <- factor(label, levels = GE_NAMES)
  metadata[[margin_col]] <- margin
  metadata[[paste0(label_col, "_max_z")]] <- maximum
  metadata
}

safe_sample_rows <- function(data, n) {
  if (nrow(data) <= n) {
    return(data)
  }

  data[sample.int(nrow(data), size = n, replace = FALSE), , drop = FALSE]
}

balanced_sample <- function(data, group_col, max_per_group) {
  split(data, data[[group_col]], drop = TRUE) |>
    lapply(safe_sample_rows, n = max_per_group) |>
    dplyr::bind_rows()
}

nearest_centroid <- function(score_matrix, centroids) {
  x2 <- rowSums(score_matrix^2)
  c2 <- rowSums(centroids^2)
  distance_squared <- outer(x2, c2, "+") - 2 * score_matrix %*% t(centroids)
  distance_squared[distance_squared < 0] <- 0

  nearest_index <- max.col(-distance_squared, ties.method = "first")

  tibble::tibble(
    nearest_label = rownames(centroids)[nearest_index],
    distance = sqrt(distance_squared[cbind(seq_len(nrow(score_matrix)), nearest_index)])
  )
}

add_metadata_to_seurat <- function(object, metadata, columns) {
  metadata_index <- match(colnames(object), metadata$spot_id)

  if (any(is.na(metadata_index))) {
    stop("Could not match all Seurat barcodes back to the analysis metadata.")
  }

  for (column in columns) {
    object[[column]] <- metadata[[column]][metadata_index]
  }

  object
}

normalize_timepoint <- function(x) {
  x <- stringr::str_to_lower(stringr::str_trim(as.character(x)))

  dplyr::case_when(
    stringr::str_detect(x, "^(pre|baseline|before)") ~ "Pre",
    stringr::str_detect(x, "^(post|after|on[- ]?treatment|residual)") ~ "Post",
    TRUE ~ NA_character_
  )
}

permutation_mean_difference <- function(
    values,
    groups,
    max_exact = 50000,
    n_permutations = 20000,
    seed = 123) {
  keep <- is.finite(values) & groups %in% c("Pre", "Post")
  values <- values[keep]
  groups <- groups[keep]

  n_pre <- sum(groups == "Pre")
  n_post <- sum(groups == "Post")

  if (n_pre < 2 || n_post < 2) {
    return(list(effect = NA_real_, p_value = NA_real_, method = NA_character_))
  }

  observed <- mean(values[groups == "Post"]) - mean(values[groups == "Pre"])
  n_total <- length(values)
  number_assignments <- choose(n_total, n_pre)

  difference_for_pre <- function(pre_index) {
    post_index <- setdiff(seq_len(n_total), pre_index)
    mean(values[post_index]) - mean(values[pre_index])
  }

  if (number_assignments <= max_exact) {
    assignments <- combn(seq_len(n_total), n_pre)
    null_effects <- apply(assignments, 2, difference_for_pre)
    p_value <- mean(abs(null_effects) >= abs(observed) - sqrt(.Machine$double.eps))
    method <- "exact label permutation"
  } else {
    set.seed(seed)
    null_effects <- replicate(
      n_permutations,
      difference_for_pre(sample.int(n_total, n_pre, replace = FALSE))
    )
    p_value <- (1 + sum(abs(null_effects) >= abs(observed))) /
      (n_permutations + 1)
    method <- paste0(n_permutations, " label permutations")
  }

  list(effect = observed, p_value = p_value, method = method)
}

cliffs_delta <- function(post, pre) {
  post <- post[is.finite(post)]
  pre <- pre[is.finite(pre)]

  if (length(post) == 0 || length(pre) == 0) {
    return(NA_real_)
  }

  comparisons <- outer(post, pre, FUN = "-")
  mean(comparisons > 0) - mean(comparisons < 0)
}

top_ge_jaccard <- function(pre_values, post_values, n_top = 3) {
  top_pre <- names(sort(pre_values, decreasing = TRUE))[seq_len(n_top)]
  top_post <- names(sort(post_values, decreasing = TRUE))[seq_len(n_top)]
  length(intersect(top_pre, top_post)) / length(union(top_pre, top_post))
}

4. Load the ADC Visium object and audit GE metadata

if (!file.exists(CFG$adc_rds)) {
  stop("Edit CFG$adc_rds so it points to the integrated ADC Visium RDS file.")
}

adc <- readRDS(CFG$adc_rds)
adc_metadata <- extract_seurat_metadata(adc)

if (!CFG$sample_col %in% colnames(adc_metadata)) {
  stop("ADC sample column is absent: ", CFG$sample_col)
}

adc_metadata$sample_id <- as.character(adc_metadata[[CFG$sample_col]])

if (any(is.na(adc_metadata$sample_id) | !nzchar(adc_metadata$sample_id))) {
  stop("At least one ADC spot has a missing sample identifier.")
}

adc_ge_map <- resolve_ge_columns(
  metadata = adc_metadata,
  explicit = CFG$adc_ge_cols,
  object_name = "ADC Visium object"
)

adc_metadata <- add_canonical_ge_scores(
  metadata = adc_metadata,
  ge_map = adc_ge_map,
  object_name = "ADC Visium object"
)

ge_column_audit <- tibble::tibble(
  GE = names(adc_ge_map),
  source_column = unname(adc_ge_map),
  minimum = vapply(adc_metadata[, GE_RAW_COLS], min, numeric(1), na.rm = TRUE),
  median = vapply(adc_metadata[, GE_RAW_COLS], median, numeric(1), na.rm = TRUE),
  maximum = vapply(adc_metadata[, GE_RAW_COLS], max, numeric(1), na.rm = TRUE),
  fraction_zero = vapply(adc_metadata[, GE_RAW_COLS], function(x) mean(x == 0, na.rm = TRUE), numeric(1)),
  fraction_missing = vapply(adc_metadata[, GE_RAW_COLS], function(x) mean(is.na(x)), numeric(1))
)

print(ge_column_audit)

## # A tibble: 10 × 7
##    GE    source_column minimum median maximum fraction_zero fraction_missing
##    <chr> <chr>           <dbl>  <dbl>   <dbl>         <dbl>            <dbl>
##  1 GE1   raw_GE1             0 0.0629   0.282         0.349                0
##  2 GE2   raw_GE2             0 0.0914   0.343         0.333                0
##  3 GE3   raw_GE3             0 0.0819   0.382         0.336                0
##  4 GE4   raw_GE4             0 0.0416   0.186         0.492                0
##  5 GE5   raw_GE5             0 0.0607   0.330         0.379                0
##  6 GE6   raw_GE6             0 0.0710   0.314         0.350                0
##  7 GE7   raw_GE7             0 0.0805   0.303         0.340                0
##  8 GE8   raw_GE8             0 0.0912   0.345         0.333                0
##  9 GE9   raw_GE9             0 0.0952   0.289         0.320                0
## 10 GE10  raw_GE10            0 0.0821   0.293         0.336                0

readr::write_csv(
  ge_column_audit,
  file.path(CFG$output_dir, "ADC_GE_column_audit.csv")
)

The original labeling was inconsistent and did not have the correct sample ID (ADC#) or patient ID. This section will address the cross-walk and ensure the labeling is accurate.

clinical_csv = "work/InternalMedicine/s239947/ADC/adc_clinical_metadata.csv"

adc_metadata$sample_id <- as.character(adc_metadata[[CFG$sample_col]])

# Read the biopsy/patient crosswalk
clinical_crosswalk <- readr::read_csv(
  CFG$clinical_csv,
  show_col_types = FALSE
) |>
  dplyr::transmute(
    source_sample_id = trimws(as.character(sample_id)),
    adc_sample_id    = trimws(as.character(xenium_id)),
    patient_id       = trimws(as.character(patient_id))
  )

# Confirm that each identifier maps uniquely
if (anyDuplicated(clinical_crosswalk$source_sample_id)) {
  stop("The clinical file contains duplicated source sample IDs.")
}

if (anyDuplicated(clinical_crosswalk$adc_sample_id)) {
  stop("The clinical file contains duplicated ADC sample IDs.")
}

# Preserve the identifier originally stored in the Seurat object
adc_metadata$sample_id_original <-
  trimws(as.character(adc_metadata[[CFG$sample_col]]))

# First try matching existing ADC### identifiers
crosswalk_index <- match(
  adc_metadata$sample_id_original,
  clinical_crosswalk$adc_sample_id
)

# For remaining samples, try matching the SU/legacy identifiers
needs_source_match <- is.na(crosswalk_index)

crosswalk_index[needs_source_match] <- match(
  adc_metadata$sample_id_original[needs_source_match],
  clinical_crosswalk$source_sample_id
)

# Stop rather than silently retaining an unmapped sample
unmapped_ids <- sort(unique(
  adc_metadata$sample_id_original[is.na(crosswalk_index)]
))

if (length(unmapped_ids) > 0) {
  stop(
    "These Visium sample IDs were not found in the clinical crosswalk: ",
    paste(unmapped_ids, collapse = ", ")
  )
}

# Add standardized identifiers
adc_metadata$sample_id <-
  clinical_crosswalk$adc_sample_id[crosswalk_index]

adc_metadata$patient_id <-
  clinical_crosswalk$patient_id[crosswalk_index]

adc_metadata$source_sample_id <-
  clinical_crosswalk$source_sample_id[crosswalk_index]

sample_id_crosswalk_audit <- adc_metadata |>
  dplyr::distinct(
    sample_id_original,
    source_sample_id,
    sample_id,
    patient_id
  ) |>
  dplyr::arrange(sample_id)

print(sample_id_crosswalk_audit)

##    sample_id_original source_sample_id sample_id patient_id
## 1      SU22_18324_ADC   SU22_18324_ADC    ADC001     PAT_01
## 2      SU22_27508_ADC   SU22_27508_ADC    ADC002     PAT_02
## 3      SU22_29416_ADC   SU22_29416_ADC    ADC003     PAT_03
## 4      SU23_33192_ADC   SU23_33192_ADC    ADC004     PAT_04
## 5      SU21_18999_ADC   SU21_18999_ADC    ADC005     PAT_05
## 6      SU21_05811_ADC   SU21_05811_ADC    ADC007     PAT_07
## 7              ADC008   SU19_18387_ADC    ADC008     PAT_05
## 8              ADC009   SU21_36102_ADC    ADC009     PAT_08
## 9              ADC010   SU23_16186_ADC    ADC010     PAT_08
## 10             ADC011   SU24_19593_ADC    ADC011     PAT_09
## 11             ADC014   SU25_20645_ADC    ADC014     PAT_09
## 12             ADC018   SU23_09922_ADC    ADC018     PAT_10
## 13          44484_ADC        44484_ADC    ADC019     PAT_11
## 14          03119_ADC        03119_ADC    ADC020     PAT_12
## 15          37834_ADC        37834_ADC    ADC021     PAT_13

4.1 Filter to cancer-epithelial-enriched spots

tumor_filter <- derive_tumor_flag(adc_metadata, CFG)
adc_metadata$is_tumor_spot <- tumor_filter$flag
adc_metadata$passes_primary_tumor_filter <- tumor_filter$primary_flag
adc_metadata$has_any_nonzero_ge <- tumor_filter$score_flag

tumor_filter_audit <- tibble::tibble(
  metric = c(
    "total_spots",
    "primary_tumor_filter_positive",
    "any_nonzero_GE",
    "final_retained_tumor_spots",
    "primary_positive_but_all_GE_zero",
    "primary_negative_but_any_GE_nonzero"
  ),
  value = c(
    nrow(adc_metadata),
    sum(adc_metadata$passes_primary_tumor_filter),
    sum(adc_metadata$has_any_nonzero_ge),
    sum(adc_metadata$is_tumor_spot),
    sum(adc_metadata$passes_primary_tumor_filter & !adc_metadata$has_any_nonzero_ge),
    sum(!adc_metadata$passes_primary_tumor_filter & adc_metadata$has_any_nonzero_ge)
  ),
  filter_criterion = tumor_filter$criterion
)

tumor_filter_audit

readr::write_csv(
  tumor_filter_audit,
  file.path(CFG$output_dir, "ADC_tumor_spot_filter_audit.csv")
)

adc_tumor <- adc_metadata |>
  dplyr::filter(is_tumor_spot)

spots_per_sample <- adc_tumor |>
  dplyr::count(sample_id, name = "n_tumor_spots") |>
  dplyr::arrange(n_tumor_spots)

spots_per_sample

readr::write_csv(
  spots_per_sample,
  file.path(CFG$output_dir, "ADC_tumor_spots_per_sample.csv")
)

if (nrow(adc_tumor) == 0) {
  stop("The tumor-spot filter retained zero spots. Review the configured filter.")
}

if (any(!stats::complete.cases(adc_tumor[, GE_RAW_COLS, drop = FALSE]))) {
  stop(
    "At least one retained tumor spot has a missing GE score. Resolve missing ",
    "scores before scaling and clustering."
  )
}

5. Build or load sample-level clinical metadata

The preferred clinical metadata file has one row per biopsy and these columns:

• sample_id
• patient_id
• timepoint (Pre or Post)
• pair_id for matched biopsies
• adc_drug
• metastatic_site
• response_group (Responder or Nonresponder, with a documented definition)
• best_response, PFS_days, and time_to_next_treatment_days, when available
• biopsy date, ADC start/stop dates, number of cycles, and time from last dose
• receptor status (ER, PR, HER2) at the biopsied lesion
• ADC target measurement such as TROP2 or HER2 IHC/H-score, when available
• tissue/slide batch and relevant QC variables

If CFG$clinical_csv is absent, this section creates a template. Analyses that need missing fields will be skipped rather than guessing them.

sample_metadata <- tibble::tibble(
  sample_id = sort(unique(adc_tumor$sample_id))
)

object_field_map <- c(
  patient_id = CFG$patient_col,
  timepoint = CFG$timepoint_col,
  adc_drug = CFG$adc_drug_col,
  metastatic_site = CFG$metastatic_site_col,
  response_group = CFG$response_group_col,
  batch = CFG$batch_col
)

for (target_field in names(object_field_map)) {
  source_field <- object_field_map[[target_field]]

  if (is_provided(source_field) && source_field %in% colnames(adc_metadata)) {
    field_values <- adc_metadata |>
      dplyr::transmute(
        sample_id,
        value = as.character(.data[[source_field]])
      ) |>
      dplyr::filter(!is.na(value), nzchar(value)) |>
      dplyr::distinct()

    duplicated_samples <- field_values |>
      dplyr::count(sample_id) |>
      dplyr::filter(n > 1)

    if (nrow(duplicated_samples) > 0) {
      stop(
        "Metadata field ", source_field,
        " has multiple values within at least one sample. Resolve this before analysis."
      )
    }

    colnames(field_values)[2] <- target_field
    sample_metadata <- dplyr::left_join(sample_metadata, field_values, by = "sample_id")
  }
}

if (is_provided(CFG$clinical_csv)) {
  if (!file.exists(CFG$clinical_csv)) {
    stop("Configured clinical metadata file does not exist: ", CFG$clinical_csv)
  }

clinical_metadata <- readr::read_csv(
  CFG$clinical_csv,
  show_col_types = FALSE
) |>
  dplyr::mutate(
    source_sample_id = trimws(as.character(sample_id)),
    adc_sample_id    = trimws(as.character(xenium_id)),
    sample_id        = adc_sample_id,
    patient_id       = trimws(as.character(patient_id))
  ) |>
  dplyr::select(-dplyr::any_of("xenium_id"))
}

known_pairs <- tibble::tribble(
  ~sample_id, ~known_pair_id,
  "ADC005",  "PAIR_01",
  "ADC008",  "PAIR_01",
  "ADC009",  "PAIR_02",
  "ADC010",  "PAIR_02",
  "ADC017",  "PAIR_03",
  "ADC018",  "PAIR_03"
)

sample_metadata <- sample_metadata |>
  dplyr::left_join(known_pairs, by = "sample_id")

if (!"pair_id" %in% colnames(sample_metadata)) {
  sample_metadata$pair_id <- sample_metadata$known_pair_id
} else {
  sample_metadata$pair_id <- dplyr::coalesce(
    as.character(sample_metadata$pair_id),
    sample_metadata$known_pair_id
  )
}

sample_metadata <- sample_metadata |>
  dplyr::select(-known_pair_id)

template_fields <- c(
  "patient_id",
  "timepoint",
  "pair_id",
  "adc_drug",
  "metastatic_site",
  "response_group",
  "best_response",
  "PFS_days",
  "time_to_next_treatment_days",
  "biopsy_date",
  "adc_start_date",
  "adc_stop_date",
  "days_from_last_adc_dose",
  "adc_cycles_before_biopsy",
  "ER_status",
  "PR_status",
  "HER2_status",
  "ADC_target_name",
  "ADC_target_IHC_Hscore",
  "batch",
  "same_lesion_as_pair",
  "notes"
)

if (is_provided(CFG$clinical_csv)) {

  if (!file.exists(CFG$clinical_csv)) {
    stop("Clinical metadata file does not exist: ", CFG$clinical_csv)
  }

  clinical_metadata_raw <- readr::read_csv(
    CFG$clinical_csv,
    show_col_types = FALSE
  )

  # Reformat the available clinical metadata
  clinical_metadata <- clinical_metadata_raw |>
    dplyr::transmute(
      # Preserve the original SU/legacy identifier
      source_sample_id = trimws(as.character(sample_id)),

      # Use ADC### as the canonical biopsy/sample identifier
      sample_id = trimws(as.character(xenium_id)),

      # Patient-level identifier
      patient_id = trimws(as.character(patient_id)),

      sample_name = trimws(as.character(`Sample Name`)),
      adc_drug = trimws(as.character(adc_drug)),
      timepoint = trimws(as.character(timepoint)),
      response_group = trimws(as.character(response_group)),
      notes = trimws(as.character(`Other Notes`))
    )

  if (anyDuplicated(clinical_metadata$sample_id)) {
    stop("Clinical metadata contains duplicated ADC sample IDs.")
  }

  # Check whether every Visium sample has clinical metadata
  unmatched_samples <- sample_metadata |>
    dplyr::select(sample_id) |>
    dplyr::anti_join(
      clinical_metadata |> dplyr::select(sample_id),
      by = "sample_id"
    )

  if (nrow(unmatched_samples) > 0) {
    stop(
      "Clinical metadata are missing for: ",
      paste(unmatched_samples$sample_id, collapse = ", ")
    )
  }

  # Join the clinical information
  external_fields <- setdiff(
    colnames(clinical_metadata),
    "sample_id"
  )

  sample_metadata <- sample_metadata |>
    dplyr::select(-dplyr::any_of(external_fields)) |>
    dplyr::left_join(
      clinical_metadata,
      by = "sample_id"
    )
}

for (field in template_fields) {
  if (!field %in% colnames(sample_metadata)) {
    sample_metadata[[field]] <- NA_character_
  }
}

sample_metadata$timepoint <- factor(
  normalize_timepoint(sample_metadata$timepoint),
  levels = c("Pre", "Post")
)

sample_metadata <- sample_metadata |>
  dplyr::left_join(spots_per_sample, by = "sample_id") |>
  dplyr::arrange(timepoint, sample_id)

sample_metadata

readr::write_csv(
  sample_metadata,
  file.path(CFG$output_dir, "ADC_clinical_metadata_template_or_audit.csv")
)

6. ADC-internal GE standardization

For the pre/post analysis, each GE is standardized using only the ADC tumor spots. Each biopsy receives equal total weight when estimating the mean and standard deviation, so a slide with many spots does not define the scale more strongly than a slide with fewer spots.

adc_raw_matrix <- as.matrix(adc_tumor[, GE_RAW_COLS, drop = FALSE])
colnames(adc_raw_matrix) <- GE_NAMES

adc_moments <- sample_balanced_moments(
  score_matrix = adc_raw_matrix,
  sample_id = adc_tumor$sample_id
)

adc_tumor <- apply_z_scaling(
  metadata = adc_tumor,
  raw_cols = GE_RAW_COLS,
  moments = adc_moments,
  suffix = "_z_adc"
)

adc_tumor <- add_dominant_ge(
  metadata = adc_tumor,
  z_cols = GE_Z_ADC_COLS,
  label_col = "GE_label_adc",
  margin_col = "GE_label_adc_margin"
)

adc_scaling_parameters <- tibble::tibble(
  GE = GE_NAMES,
  sample_balanced_mean = unname(adc_moments$mean),
  sample_balanced_sd = unname(adc_moments$sd)
)

adc_scaling_parameters

readr::write_csv(
  adc_scaling_parameters,
  file.path(CFG$output_dir, "ADC_GE_scaling_parameters.csv")
)

7. Optional Xu reference anchoring

This section compares the datasets only in their common GE1-GE10 score space. It does not jointly integrate the full Xu scRNA-seq matrix with the Visium expression matrix.

Reference-anchored scores are useful for mapping and novelty screening. The ADC-internal scores from Section 6 remain primary for pre/post comparisons, because an apparent shift between raw Xu and Visium UCell scores could reflect platform and spot-mixture differences.

xu_available <- is_provided(CFG$xu_rds) && file.exists(CFG$xu_rds)

if (xu_available) {
  xu <- readRDS(CFG$xu_rds)
  xu_metadata <- extract_seurat_metadata(xu)

  # Identify the Xu UCell columns explicitly
  xu_ge_map <- stats::setNames(
    vapply(seq_along(GE_NAMES), function(i) {

    hits <- grep(
      pattern = paste0("^GE", i, "_UCell$"),
      x = colnames(xu_metadata),
      value = TRUE,
      ignore.case = TRUE
    )

    if (length(hits) != 1) {
      stop(
        "Expected exactly one Xu column for GE", i,
        "_UCell, but found: ",
        ifelse(length(hits) == 0, "none", paste(hits, collapse = ", "))
      )
    }

    hits
  }, character(1)),
  GE_NAMES
)

print(xu_ge_map)

# Explicitly create GE1_raw through GE10_raw
for (i in seq_along(GE_NAMES)) {

  source_column <- unname(xu_ge_map[i])
  target_column <- GE_RAW_COLS[i]

  source_values <- xu_metadata[[source_column]]

  if (is.factor(source_values)) {
    source_values <- as.character(source_values)
  }

  xu_metadata[[target_column]] <-
    suppressWarnings(as.numeric(source_values))
}

# Confirm that all canonical columns were created
missing_xu_raw <- setdiff(
  GE_RAW_COLS,
  colnames(xu_metadata)
)

if (length(missing_xu_raw) > 0) {
  stop(
    "The following Xu canonical columns were not created: ",
    paste(missing_xu_raw, collapse = ", ")
  )
}

print(
  xu_metadata |>
    dplyr::select(dplyr::all_of(GE_RAW_COLS)) |>
    summary()
)

  if (CFG$xu_object_is_cancer_epithelial) {
    xu_tumor <- xu_metadata
  } else {
    if (!is_provided(CFG$xu_epithelial_col) ||
        !CFG$xu_epithelial_col %in% colnames(xu_metadata)) {
      stop(
        "Provide a valid Xu cancer-epithelial annotation column, or use an ",
        "already-subset Xu cancer epithelial object."
      )
    }

    xu_tumor <- xu_metadata |>
      dplyr::filter(
        as.character(.data[[CFG$xu_epithelial_col]]) %in% CFG$xu_cancer_values
      )
  }

  xu_raw_matrix <- as.matrix(xu_tumor[, GE_RAW_COLS, drop = FALSE])
  colnames(xu_raw_matrix) <- GE_NAMES
  xu_moments <- ordinary_moments(xu_raw_matrix)

  xu_tumor <- apply_z_scaling(
    metadata = xu_tumor,
    raw_cols = GE_RAW_COLS,
    moments = xu_moments,
    suffix = "_z_ref"
  ) |>
    add_dominant_ge(
      z_cols = GE_Z_REF_COLS,
      label_col = "GE_label_ref",
      margin_col = "GE_label_ref_margin"
    )

  adc_tumor <- apply_z_scaling(
    metadata = adc_tumor,
    raw_cols = GE_RAW_COLS,
    moments = xu_moments,
    suffix = "_z_ref"
  ) |>
    add_dominant_ge(
      z_cols = GE_Z_REF_COLS,
      label_col = "GE_label_ref",
      margin_col = "GE_label_ref_margin"
    )

  xu_reference_parameters <- tibble::tibble(
    GE = GE_NAMES,
    xu_mean = unname(xu_moments$mean),
    xu_sd = unname(xu_moments$sd)
  )

  readr::write_csv(
    xu_reference_parameters,
    file.path(CFG$output_dir, "Xu_GE_reference_parameters.csv")
  )

  # Xu GE centroids and ADC nearest-centroid assignment -----------------------
  xu_centroid_table <- xu_tumor |>
    dplyr::filter(!is.na(GE_label_ref)) |>
    dplyr::group_by(GE_label_ref) |>
    dplyr::summarise(
      dplyr::across(dplyr::all_of(GE_Z_REF_COLS), mean, na.rm = TRUE),
      .groups = "drop"
    )

  xu_centroids <- as.matrix(xu_centroid_table[, GE_Z_REF_COLS, drop = FALSE])
  rownames(xu_centroids) <- as.character(xu_centroid_table$GE_label_ref)

  xu_nearest <- nearest_centroid(
    as.matrix(xu_tumor[, GE_Z_REF_COLS, drop = FALSE]),
    xu_centroids
  )

  adc_nearest <- nearest_centroid(
    as.matrix(adc_tumor[, GE_Z_REF_COLS, drop = FALSE]),
    xu_centroids
  )

  xu_tumor$GE_nearest_xu_centroid <- xu_nearest$nearest_label
  xu_tumor$GE_xu_centroid_distance <- xu_nearest$distance
  adc_tumor$GE_nearest_xu_centroid <- adc_nearest$nearest_label
  adc_tumor$GE_xu_centroid_distance <- adc_nearest$distance

  distance_thresholds <- xu_tumor |>
    dplyr::group_by(GE_nearest_xu_centroid) |>
    dplyr::summarise(
      distance_threshold = stats::quantile(
        GE_xu_centroid_distance,
        probs = CFG$out_of_reference_quantile,
        na.rm = TRUE
      ),
      .groups = "drop"
    )

  adc_tumor <- adc_tumor |>
    dplyr::left_join(distance_thresholds, by = "GE_nearest_xu_centroid") |>
    dplyr::mutate(
      GE_out_of_xu_reference = GE_xu_centroid_distance > distance_threshold
    )

  readr::write_csv(
    xu_centroid_table,
    file.path(CFG$output_dir, "Xu_GE_centroids.csv")
  )

  # Balanced joint PCA for visualization only --------------------------------
  xu_display <- balanced_sample(
    xu_tumor,
    group_col = "GE_label_ref",
    max_per_group = CFG$max_xu_cells_per_ge_for_display
  )

  adc_display_ref <- balanced_sample(
    adc_tumor,
    group_col = "sample_id",
    max_per_group = CFG$max_spots_per_sample_for_clustering
  )

  joint_matrix <- rbind(
    as.matrix(xu_display[, GE_Z_REF_COLS, drop = FALSE]),
    as.matrix(adc_display_ref[, GE_Z_REF_COLS, drop = FALSE])
  )

  joint_pca <- stats::prcomp(joint_matrix, center = FALSE, scale. = FALSE)

  joint_coordinates <- tibble::tibble(
    PC1 = joint_pca$x[, 1],
    PC2 = joint_pca$x[, 2],
    source = c(
      rep("Xu cancer epithelial cells", nrow(xu_display)),
      rep("ADC Visium tumor spots", nrow(adc_display_ref))
    ),
    GE = c(
      as.character(xu_display$GE_label_ref),
      as.character(adc_display_ref$GE_label_ref)
    )
  )

  p_joint_pca <- ggplot(joint_coordinates, aes(PC1, PC2, color = GE)) +
    geom_point(alpha = 0.35, size = 0.6) +
    facet_wrap(~source) +
    scale_color_brewer(palette = "Paired", drop = FALSE) +
    theme_classic() +
    labs(
      title = "Xu and ADC observations in shared GE-score space",
      subtitle = "Visualization only; not a whole-transcriptome integration"
    )

  print(p_joint_pca)

  ggsave(
    file.path(CFG$output_dir, "Xu_ADC_shared_GE_PCA.pdf"),
    p_joint_pca,
    width = 10,
    height = 5
  )
} else {
  message(
    "Xu reference object was not provided. ADC-only sections will run; ",
    "reference anchoring and out-of-reference screening will be skipped."
  )
}

##          GE1          GE2          GE3          GE4          GE5          GE6 
##  "GE1_UCell"  "GE2_UCell"  "GE3_UCell"  "GE4_UCell"  "GE5_UCell"  "GE6_UCell" 
##          GE7          GE8          GE9         GE10 
##  "GE7_UCell"  "GE8_UCell"  "GE9_UCell" "GE10_UCell" 
##     raw_GE1           raw_GE2            raw_GE3           raw_GE4       
##  Min.   :0.00000   Min.   :0.008924   Min.   :0.00000   Min.   :0.00000  
##  1st Qu.:0.07447   1st Qu.:0.105247   1st Qu.:0.07023   1st Qu.:0.03724  
##  Median :0.09988   Median :0.129784   Median :0.09063   Median :0.05252  
##  Mean   :0.12023   Mean   :0.135443   Mean   :0.10101   Mean   :0.05609  
##  3rd Qu.:0.13133   3rd Qu.:0.159093   3rd Qu.:0.10983   3rd Qu.:0.06700  
##  Max.   :0.50611   Max.   :0.353101   Max.   :0.45526   Max.   :0.52068  
##     raw_GE5          raw_GE6          raw_GE7          raw_GE8        
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.009155  
##  1st Qu.:0.1013   1st Qu.:0.1603   1st Qu.:0.1019   1st Qu.:0.125212  
##  Median :0.2716   Median :0.2062   Median :0.1388   Median :0.154753  
##  Mean   :0.2380   Mean   :0.2223   Mean   :0.1482   Mean   :0.167139  
##  3rd Qu.:0.3537   3rd Qu.:0.2563   3rd Qu.:0.1793   3rd Qu.:0.193093  
##  Max.   :0.5511   Max.   :0.7196   Max.   :0.5878   Max.   :0.535678  
##     raw_GE9           raw_GE10     
##  Min.   :0.00000   Min.   :0.0000  
##  1st Qu.:0.07649   1st Qu.:0.1236  
##  Median :0.09805   Median :0.1536  
##  Mean   :0.10696   Mean   :0.1588  
##  3rd Qu.:0.12898   3rd Qu.:0.1894  
##  Max.   :0.39365   Max.   :0.5050

8. Secondary unsupervised clustering of ADC GE profiles

Dominant-GE labels reproduce the published Xu approach. This secondary analysis asks whether ADC tumor spots form broader or mixed GE states. Clustering is fit on a sample-balanced subset so a large slide does not dominate. All retained spots are then assigned to the nearest fitted cluster centroid.

adc_cluster_fit <- balanced_sample(
  adc_tumor,
  group_col = "sample_id",
  max_per_group = CFG$max_spots_per_sample_for_clustering
)

cluster_matrix <- as.matrix(adc_cluster_fit[, GE_Z_ADC_COLS, drop = FALSE])

silhouette_data <- safe_sample_rows(
  adc_cluster_fit,
  n = CFG$silhouette_max_n
)
silhouette_matrix <- as.matrix(silhouette_data[, GE_Z_ADC_COLS, drop = FALSE])

valid_k <- CFG$k_range[
  CFG$k_range >= 2 & CFG$k_range < nrow(silhouette_matrix)
]

if (length(valid_k) == 0) {
  stop("Too few tumor spots for the requested clustering range.")
}

silhouette_results <- purrr::map_dfr(valid_k, function(k) {
  fit <- stats::kmeans(silhouette_matrix, centers = k, nstart = 50)
  sil <- cluster::silhouette(fit$cluster, stats::dist(silhouette_matrix))

  tibble::tibble(
    k = k,
    mean_silhouette = mean(sil[, "sil_width"])
  )
})

selected_k <- if (is.null(CFG$k_ge_clusters)) {
  silhouette_results$k[which.max(silhouette_results$mean_silhouette)]
} else {
  CFG$k_ge_clusters
}

if (length(selected_k) != 1 || !is.finite(selected_k) ||
    selected_k < 2 || selected_k >= nrow(cluster_matrix)) {
  stop("CFG$k_ge_clusters must be one valid integer between 2 and n_spots - 1.")
}

silhouette_results

selected_k

## [1] 2

readr::write_csv(
  silhouette_results,
  file.path(CFG$output_dir, "ADC_GE_cluster_silhouette_screen.csv")
)

final_kmeans <- stats::kmeans(
  cluster_matrix,
  centers = selected_k,
  nstart = 100
)

all_assignments <- nearest_centroid(
  as.matrix(adc_tumor[, GE_Z_ADC_COLS, drop = FALSE]),
  final_kmeans$centers
)

dominant_center_ge <- GE_NAMES[
  max.col(final_kmeans$centers, ties.method = "first")
]

cluster_names <- paste0(
  "C",
  seq_len(selected_k),
  "_",
  dominant_center_ge
)

adc_tumor$GE_profile_cluster <- factor(
  cluster_names[as.integer(all_assignments$nearest_label)],
  levels = cluster_names
)
adc_tumor$GE_profile_cluster_distance <- all_assignments$distance

cluster_centers <- as.data.frame(final_kmeans$centers)
colnames(cluster_centers) <- GE_NAMES
cluster_centers$cluster <- cluster_names

readr::write_csv(
  cluster_centers,
  file.path(CFG$output_dir, "ADC_GE_profile_cluster_centers.csv")
)

p_silhouette <- ggplot(silhouette_results, aes(k, mean_silhouette)) +
  geom_line() +
  geom_point(size = 2) +
  geom_vline(xintercept = selected_k, linetype = 2, color = "firebrick") +
  scale_x_continuous(breaks = valid_k) +
  theme_classic() +
  labs(
    title = "Selection of the number of GE-profile clusters",
    y = "Mean silhouette width"
  )

print(p_silhouette)

ggsave(
  file.path(CFG$output_dir, "ADC_GE_cluster_silhouette_screen.pdf"),
  p_silhouette,
  width = 6,
  height = 4
)

cluster_center_matrix <- as.matrix(
  cluster_centers[, GE_NAMES, drop = FALSE]
)

rownames(cluster_center_matrix) <- cluster_centers$cluster

pheatmap::pheatmap(
  cluster_center_matrix,
  cluster_rows = TRUE,
  cluster_cols = FALSE,
  border_color = NA,
  main = "Mean ADC-internal GE z-score by cluster",
  filename = file.path(CFG$output_dir, "ADC_GE_profile_cluster_heatmap.pdf"),
  width = 8,
  height = max(4, selected_k * 0.45)
)

pheatmap::pheatmap(
  cluster_center_matrix,
  cluster_rows = TRUE,
  cluster_cols = FALSE,
  border_color = NA,
  main = "Mean ADC-internal GE z-score by cluster"
)

8.1 PCA of ADC tumor spots

adc_pca <- stats::prcomp(cluster_matrix, center = FALSE, scale. = FALSE)

adc_pca_coordinates <- tibble::tibble(
  PC1 = adc_pca$x[, 1],
  PC2 = adc_pca$x[, 2],
  sample_id = adc_cluster_fit$sample_id,
  dominant_GE = adc_cluster_fit$GE_label_adc,
  GE_profile_cluster = factor(
    cluster_names[final_kmeans$cluster],
    levels = cluster_names
  )
)

p_adc_pca_ge <- ggplot(adc_pca_coordinates, aes(PC1, PC2, color = dominant_GE)) +
  geom_point(alpha = 0.45, size = 0.7) +
  scale_color_brewer(palette = "Paired", drop = FALSE) +
  theme_classic() +
  labs(title = "ADC tumor spots colored by dominant GE")

p_adc_pca_cluster <- ggplot(
  adc_pca_coordinates,
  aes(PC1, PC2, color = GE_profile_cluster)
) +
  geom_point(alpha = 0.45, size = 0.7) +
  theme_classic() +
  labs(title = "ADC tumor spots colored by GE-profile cluster")

p_adc_pca_ge + p_adc_pca_cluster

ggsave(
  file.path(CFG$output_dir, "ADC_GE_PCA.pdf"),
  p_adc_pca_ge + p_adc_pca_cluster,
  width = 12,
  height = 5
)

9. Return GE states to the Seurat object

Non-tumor spots receive NA for derived GE labels and clusters.

derived_columns <- c(
  GE_Z_ADC_COLS,
  "GE_label_adc",
  "GE_label_adc_margin",
  "GE_label_adc_max_z",
  "GE_profile_cluster",
  "GE_profile_cluster_distance"
)

if (xu_available) {
  derived_columns <- c(
    derived_columns,
    GE_Z_REF_COLS,
    "GE_label_ref",
    "GE_label_ref_margin",
    "GE_label_ref_max_z",
    "GE_nearest_xu_centroid",
    "GE_xu_centroid_distance",
    "GE_out_of_xu_reference"
  )
}

derived_columns <- intersect(derived_columns, colnames(adc_tumor))

adc_metadata_augmented <- adc_metadata

for (column in derived_columns) {
  values <- adc_tumor[[column]]
  names(values) <- adc_tumor$spot_id
  adc_metadata_augmented[[column]] <- values[adc_metadata_augmented$spot_id]
}

adc <- add_metadata_to_seurat(
  object = adc,
  metadata = adc_metadata_augmented,
  columns = c("is_tumor_spot", derived_columns)
)

if (isTRUE(CFG$save_augmented_rds)) {
  saveRDS(
    adc,
    file.path(CFG$output_dir, "ADC_Visium_with_GE_states.rds")
  )
}

10. Sample-level GE summaries

Two complementary summaries are generated:

1. median standardized GE score per biopsy;
2. proportion of tumor spots assigned to each dominant GE.

The score analysis captures continuous shifts. The composition analysis captures changes in the prevalence of GE states.

dominant_label_col <- if (xu_available) "GE_label_ref" else "GE_label_adc"

sample_ge_scores <- adc_tumor |>
  dplyr::select(sample_id, dplyr::all_of(GE_Z_ADC_COLS)) |>
  tidyr::pivot_longer(
    cols = dplyr::all_of(GE_Z_ADC_COLS),
    names_to = "GE",
    values_to = "z_score"
  ) |>
  dplyr::mutate(GE = stringr::str_remove(GE, "_z_adc$")) |>
  dplyr::group_by(sample_id, GE) |>
  dplyr::summarise(
    median_score = median(z_score, na.rm = TRUE),
    mean_score = mean(z_score, na.rm = TRUE),
    q25 = stats::quantile(z_score, 0.25, na.rm = TRUE),
    q75 = stats::quantile(z_score, 0.75, na.rm = TRUE),
    n_tumor_spots = dplyr::n(),
    .groups = "drop"
  )

sample_ge_proportions <- adc_tumor |>
  dplyr::transmute(
    sample_id,
    GE = as.character(.data[[dominant_label_col]])
  ) |>
  dplyr::filter(!is.na(GE)) |>
  dplyr::count(sample_id, GE, name = "n") |>
  tidyr::complete(sample_id, GE = GE_NAMES, fill = list(n = 0)) |>
  dplyr::group_by(sample_id) |>
  dplyr::mutate(proportion = n / sum(n)) |>
  dplyr::ungroup()

sample_ge_entropy <- sample_ge_proportions |>
  dplyr::group_by(sample_id) |>
  dplyr::summarise(
    GE_entropy = {
      p <- proportion[proportion > 0]
      -sum(p * log(p)) / log(length(GE_NAMES))
    },
    dominant_GE_fraction = max(proportion),
    .groups = "drop"
  )

sample_metadata_analysis <- sample_metadata |>
  dplyr::left_join(sample_ge_entropy, by = "sample_id")

readr::write_csv(
  sample_ge_scores,
  file.path(CFG$output_dir, "ADC_sample_GE_scores.csv")
)

readr::write_csv(
  sample_ge_proportions,
  file.path(CFG$output_dir, "ADC_sample_GE_proportions.csv")
)

readr::write_csv(
  sample_metadata_analysis,
  file.path(CFG$output_dir, "ADC_sample_GE_heterogeneity.csv")
)

sample_metadata_analysis |>
  dplyr::select(
    sample_id,
    timepoint,
    adc_drug,
    metastatic_site,
    n_tumor_spots,
    GE_entropy,
    dominant_GE_fraction
  )

10.1 Sample-by-GE heatmap

sample_score_matrix <- sample_ge_scores |>
  dplyr::select(sample_id, GE, median_score) |>
  tidyr::pivot_wider(names_from = GE, values_from = median_score) |>
  dplyr::select(sample_id, dplyr::all_of(GE_NAMES)) |>
  tibble::column_to_rownames("sample_id") |>
  as.matrix()

annotation_fields <- c("timepoint", "adc_drug", "metastatic_site", "response_group")

sample_annotation <- sample_metadata_analysis |>
  dplyr::filter(sample_id %in% rownames(sample_score_matrix)) |>
  dplyr::select(sample_id, dplyr::any_of(annotation_fields)) |>
  tibble::column_to_rownames("sample_id")

sample_annotation <- sample_annotation[
  rownames(sample_score_matrix),
  ,
  drop = FALSE
]

keep_annotation <- vapply(
  sample_annotation,
  function(x) any(!is.na(x) & nzchar(as.character(x))),
  logical(1)
)

sample_annotation <- sample_annotation[, keep_annotation, drop = FALSE]

pheatmap::pheatmap(
  sample_score_matrix,
  cluster_rows = TRUE,
  cluster_cols = FALSE,
  annotation_row = if (ncol(sample_annotation) > 0) sample_annotation else NULL,
  border_color = NA,
  main = "Median ADC-internal GE score per biopsy",
  filename = file.path(CFG$output_dir, "ADC_sample_GE_score_heatmap.pdf"),
  width = 9,
  height = max(5, nrow(sample_score_matrix) * 0.30)
)

pheatmap::pheatmap(
  sample_score_matrix,
  cluster_rows = TRUE,
  cluster_cols = FALSE,
  annotation_row = if (ncol(sample_annotation) > 0) sample_annotation else NULL,
  border_color = NA,
  main = "Median ADC-internal GE score per biopsy"
)

10.2 GE-state composition by sample

sample_order <- sample_metadata_analysis |>
  dplyr::arrange(timepoint, adc_drug, metastatic_site, sample_id) |>
  dplyr::pull(sample_id)

p_ge_composition <- sample_ge_proportions |>
  dplyr::left_join(
    sample_metadata_analysis |>
      dplyr::select(sample_id, timepoint, adc_drug, metastatic_site),
    by = "sample_id"
  ) |>
  dplyr::mutate(
    sample_id = factor(sample_id, levels = sample_order),
    GE = factor(GE, levels = GE_NAMES)
  ) |>
  ggplot(aes(sample_id, proportion, fill = GE)) +
  geom_col(width = 0.85) +
  scale_fill_brewer(palette = "Paired", drop = FALSE) +
  scale_y_continuous(labels = scales::percent) +
  facet_grid(
    rows = vars(timepoint),
    scales = "free_y",
    space = "free_y"
  ) +
  coord_flip() +
  theme_classic() +
  labs(
    x = NULL,
    y = "Fraction of tumor spots",
    title = "GE-state composition of each biopsy"
  )

print(p_ge_composition)

ggsave(
  file.path(CFG$output_dir, "ADC_sample_GE_composition.pdf"),
  p_ge_composition,
  width = 9,
  height = max(6, length(sample_order) * 0.35)
)

10.3 Xu out-of-reference fraction by sample

An out-of-reference spot is farther from its nearest Xu GE centroid than the configured percentile of Xu cells assigned to that centroid. Because Visium spots are mixtures and Xu observations are single cells, this is a candidate novel/residual-state flag, not proof of a new state.

if (xu_available) {
  sample_out_of_reference <- adc_tumor |>
    dplyr::group_by(sample_id) |>
    dplyr::summarise(
      out_of_reference_fraction = mean(GE_out_of_xu_reference, na.rm = TRUE),
      median_xu_centroid_distance = median(GE_xu_centroid_distance, na.rm = TRUE),
      .groups = "drop"
    ) |>
    dplyr::left_join(sample_metadata_analysis, by = "sample_id")

  print(sample_out_of_reference)

  readr::write_csv(
    sample_out_of_reference,
    file.path(CFG$output_dir, "ADC_sample_out_of_Xu_reference.csv")
  )
}

## # A tibble: 15 × 30
##    sample_id out_of_reference_fraction median_xu_centroid_distance pair_id
##    <chr>                         <dbl>                       <dbl> <chr>  
##  1 ADC001                       0.167                         2.46 <NA>   
##  2 ADC002                       0.264                         3.02 <NA>   
##  3 ADC003                       0.115                         3.01 <NA>   
##  4 ADC004                       0.0159                        2.24 <NA>   
##  5 ADC005                       0.369                         3.16 PAIR_01
##  6 ADC007                       0.536                         3.49 <NA>   
##  7 ADC008                       0.837                         4.24 PAIR_01
##  8 ADC009                       0.549                         3.55 PAIR_02
##  9 ADC010                       0.134                         2.38 PAIR_02
## 10 ADC011                       0.0124                        2.17 <NA>   
## 11 ADC014                       0.106                         2.31 <NA>   
## 12 ADC018                       0.582                         3.53 PAIR_03
## 13 ADC019                       0.0635                        1.97 <NA>   
## 14 ADC020                       0.0432                        1.99 <NA>   
## 15 ADC021                       0.469                         3.29 <NA>   
## # ℹ 26 more variables: source_sample_id <chr>, patient_id <chr>,
## #   sample_name <chr>, adc_drug <chr>, timepoint <fct>, response_group <chr>,
## #   notes <chr>, metastatic_site <chr>, best_response <chr>, PFS_days <chr>,
## #   time_to_next_treatment_days <chr>, biopsy_date <chr>, adc_start_date <chr>,
## #   adc_stop_date <chr>, days_from_last_adc_dose <chr>,
## #   adc_cycles_before_biopsy <chr>, ER_status <chr>, PR_status <chr>,
## #   HER2_status <chr>, ADC_target_name <chr>, ADC_target_IHC_Hscore <chr>, …

11. Exploratory pre-ADC versus post-ADC comparison

This analysis runs only after timepoint is available. For each GE it reports:

• the sample-level mean difference (post minus pre) in median GE score;
• Cliff’s delta as a rank-based effect size;
• an exact label-permutation p-value when the number of assignments is small;
• Benjamini-Hochberg FDR across the 10 GEs.

This cross-sectional comparison is vulnerable to confounding by ADC drug, metastatic site, tumor subtype, and biopsy timing. Effect sizes and consistency are more informative than a binary p-value in this 6-pre/11-post cohort.

has_timepoint <- sum(sample_metadata_analysis$timepoint == "Pre", na.rm = TRUE) >= 2 &&
  sum(sample_metadata_analysis$timepoint == "Post", na.rm = TRUE) >= 2

if (has_timepoint) {
  prepost_score_data <- sample_ge_scores |>
    dplyr::left_join(
      sample_metadata_analysis |>
        dplyr::select(
          sample_id,
          patient_id,
          pair_id,
          timepoint,
          adc_drug,
          metastatic_site,
          response_group
        ),
      by = "sample_id"
    ) |>
    dplyr::filter(timepoint %in% c("Pre", "Post"))

  prepost_ge_results <- prepost_score_data |>
    dplyr::group_by(GE) |>
    dplyr::group_modify(function(data, key) {
      permutation <- permutation_mean_difference(
        values = data$median_score,
        groups = data$timepoint
      )

      pre_values <- data$median_score[data$timepoint == "Pre"]
      post_values <- data$median_score[data$timepoint == "Post"]

      tibble::tibble(
        n_pre = length(pre_values),
        n_post = length(post_values),
        mean_pre = mean(pre_values),
        mean_post = mean(post_values),
        median_pre = median(pre_values),
        median_post = median(post_values),
        mean_difference_post_minus_pre = permutation$effect,
        cliffs_delta_post_vs_pre = cliffs_delta(post_values, pre_values),
        p_value = permutation$p_value,
        permutation_method = permutation$method,
        prevalence_pre_above_ADC_center = mean(pre_values > 0),
        prevalence_post_above_ADC_center = mean(post_values > 0)
      )
    }) |>
    dplyr::ungroup() |>
    dplyr::mutate(
      FDR = p.adjust(p_value, method = "BH"),
      stable_high_candidate =
        abs(mean_difference_post_minus_pre) <= CFG$conservation_margin_z &
        prevalence_pre_above_ADC_center >= 0.5 &
        prevalence_post_above_ADC_center >= 0.5
    ) |>
    dplyr::arrange(FDR, dplyr::desc(abs(mean_difference_post_minus_pre)))

  prepost_ge_results

  readr::write_csv(
    prepost_ge_results,
    file.path(CFG$output_dir, "ADC_pre_vs_post_GE_score_results.csv")
  )

  paired_line_data <- prepost_score_data |>
    dplyr::filter(!is.na(pair_id), nzchar(pair_id))

  p_prepost_scores <- ggplot(
    prepost_score_data,
    aes(timepoint, median_score, color = timepoint)
  ) +
    geom_boxplot(outlier.shape = NA, width = 0.55, alpha = 0.15) +
    geom_line(
      data = paired_line_data,
      aes(group = pair_id),
      color = "grey55",
      linewidth = 0.45
    ) +
    geom_point(
      position = position_jitter(width = 0.06, height = 0),
      size = 2
    ) +
    facet_wrap(~GE, ncol = 5, scales = "free_y") +
    scale_color_manual(values = c(Pre = "#2166AC", Post = "#B2182B")) +
    theme_classic() +
    theme(legend.position = "bottom") +
    labs(
      x = NULL,
      y = "Median ADC-internal GE z-score",
      title = "Biopsy-level GE scores before and after ADC treatment",
      subtitle = "Grey lines connect known matched pairs"
    )

  print(p_prepost_scores)

  ggsave(
    file.path(CFG$output_dir, "ADC_pre_vs_post_GE_scores.pdf"),
    p_prepost_scores,
    width = 13,
    height = 7
  )
} else {
  message(
    "Pre/post comparison skipped. Add standardized Pre/Post labels to the ",
    "clinical metadata and rerun."
  )
}

11.1 Compositional pre/post analysis

GE-state proportions sum to one, so they are compositional. This section adds 0.5 to each GE count and applies a centered log-ratio (CLR) transform before the same biopsy-level permutation comparison.

sample_ge_clr <- sample_ge_proportions |>
  dplyr::group_by(sample_id) |>
  dplyr::mutate(
    adjusted_proportion = (n + 0.5) / (sum(n) + 0.5 * length(GE_NAMES)),
    clr = log(adjusted_proportion) - mean(log(adjusted_proportion))
  ) |>
  dplyr::ungroup() |>
  dplyr::left_join(
    sample_metadata_analysis |>
      dplyr::select(sample_id, timepoint, adc_drug, metastatic_site, pair_id),
    by = "sample_id"
  )

if (has_timepoint) {
  prepost_composition_results <- sample_ge_clr |>
    dplyr::filter(timepoint %in% c("Pre", "Post")) |>
    dplyr::group_by(GE) |>
    dplyr::group_modify(function(data, key) {
      permutation <- permutation_mean_difference(data$clr, data$timepoint)

      tibble::tibble(
        n_pre = sum(data$timepoint == "Pre"),
        n_post = sum(data$timepoint == "Post"),
        mean_clr_difference_post_minus_pre = permutation$effect,
        p_value = permutation$p_value,
        permutation_method = permutation$method
      )
    }) |>
    dplyr::ungroup() |>
    dplyr::mutate(FDR = p.adjust(p_value, method = "BH")) |>
    dplyr::arrange(FDR)

  print(prepost_composition_results)

  readr::write_csv(
    prepost_composition_results,
    file.path(CFG$output_dir, "ADC_pre_vs_post_GE_composition_results.csv")
  )
}

## # A tibble: 10 × 7
##    GE    n_pre n_post mean_clr_difference_pos…¹ p_value permutation_method   FDR
##    <chr> <int>  <int>                     <dbl>   <dbl> <chr>              <dbl>
##  1 GE9       5     10                    -1.53   0.0196 exact label permu… 0.196
##  2 GE10      5     10                     0.597  0.0966 exact label permu… 0.317
##  3 GE2       5     10                     1.28   0.169  exact label permu… 0.317
##  4 GE4       5     10                     0.984  0.110  exact label permu… 0.317
##  5 GE5       5     10                    -0.626  0.184  exact label permu… 0.317
##  6 GE7       5     10                    -0.665  0.190  exact label permu… 0.317
##  7 GE6       5     10                     0.217  0.554  exact label permu… 0.792
##  8 GE8       5     10                    -0.286  0.666  exact label permu… 0.833
##  9 GE3       5     10                     0.174  0.759  exact label permu… 0.843
## 10 GE1       5     10                    -0.140  0.868  exact label permu… 0.868
## # ℹ abbreviated name: ¹​mean_clr_difference_post_minus_pre

12. Matched-pair change and conservation

With three pairs, a two-sided paired significance test cannot provide strong evidence. The useful outputs are the direction and magnitude of each change, agreement across pairs, and similarity of the entire 10-GE profile.

Calling a GE “conserved” requires a prespecified biological equivalence margin and more matched pairs. Here, conservation_margin_z is only a screening threshold for stable candidates.

if (has_timepoint) {
  matched_score_data <- sample_ge_scores |>
    dplyr::left_join(
      sample_metadata_analysis |>
        dplyr::select(sample_id, pair_id, timepoint, same_lesion_as_pair),
      by = "sample_id"
    ) |>
    dplyr::filter(
      !is.na(pair_id),
      nzchar(pair_id),
      timepoint %in% c("Pre", "Post")
    )

  matched_ge_deltas <- matched_score_data |>
    dplyr::select(pair_id, sample_id, timepoint, GE, median_score) |>
    tidyr::pivot_wider(
      id_cols = c(pair_id, GE),
      names_from = timepoint,
      values_from = median_score
    ) |>
    dplyr::filter(!is.na(Pre), !is.na(Post)) |>
    dplyr::mutate(
      delta_post_minus_pre = Post - Pre,
      absolute_delta = abs(delta_post_minus_pre),
      stable_within_margin = absolute_delta <= CFG$conservation_margin_z
    )

  matched_ge_summary <- matched_ge_deltas |>
    dplyr::group_by(GE) |>
    dplyr::summarise(
      n_pairs = dplyr::n(),
      median_delta = median(delta_post_minus_pre),
      median_absolute_delta = median(absolute_delta),
      fraction_increased = mean(delta_post_minus_pre > 0),
      fraction_stable_within_margin = mean(stable_within_margin),
      .groups = "drop"
    ) |>
    dplyr::arrange(median_absolute_delta)

  matched_profile_similarity <- matched_score_data |>
    dplyr::group_by(pair_id) |>
    dplyr::group_modify(function(data, key) {
      pre <- data |>
        dplyr::filter(timepoint == "Pre") |>
        dplyr::arrange(match(GE, GE_NAMES))

      post <- data |>
        dplyr::filter(timepoint == "Post") |>
        dplyr::arrange(match(GE, GE_NAMES))

      if (nrow(pre) != 10 || nrow(post) != 10) {
        return(tibble::tibble(
          spearman_GE_profile = NA_real_,
          pearson_GE_profile = NA_real_,
          top3_GE_jaccard = NA_real_
        ))
      }

      pre_values <- setNames(pre$median_score, pre$GE)
      post_values <- setNames(post$median_score, post$GE)

      tibble::tibble(
        spearman_GE_profile = stats::cor(pre_values, post_values, method = "spearman"),
        pearson_GE_profile = stats::cor(pre_values, post_values, method = "pearson"),
        top3_GE_jaccard = top_ge_jaccard(pre_values, post_values, n_top = 3)
      )
    }) |>
    dplyr::ungroup()

  matched_ge_deltas
  matched_ge_summary
  matched_profile_similarity

  readr::write_csv(
    matched_ge_deltas,
    file.path(CFG$output_dir, "ADC_matched_pair_GE_deltas.csv")
  )

  readr::write_csv(
    matched_ge_summary,
    file.path(CFG$output_dir, "ADC_matched_pair_GE_stability_summary.csv")
  )

  readr::write_csv(
    matched_profile_similarity,
    file.path(CFG$output_dir, "ADC_matched_pair_GE_profile_similarity.csv")
  )

  p_matched_delta <- ggplot(
    matched_ge_deltas,
    aes(GE, delta_post_minus_pre, group = pair_id, color = pair_id)
  ) +
    geom_hline(yintercept = 0, linetype = 2, color = "grey50") +
    geom_line(linewidth = 0.6) +
    geom_point(size = 2) +
    theme_classic() +
    labs(
      x = NULL,
      y = "Post minus pre median GE z-score",
      title = "Within-pair changes across the 10 GEs"
    )

  print(p_matched_delta)

  ggsave(
    file.path(CFG$output_dir, "ADC_matched_pair_GE_deltas.pdf"),
    p_matched_delta,
    width = 9,
    height = 5
  )
}

13. Exploratory association with clinical response

This section uses pretreatment biopsies only when asking about a predictive signal. With approximately six pretreatment samples, it is a hypothesis screen, not a trained or validated biomarker model. Do not fit a multivariable machine- learning classifier to this cohort.

The response definition and its time relationship to each biopsy must be fixed before unblinding the GE results.

has_response <- "response_group" %in% colnames(sample_metadata_analysis) &&
  sum(sample_metadata_analysis$response_group %in% c("Y", "N"), na.rm = TRUE) >= 4

if (has_response) {
  pretreatment_response_data <- sample_ge_scores |>
    dplyr::left_join(
      sample_metadata_analysis |>
        dplyr::select(sample_id, timepoint, response_group, adc_drug, metastatic_site),
      by = "sample_id"
    ) |>
    dplyr::filter(
      timepoint == "Pre",
      response_group %in% c("Y", "N")
    )

  if (dplyr::n_distinct(pretreatment_response_data$sample_id) >= 2 &&
      dplyr::n_distinct(pretreatment_response_data$response_group) == 2) {
    pretreatment_response_results <- pretreatment_response_data |>
      dplyr::group_by(GE) |>
      dplyr::summarise(
        n_responder = sum(response_group == "Y"),
        n_nonresponder = sum(response_group == "N"),
        mean_responder = mean(median_score[response_group == "Y"]),
        mean_nonresponder = mean(median_score[response_group == "N"]),
        mean_difference_responder_minus_nonresponder =
          mean_responder - mean_nonresponder,
        cliffs_delta_responder_vs_nonresponder = cliffs_delta(
          median_score[response_group == "Y"],
          median_score[response_group == "N"]
        ),
        .groups = "drop"
      ) |>
      dplyr::arrange(dplyr::desc(abs(mean_difference_responder_minus_nonresponder)))

    pretreatment_response_results

    readr::write_csv(
      pretreatment_response_results,
      file.path(CFG$output_dir, "ADC_pretreatment_GE_response_screen.csv")
    )

    p_response <- ggplot(
      pretreatment_response_data,
      aes(response_group, median_score, color = response_group)
    ) +
      geom_boxplot(outlier.shape = NA, alpha = 0.15) +
      geom_point(size = 2, position = position_jitter(width = 0.06)) +
      facet_wrap(~GE, ncol = 5, scales = "free_y") +
      theme_classic() +
      theme(legend.position = "none") +
      labs(
        x = NULL,
        y = "Pretreatment median GE z-score",
        title = "Hypothesis-generating pretreatment response screen"
      )

    print(p_response)
  } else {
    message("Too few pretreatment samples with both response groups for the response screen.")
  }
} else {
  message(
    "Response screen skipped. Add a prespecified Responder/Nonresponder field ",
    "and document the response definition."
  )
}

14. Spatial visualization

This section plots the dominant GE and unsupervised GE-profile cluster on each Visium image. Review these maps next to H&E and tumor annotations. A cluster confined to tissue edges, necrosis, low-RNA regions, or a single slide should be treated as a technical candidate until validated.

adc <- SeuratObject::UpdateSeuratObject(adc)

spatial_label <- if (xu_available) {
  "GE_label_ref"
} else {
  "GE_label_adc"
}

required_spatial_fields <- c(
  spatial_label,
  "GE_profile_cluster"
)

missing_spatial_fields <- setdiff(
  required_spatial_fields,
  colnames(adc[[]])
)

if (length(missing_spatial_fields) > 0) {
  stop(
    "The following GE metadata were not added back to adc: ",
    paste(missing_spatial_fields, collapse = ", "),
    ". Rerun section 9."
  )
}

test_image <- Seurat::Images(adc)[1]

Seurat::SpatialDimPlot(
  adc,
  images = test_image,
  group.by = spatial_label,
  pt.size.factor = 5,
  image.alpha = 0.5
)

if (length(Seurat::Images(adc)) > 0) {
  spatial_label <- if (xu_available) "GE_label_ref" else "GE_label_adc"

  for (image_name in Seurat::Images(adc)) {
    image_id = adc$sample
    safe_image_name <- stringr::str_replace_all(image_name, "[^A-Za-z0-9_-]", "_")

    p_label <- Seurat::SpatialDimPlot(
      adc,
      images = image_name,
      group.by = spatial_label,
      pt.size.factor = 5,
      image.alpha = 0.5
    ) +
      ggtitle(paste(image_id, "- dominant GE"))

    p_cluster <- Seurat::SpatialDimPlot(
      adc,
      images = image_name,
      group.by = "GE_profile_cluster",
      pt.size.factor = 5,
      image.alpha = 0.5
    ) +
      ggtitle(paste(image_id, "- GE profile cluster"))

    print(p_label + p_cluster)

    ggsave(
      file.path(
        CFG$output_dir,
        paste0("spatial_GE_states_", safe_image_name, ".pdf")
      ),
      p_label + p_cluster,
      width = 12,
      height = 6
    )
  }
}

15. How to interpret the GEs biologically

The Xu paper provides useful anchors, but these are not exclusive labels:

• GE4: cell-cycle/proliferation enrichment.
• GE2 and GE3: epithelial-mesenchymal transition features; GE3 was largely basal-like in the Xu atlas.
• GE1, GE6, GE7, and GE9: estrogen-response features; GE1/GE7 were mainly luminal and GE9 was mainly HER2-like in the Xu atlas.
• GE5 and GE10: immune/allograft-rejection and complement features; these were observed across molecular subtypes.
• GE1 and GE6: associated experimentally with resistance to NK-cell killing in the Xu study.

For an ADC resistance study, a post-treatment increase in a GE is a starting hypothesis, not automatically an ADC-resistance mechanism. The increase may reflect selection of a pre-existing state, induction by treatment, metastatic site, subtype, tumor purity, or sampling of a different lesion.

16. Recommended next analyses

16.1 Clinical metadata audit before formal inference

Before interpreting pre/post differences, confirm:

1. Which ADC each patient received.
2. Whether each biopsy was obtained before treatment, on treatment, or after progression and how many days/cycles separated treatment from biopsy.
3. Whether paired samples are from the same metastatic lesion and organ.
4. Objective response, best response, duration of response, PFS, and the exact assessment window.
5. ER/PR/HER2 status and ADC-target protein expression in the biopsied lesion.
6. Prior lines of therapy, intervening therapies, and tissue/slide batch.
7. Tumor fraction and spatial QC for every biopsy.

With 17 biopsies, avoid a single model containing timepoint, ADC, metastatic site, subtype, response, and batch simultaneously. First visualize the confounding structure. Then use one prespecified primary contrast and a small number of sensitivity analyses.

16.2 GE-stratified pseudobulk differential expression

After identifying a reproducible post-ADC GE state, perform differential expression using raw RNA counts aggregated by sample, preferably within the same GE state. This separates a change in state abundance from a transcriptional change inside that state. Do not use integrated expression values and do not treat spots as biological replicates.

Useful resistance themes include:

• ADC target expression and internalization;
• lysosomal trafficking and endosomal processing;
• payload mechanism, DNA-damage response, apoptosis, and drug efflux;
• EMT, stem-like states, proliferation, interferon response, and antigen presentation;
• stress, hypoxia, and metabolic adaptation.

These should be tested as prespecified pathway scores or sample-level pseudobulk gene-set enrichment, followed by gene-level inspection.

16.3 Spatial selection versus spatial induction

For matched cases, compare whether a post-ADC GE state:

• was already present as a rare spatial region in the pretreatment biopsy;
• expands after therapy;
• appears near immune-excluded, necrotic, hypoxic, or stromal-rich regions;
• forms spatially coherent patches rather than isolated noisy spots.

Moran’s I, neighborhood enrichment, boundary analysis, and distance-to-cell-type maps are appropriate next steps. These should use the original Visium coordinates and tumor masks.

16.4 Cell-cell interaction analysis

The Xu atlas connected GE states to immune interactions. In the ADC cohort, test whether post-treatment GE1/GE6 regions show altered NK-cell proximity or inhibitory ligand-receptor patterns, and whether GE5 regions co-localize with CD8 T-cell signals. Visium deconvolution is required before assigning sender and receiver cell types. Validate the most important findings in Xenium, where single-cell distances are observed directly.

16.5 Perturbation analysis

Use the observational results to nominate perturbations rather than training a large prediction model on 17 biopsies. A practical sequence is:

1. identify post-ADC-enriched or conserved GE states;
2. identify within-state pseudobulk genes and pathways;
3. prioritize regulators with pathway, ligand-receptor, and target-mechanism support;
4. compare predicted perturbation direction with public breast cancer cell-line or organoid perturbation data;
5. test a short list in ADC-sensitive and ADC-resistant models;
6. ask whether perturbation reverses the GE state and restores ADC sensitivity.

16.6 Cross-modality validation

Use matched Xenium or scRNA-seq to determine whether Visium GE mixtures correspond to true cancer-cell states. Recommended validation includes:

• the same GE scoring and dominant-GE assignment at single-cell resolution;
• concordance of sample-level GE composition between Visium and Xenium;
• mapping Visium out-of-reference regions to Xenium cancer-cell clusters;
• protein validation of the ADC target and selected state markers by IHC/IF;
• validation across ADC drugs and metastatic sites rather than pooling them blindly.

17. Save compact analysis tables and session information

The full spot-level table can be large, so this notebook saves only identifiers, GE scores, labels, and cluster results in a compressed CSV.

spot_output_columns <- c(
  "spot_id",
  "sample_id",
  GE_RAW_COLS,
  GE_Z_ADC_COLS,
  "GE_label_adc",
  "GE_label_adc_margin",
  "GE_profile_cluster",
  "GE_profile_cluster_distance"
)

if (xu_available) {
  spot_output_columns <- c(
    spot_output_columns,
    GE_Z_REF_COLS,
    "GE_label_ref",
    "GE_label_ref_margin",
    "GE_nearest_xu_centroid",
    "GE_xu_centroid_distance",
    "GE_out_of_xu_reference"
  )
}

readr::write_csv(
  adc_tumor |>
    dplyr::select(dplyr::all_of(intersect(spot_output_columns, colnames(adc_tumor)))),
  file.path(CFG$output_dir, "ADC_tumor_spot_GE_results.csv.gz")
)

writeLines(
  capture.output(sessionInfo()),
  file.path(CFG$output_dir, "sessionInfo.txt")
)

18. Decision rules for the next notebook

Carry a GE state forward only if several forms of evidence agree:

1. a sample-level pre/post effect, not merely a spot-level p-value;
2. a consistent direction in at least two of the three matched pairs, or a clear case-specific biological rationale;
3. no obvious explanation by metastatic site, ADC drug, tumor fraction, or technical batch;
4. spatial coherence on the tissue section;
5. a plausible pathway or target mechanism supported by raw-count pseudobulk;
6. validation in Xenium/scRNA-seq or protein imaging when available.

The immediate next input needed is the completed clinical metadata table. The most consequential fields are timepoint, patient_id, adc_drug, metastatic_site, biopsy timing relative to ADC exposure, receptor/target status, and a prespecified response definition.