<!DOCTYPE html>

ADC Cohort — Blinded Functional State Freeze

Purpose

This notebook defines and freezes blinded functional cell-state features in the metastatic breast cancer ADC Xenium cohort after broad cell identities were validated and frozen in Notebook 1.

The objectives are to:

Load the frozen cell-identity object generated in Notebook 1. Define biologically interpretable malignant-cell programs measurable within the targeted Xenium panel. Define cytotoxic and immunoregulatory programs within validated immune populations. Quantify functional-state scores by cell and specimen without using treatment or clinical outcome information. Evaluate whether state scores are driven disproportionately by individual specimens. Freeze continuous functional-state features for later spatial neighborhood analysis.

No ADC identity, pre/post-treatment status, response, progression, survival, or other clinical outcome information is used in this notebook.

Functional programs will primarily be retained as continuous scores. Discrete cell-state labels will only be created when marker support and cell abundance are sufficient.

```{r} # ============================================================ # 01. SETUP # ============================================================ library(Seurat) library(dplyr) library(tidyr) library(ggplot2) library(patchwork) set.seed(1234)

————————————————————

Input

————————————————————

XENIUM_FILE <- file.path(“01_ADC_Audit_Freeze”,“objects”,“ADC_Xenium_Blinded_Cell_Identity_Frozen.rds”)

————————————————————

Output

————————————————————

OUTPUT_DIR <- “02_ADC_Blinded_Functional_State_Freeze”

FIGURE_DIR <- file.path(OUTPUT_DIR, “figures”) TABLE_DIR <- file.path(OUTPUT_DIR, “tables”) OBJECT_DIR <- file.path(OUTPUT_DIR, “objects”)

dir.create(OUTPUT_DIR, recursive = TRUE, showWarnings = FALSE) dir.create(FIGURE_DIR, showWarnings = FALSE) dir.create(TABLE_DIR, showWarnings = FALSE) dir.create(OBJECT_DIR, showWarnings = FALSE)


#Load frozen cell identities

```{r}
# ============================================================
# 02. LOAD FROZEN OBJECT
# ============================================================

xenium <- readRDS("/work/InternalMedicine/s239947/ADC/01_ADC_Audit_Freeze/objects/ADC_Xenium_Blinded_Cell_Identity_Frozen.rds")

xenium

stopifnot(
  "cell_type_final" %in% colnames(xenium@meta.data),
  "broad_lineage" %in% colnames(xenium@meta.data)
)

DefaultAssay(xenium) <- "SCT"

table(xenium$cell_type_final)
table(xenium$broad_lineage)

The broad cell identities established in Notebook 1 are treated as frozen and will not be modified in this notebook.

  1. Define candidate functional programs

The candidate programs below were selected to represent biologically relevant malignant and immune phenotypes potentially relevant to ADC response. Their relationship to treatment response will not be examined until after blinded feature definition is complete.

```{r} # ============================================================ # 03. CANDIDATE FUNCTIONAL PROGRAMS # ============================================================

candidate_programs <- list(

———————————————————- # MALIGNANT-CELL PROGRAMS

———————————————————-

Proliferation = c( “MKI67”, “TOP2A”, “PCNA”, “CCNB1”, “CCNA2”, “CDK1”, “UBE2C” ),

EMT_Plasticity = c( “VIM”, “ZEB1”, “ZEB2”, “SNAI1”, “SNAI2”, “TWIST1”, “FN1”, “ITGA5” ),

Interferon_Response = c( “STAT1”, “IRF1”, “ISG15”, “IFIT1”, “IFIT3”, “MX1”, “CXCL10” ),

Antigen_Presentation = c( “B2M”, “HLA-A”, “HLA-B”, “HLA-C”, “TAP1”, “TAP2” ),

DNA_Damage_Stress = c( “TP53”, “CDKN1A”, “GADD45A”, “ATM”, “ATR”, “CHEK1”, “CHEK2”, “PARP1” ),

Survival_Apoptosis = c( “BCL2”, “BCL2L1”, “MCL1”, “BAX”, “CASP3” ),

Hypoxia_Stress = c( “HIF1A”, “VEGFA”, “CA9”, “JUN”, “FOS”, “ATF3” ),

# ———————————————————- # IMMUNE PROGRAMS # ———————————————————-

Cytotoxicity = c( “NKG7”, “GNLY”, “PRF1”, “GZMB” ),

T_Cell_Inhibitory = c( “PDCD1”, “LAG3”, “TIGIT”, “HAVCR2” ),

Myeloid_Suppressive = c( “CD163”, “MRC1”, “ARG1”, “IL10”, “FCER1G”, “TYROBP” ) )

# 4. Audit gene availability

Not all genes in each candidate program are expected to be present in the targeted Xenium panel. We therefore explicitly record which genes are available before calculating any scores.

```{r} # ============================================================ # 04. GENE AVAILABILITY AUDIT # ============================================================

