<!DOCTYPE html>
ADC Cohort — Blinded Functional State Freeze
Chan Lab
r Sys.Date()
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)
XENIUM_FILE <- file.path(“01_ADC_Audit_Freeze”,“objects”,“ADC_Xenium_Blinded_Cell_Identity_Frozen.rds”)
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.
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(
———————————————————-
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 availabilityNot 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 <- 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)
}
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@meta.data )[ 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 xenium@meta.data[[score_name]] <- NA_real_
scores <- score_program( object = xenium, genes = program_gene_sets[[program_name]], cells = malignant_cells )
xenium@meta.data[ names(scores), score_name ] <- scores }
#Confirm creation malignant_score_columns <- paste0( malignant_program_names, “_score” )
malignant_score_columns
head( xenium@meta.data[ 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
)
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 <- xenium@meta.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( xenium@meta.data[ 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)
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@meta.data )[ 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 )
xenium@meta.data[ names(cytotoxic_scores), “Cytotoxicity_score”] <- cytotoxic_scores
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) <- “cell_type_final”</p> <p>p_cytotoxicity <- 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 )
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@meta.data )[ 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 )
xenium@meta.data[ 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(xenium@meta.data) ) {
t_inhibitory_summary <- xenium@meta.data %>% 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 = "drop" )write.csv( t_inhibitory_summary, file.path( TABLE_DIR, “Tcell_inhibitory_score_by_sample.csv” ), row.names = FALSE ) }
print(t_inhibitory_summary)
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@meta.data )[ 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 )
xenium@meta.data[ names(myeloid_scores), “Myeloid_Suppressive_score” ] <- myeloid_scores }
Summarize.```{r} if ( “Myeloid_Suppressive_score” %in% colnames(xenium@meta.data) ) {
myeloid_summary <- xenium@meta.data %>% 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 = "drop" )write.csv( myeloid_summary, file.path( TABLE_DIR, “myeloid_suppressive_score_by_sample.csv” ), row.names = FALSE ) }
print(myeloid_summary)
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 <- xenium@meta.data %>% 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 <- xenium\)sample[ rownames(target_expression) ]
target_expression\(cell_type_final <- 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
)
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 metadataSaving 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 <- xenium@meta.data %>% select( any_of( functional_metadata_columns ) )
write.csv( functional_state_metadata, file.path( TABLE_DIR, “functional_state_metadata.csv” ), row.names = TRUE )
:::
::::
:::::
::::::
:::::::
::::::::
:::::::::
::::::::::
:::::::::::
::::::::::::
:::::::::::::::::::::
::::::::::::::::::::::::
:::::::::::::::::::::::::