program_availability <- bind_rows( lapply( names(candidate_programs), function(program_name) {

  genes &lt;- candidate_programs[[program_name]]

  data.frame(
    program = program_name,
    gene = genes,
    available = genes %in% rownames(xenium)
  )
}

) )

program_availability

write.csv( program_availability, file.path( TABLE_DIR, “functional_program_gene_availability.csv” ), row.names = FALSE )

Summarize gene availability for each program.

```{r} program_availability_summary <- program_availability %>% group_by(program) %>% summarise( genes_requested = n(), genes_available = sum(available), fraction_available = mean(available), available_genes = paste( gene[available], collapse = “;” ), .groups = “drop” )

program_availability_summary

write.csv( program_availability_summary, file.path( TABLE_DIR, “functional_program_availability_summary.csv” ), row.names = FALSE )


# 5. Define supported programs

A multigene functional program will be retained only when at least three measured genes are available.

Single genes with specific mechanistic importance, such as `FCGR3A`, `TACSTD2`, and `ERBB2`, will be retained separately rather than described as pathway scores.

```{r}
# ============================================================
# 05. DEFINE SUPPORTED PROGRAMS
# ============================================================

MIN_PROGRAM_GENES <- 3

supported_programs <- program_availability_summary %>%
  filter(
    genes_available >= MIN_PROGRAM_GENES
  ) %>%
  pull(program)

supported_programs

Construct final measured gene sets.

```{r} program_gene_sets <- lapply( supported_programs, function(program_name) {

program_availability %>%
  filter(
    program == program_name,
    available
  ) %>%
  pull(gene)

} )

names(program_gene_sets) <- supported_programs

program_gene_sets

#Save exactly which genes will contribute to the final scores final_program_gene_table <- bind_rows( lapply( names(program_gene_sets), function(program_name) {

  data.frame(
    program = program_name,
    gene = program_gene_sets[[program_name]]
  )
}

) )

write.csv( final_program_gene_table, file.path( TABLE_DIR, “final_functional_program_genes.csv” ), row.names = FALSE )


# 6. Define the targeted-panel scoring function

Each functional program is calculated as the mean scaled expression of its measured component genes.

For a cell \(i\) and program containing \(G\) measured genes:

\[
S_i = \frac{1}{G}\sum_{g=1}^{G} z_{ig}
\]

where \(z_{ig}\) represents scaled expression of gene \(g\) in cell \(i\).

```{r}
# ============================================================
# 06. TARGETED-PANEL SCORING FUNCTION
# ============================================================

score_program <- function(
  object,
  genes,
  cells = NULL,
  assay = "SCT",
  layer = "scale.data",
  min_genes = 3
) {

  genes_present <- intersect(
    genes,
    rownames(object[[assay]])
  )

  if (length(genes_present) < min_genes) {

    stop(
      paste0(
        "Insufficient genes for scoring: ",
        length(genes_present),
        " genes available; minimum = ",
        min_genes
      )
    )
  }

  if (is.null(cells)) {
    cells <- colnames(object)
  }

  mat <- LayerData(
    object = object,
    assay = assay,
    layer = layer,
    features = genes_present,
    cells = cells
  )

  if (inherits(mat, "Matrix")) {
    scores <- Matrix::colMeans(mat)
  } else {
    scores <- colMeans(mat)
  }

  return(scores)
}

  1. Score malignant-cell programs

Functional malignant-cell scores will only be assigned to cells previously classified as Cancer Epithelial Cells.

Other cell types will receive NA for malignant-specific scores.

```{r} # ============================================================ # 07. MALIGNANT-CELL PROGRAMS # ============================================================

malignant_program_names <- intersect( c( “Proliferation”, “EMT_Plasticity”, “Interferon_Response”, “Antigen_Presentation”, “DNA_Damage_Stress”, “Survival_Apoptosis”, “Hypoxia_Stress” ), names(program_gene_sets) )

malignant_program_names

Identify malignant cells.

```{r} malignant_cells <- rownames( )[ xenium$cell_type_final == “Cancer Epithelial Cells”]

length(malignant_cells)

Score each supported malignant program.

```{r} for (program_name in malignant_program_names) {

score_name <- paste0( program_name, “_score” )

# Initialize as NA for all cells [[score_name]] <- NA_real_

scores <- score_program( object = xenium, genes = program_gene_sets[[program_name]], cells = malignant_cells )

[ names(scores), score_name ] <- scores }

#Confirm creation malignant_score_columns <- paste0( malignant_program_names, “_score” )

malignant_score_columns

head( [ malignant_cells, c( “sample”, “cell_type_final”, malignant_score_columns ), drop = FALSE ] )


# 8. Summarize malignant programs by specimen

These summaries are descriptive only. No sample is currently classified by ADC, response, or treatment timing.

```{r}
# ============================================================
# 08. MALIGNANT STATE SUMMARY
# ============================================================

malignant_state_summary <- xenium@meta.data %>%
  filter(
    cell_type_final ==
      "Cancer Epithelial Cells"
  ) %>%
  group_by(sample) %>%
  summarise(
    cancer_cells = n(),

    across(
      all_of(malignant_score_columns),

      list(
        median = ~median(
          .x,
          na.rm = TRUE
        ),

        mean = ~mean(
          .x,
          na.rm = TRUE
        )
      )
    ),

    .groups = "drop"
  )

malignant_state_summary

#Save
write.csv(
  malignant_state_summary,
  file.path(
    TABLE_DIR,
    "malignant_state_summary_by_sample.csv"
  ),
  row.names = FALSE
)

  1. Visualize specimen-level malignant states

For visualization, use specimen-level median scores rather than plotting hundreds of thousands of individual malignant cells.

```{r} # ============================================================ # 09. MALIGNANT STATE HEATMAP # ============================================================

malignant_heatmap_data <- %>% filter( cell_type_final == “Cancer Epithelial Cells” ) %>% group_by(sample) %>% summarise( across( all_of(malignant_score_columns), ~median(.x, na.rm = TRUE) ), .groups = “drop” ) %>% pivot_longer( -sample, names_to = “program”, values_to = “median_score” )

p_malignant_heatmap <- ggplot( malignant_heatmap_data, aes( x = program, y = sample, fill = median_score ) ) + geom_tile() + labs( title = “Blinded malignant functional programs by specimen”, x = NULL, y = “Specimen”, fill = “Medianscore” ) + theme_minimal() + theme( axis.text.x = element_text( angle = 45, hjust = 1 ) )

p_malignant_heatmap

#Save ggsave( file.path( FIGURE_DIR, “Heatmap_malignant_programs_by_sample.png” ), p_malignant_heatmap, width = 10, height = 7, dpi = 300 )

ggsave( file.path( FIGURE_DIR, “Heatmap_malignant_programs_by_sample.pdf” ), p_malignant_heatmap, width = 10, height = 7 )

# 10. Evaluate correlation among malignant programs

```{r} # ============================================================ # 10. MALIGNANT PROGRAM CORRELATIONS # ============================================================

if (length(malignant_score_columns) >= 2) {

malignant_cor <- cor( [ malignant_cells, malignant_score_columns, drop = FALSE ], use = “pairwise.complete.obs”, method = “spearman” )

write.csv( malignant_cor, file.path( TABLE_DIR, “malignant_program_correlations.csv” ) )

malignant_cor_long <- as.data.frame( malignant_cor ) %>% tibble::rownames_to_column( “program_1” ) %>% pivot_longer( -program_1, names_to = “program_2”, values_to = “rho” )

p_malignant_cor <- ggplot( malignant_cor_long, aes( x = program_1, y = program_2, fill = rho ) ) + geom_tile() + labs( title = “Correlation among malignant functional programs”, x = NULL, y = NULL, fill = “Spearman” ) + theme_minimal() + theme( axis.text.x = element_text( angle = 45, hjust = 1 ) )

ggsave( file.path( FIGURE_DIR, “Heatmap_malignant_program_correlations.png” ), p_malignant_cor, width = 8, height = 7, dpi = 300 ) }

print(p_malignant_cor)

  1. Score cytotoxic lymphocyte function

Notebook 1 validated NK cells and CD8 T cells as distinct broad identities. Because both populations express cytotoxic machinery, the same continuous cytotoxic program will be calculated within both populations.

NK cells will not be subclustered because of their limited abundance.

```{r} # ============================================================ # 11. CYTOTOXIC LYMPHOCYTE FUNCTION # ============================================================

cytotoxic_cells <- rownames( )[ xenium\(cell_type_final %in% c( “NK Cells”, “CD8+ T Cells” )]</p> <p>length(cytotoxic_cells)</p> <p>#Confirm cytotoxic program is supported “Cytotoxicity” %in% names(program_gene_sets)</p> <p>program_gene_sets\)Cytotoxicity

Calculate the score.

```{r} xenium$Cytotoxicity_score <- NA_real_

cytotoxic_scores <- score_program( object = xenium, genes = program_gene_sets$Cytotoxicity, cells = cytotoxic_cells )

[ names(cytotoxic_scores), “Cytotoxicity_score”] <- cytotoxic_scores

  1. Visualize NK and CD8 functional differences

These plots verify that the intended measurements behave biologically as expected. They are not treatment-response comparisons.

{r} Idents(xenium) &lt;- “cell_type_final”</p> <p>p_cytotoxicity &lt;- VlnPlot( xenium, features = “Cytotoxicity_score”, idents = c( “NK Cells”, “CD8+ T Cells” ), pt.size = 0.3 ) + ggtitle( “Cytotoxicity program in NK and CD8 T cells” )</p> <p>p_cytotoxicity</p> <pre><code>{r} p_fcgr3a <- VlnPlot( xenium, features = “FCGR3A”, idents = c( “NK Cells”, “CD8+ T Cells” ), pt.size = 0.3 ) + ggtitle( “FCGR3A/CD16 expression in cytotoxic lymphocytes” )

p_fcgr3a

#Save plots ggsave( file.path( FIGURE_DIR, “VlnPlot_NK_CD8_cytotoxicity.png” ), p_cytotoxicity, width = 7, height = 6, dpi = 300 )

ggsave( file.path( FIGURE_DIR, “VlnPlot_NK_CD8_FCGR3A.png” ), p_fcgr3a, width = 7, height = 6, dpi = 300 )

  1. Define the T-cell inhibitory-receptor program

Expression of inhibitory receptors does not by itself establish biological T-cell exhaustion. Therefore, this feature is conservatively termed the T_Cell_Inhibitory_score.

```{r} # ============================================================ # 15. T-CELL INHIBITORY PROGRAM # ============================================================

if ( “T_Cell_Inhibitory” %in% names(program_gene_sets) ) {

t_cells <- rownames( )[ xenium\(cell_type_final %in% c( “CD8+ T Cells”, “CD4+ T Cells”, “Regulatory T Cells” ) ]</p> <p>xenium\)T_Cell_Inhibitory_score <- NA_real_

t_inhibitory_scores <- score_program( object = xenium, genes = program_gene_sets$T_Cell_Inhibitory, cells = t_cells )

[ names(t_inhibitory_scores), “T_Cell_Inhibitory_score” ] <- t_inhibitory_scores }

Summarize if the score was supported.

```{r} if ( “T_Cell_Inhibitory_score” %in% colnames() ) {

t_inhibitory_summary <- %>% filter( cell_type_final %in% c( “CD8+ T Cells”, “CD4+ T Cells”, “Regulatory T Cells” ) ) %>% group_by( sample, cell_type_final ) %>% summarise( cells = n(),

  median_inhibitory_score =
    median(
      T_Cell_Inhibitory_score,
      na.rm = TRUE
    ),

  .groups = &quot;drop&quot;
)

write.csv( t_inhibitory_summary, file.path( TABLE_DIR, “Tcell_inhibitory_score_by_sample.csv” ), row.names = FALSE ) }

print(t_inhibitory_summary)

  1. Define suppressive myeloid function

Macrophages, monocytes, and cells previously classified as MDSCs will be evaluated using a continuous suppressive-myeloid program if at least three relevant genes are measured.

The score does not itself establish functional MDSC identity.

```{r} # ============================================================ # 16. MYELOID SUPPRESSIVE PROGRAM # ============================================================

if ( “Myeloid_Suppressive” %in% names(program_gene_sets) ) {

suppressive_myeloid_cells <- rownames( )[ xenium\(cell_type_final %in% c( “Macrophages”, “Monocytes”, “MDSCs” ) ]</p> <p>xenium\)Myeloid_Suppressive_score <- NA_real_

myeloid_scores <- score_program( object = xenium, genes = program_gene_sets$Myeloid_Suppressive, cells = suppressive_myeloid_cells )

[ names(myeloid_scores), “Myeloid_Suppressive_score” ] <- myeloid_scores }

Summarize.

```{r} if ( “Myeloid_Suppressive_score” %in% colnames() ) {

myeloid_summary <- %>% filter( cell_type_final %in% c( “Macrophages”, “Monocytes”, “MDSCs” ) ) %>% group_by( sample, cell_type_final ) %>% summarise( cells = n(),

  median_suppressive_score =
    median(
      Myeloid_Suppressive_score,
      na.rm = TRUE
    ),

  .groups = &quot;drop&quot;
)

write.csv( myeloid_summary, file.path( TABLE_DIR, “myeloid_suppressive_score_by_sample.csv” ), row.names = FALSE ) }

print(myeloid_summary)

  1. Record Treg abundance without further subclustering

Because regulatory T cells are sparse, the existing Treg identity will be retained without attempting to define multiple Treg subtypes.

Their abundance and spatial relationship to tumor cells and cytotoxic lymphocytes will be evaluated in Notebook 3.

```{r} # ============================================================ # 17. TREG AUDIT # ============================================================

treg_summary <- %>% group_by(sample) %>% summarise( total_cells = n(),

Treg_cells = sum( cell_type_final == “Regulatory T Cells” ),

Treg_fraction = Treg_cells / total_cells,

.groups = “drop”

)

treg_summary

#Save write.csv( treg_summary, file.path( TABLE_DIR, “Treg_abundance_by_sample.csv” ), row.names = FALSE )


# 18. Retain ADC-target expression as blinded tumor features

The Xenium panel contains both `TACSTD2`, encoding Trop-2, and `ERBB2`, encoding HER2.

These measurements are retained as blinded malignant-cell features for later testing after ADC identity is revealed.

```{r}
# ============================================================
# 18. ADC TARGET FEATURES
# ============================================================

target_genes <- intersect(
  c(
    "TACSTD2",
    "ERBB2"
  ),
  rownames(xenium)
)

target_genes

Summarize expression among cancer epithelial cells.

```{r} target_expression <- FetchData( xenium, vars = target_genes, layer = “data” )

target_expression\(sample &lt;- xenium\)sample[ rownames(target_expression) ]

target_expression\(cell_type_final &lt;- xenium\)cell_type_final[ rownames(target_expression)]

target_summary <- target_expression %>% filter( cell_type_final == “Cancer Epithelial Cells” ) %>% group_by(sample) %>% summarise( across( all_of(target_genes),

  list(
    median = ~median(
      .x,
      na.rm = TRUE
    ),

    mean = ~mean(
      .x,
      na.rm = TRUE
    ),

    fraction_positive =
      ~mean(
        .x > 0,
        na.rm = TRUE
      )
  )
),

.groups = "drop"

)

target_summary

#Save write.csv( target_summary, file.path( TABLE_DIR, “ADC_target_expression_by_sample.csv” ), row.names = FALSE )


# 19. Generate specimen-level functional feature table

The patient/specimen—not the individual cell—will ultimately serve as the biological unit for clinical association testing.

We therefore generate blinded specimen-level summaries now.

```{r}
# ============================================================
# 19. SPECIMEN-LEVEL FUNCTIONAL FEATURES
# ============================================================

sample_cell_composition <- xenium@meta.data %>%
  group_by(sample) %>%
  summarise(
    total_cells = n(),

    cancer_cells =
      sum(
        cell_type_final ==
          "Cancer Epithelial Cells"
      ),

    NK_cells =
      sum(
        cell_type_final ==
          "NK Cells"
      ),

    CD8_cells =
      sum(
        cell_type_final ==
          "CD8+ T Cells"
      ),

    Treg_cells =
      sum(
        cell_type_final ==
          "Regulatory T Cells"
      ),

    macrophage_cells =
      sum(
        cell_type_final ==
          "Macrophages"
      ),

    MDSC_cells =
      sum(
        cell_type_final ==
          "MDSCs"
      ),

    fibroblast_cells =
      sum(
        cell_type_final ==
          "Fibroblasts"
      ),

    .groups = "drop"
  )

Merge currently available functional summaries.

{r} sample_functional_features <- sample_cell_composition %>% left_join( malignant_state_summary, by = “sample” )

Add target-expression summaries.

```{r} sample_functional_features <- sample_functional_features %>% left_join( target_summary, by = “sample” )

print(sample_functional_features)

#Save write.csv( sample_functional_features, file.path( TABLE_DIR, “blinded_sample_functional_features.csv” ), row.names = FALSE )


# 20. Functional-state freeze manifest

```{r}
# ============================================================
# 20. FREEZE MANIFEST
# ============================================================

state_score_columns <- grep(
  "_score$|FCGR3A_expression$",
  colnames(xenium@meta.data),
  value = TRUE
)

state_freeze_manifest <- data.frame(

  field = c(
    "freeze_date",
    "input_object",
    "clinical_metadata_used",
    "scoring_method",
    "minimum_genes_per_program",
    "malignant_programs",
    "NK_subclustering",
    "NK_features",
    "T_cell_feature",
    "myeloid_feature",
    "Treg_strategy",
    "ADC_target_features",
    "status"
  ),

  value = c(
    as.character(Sys.Date()),

    basename(XENIUM_FILE),

    "No",

    paste(
      "Mean scaled expression of measured genes;",
      "no AddModuleScore control-gene sampling"
    ),

    as.character(
      MIN_PROGRAM_GENES
    ),

    paste(
      malignant_program_names,
      collapse = "; "
    ),

    "Not performed because NK cells are sparse",

    paste(
      "NK identity; cytotoxicity score;",
      "FCGR3A/CD16 expression"
    ),

    ifelse(
      "T_Cell_Inhibitory_score" %in%
        colnames(xenium@meta.data),
      "Continuous inhibitory-receptor score",
      "Insufficient panel support"
    ),

    ifelse(
      "Myeloid_Suppressive_score" %in%
        colnames(xenium@meta.data),
      "Continuous suppressive-myeloid score",
      "Insufficient panel support"
    ),

    "Retain validated Treg identity; no Treg subclustering",

    paste(
      target_genes,
      collapse = "; "
    ),

    "Blinded functional features frozen"
  )
)

state_freeze_manifest

#Save
write.csv(
  state_freeze_manifest,
  file.path(
    TABLE_DIR,
    "functional_state_freeze_manifest.csv"
  ),
  row.names = FALSE
)

  1. Save frozen functional-state object

```{r} # ============================================================ # 21. SAVE FROZEN OBJECT # ============================================================

saveRDS( xenium, file = file.path( OBJECT_DIR, “ADC_Xenium_Blinded_Functional_States_Frozen.rds” ) )

# 22. Save compact functional metadata

Saving the full metadata table separately makes later auditing easier without requiring the complete Seurat object.

```{r} functional_metadata_columns <- c( “sample”, “cell_type_final”, “broad_lineage”, state_score_columns )

functional_state_metadata <- %>% select( any_of( functional_metadata_columns ) )

write.csv( functional_state_metadata, file.path( TABLE_DIR, “functional_state_metadata.csv” ), row.names = TRUE )

:::

::::

:::::

::::::

:::::::

::::::::

:::::::::

::::::::::

:::::::::::

::::::::::::

:::::::::::::::::::::

::::::::::::::::::::::::

:::::::::::::::::::::::::