Spatial multi-omics analysis of tumor-stroma boundary cell features for predicting breast cancer progression and therapy response

Abstract

Background: The tumor boundary of breast cancer represents a highly heterogeneous region. In this area, the interactions between malignant and non-malignant cells influence tumor progression, immune evasion, and drug resistance. However, the spatial transcriptional profile of the tumor boundary and its role in the prognosis and treatment response of breast cancer remain unclear.

Method: Utilizing the Cottrazm algorithm, we reconstructed the intricate boundaries and identified differentially expressed genes (DEGs) associated with these regions. Cell-cell co-positioning analysis was conducted using SpaCET, which revealed key interactions between tumor-associated macrophage (TAMs) and cancer-associated fibroblasts (CAFs). Additionally, Lasso regression analysis was employed to develop a malignant body signature (MBS), which was subsequently validated using the TCGA dataset for prognosis prediction and treatment response assessment.

Results: Our research indicates that the tumor boundary is characterized by a rich reconstruction of the extracellular matrix (ECM), immunomodulatory regulation, and the epithelial-to-mesenchymal transition (EMT), underscoring its significance in tumor progression. Spatial colocalization analysis reveals a significant interaction between CAFs and M2-like tumor-associated macrophage (TAM), which contributes to immune exclusion and drug resistance. The MBS score effectively stratifies patients into high-risk groups, with survival outcomes for patients exhibiting high MBS scores being significantly poorer. Furthermore, drug sensitivity analysis demonstrates that high-MB tumors had poor response to chemotherapy strategies, highlighting the role of the tumor boundary in modulating therapeutic efficacy.

Conclusion: Collectively, we investigate the spatial transcription group and bulk data to elucidate the characteristics of tumor boundary molecules in breast cancer. The CAF-M2 phenotype emerges as a critical determinant of immunosuppression and drug resistance, suggesting that targeting this interaction may improve treatment responses. Furthermore, the MBS serves as a novel prognostic tool and offers potential strategies for guiding personalized treatment approaches in breast cancer.

Keywords: CAF-M2 interaction; breast cancer; prognostic model; spatial transcriptomics; therapy resistance; tumor boundary.

FIGURE 1

Breast cancer tumor boundary reconstruction. (A) The Cottrazm algorithm is employed to determine the spatial distribution of malignant (Mal), boundary (Bdy), and non-malignant (nMal) regions across four breast cancer datasets. (B) Venn diagram illustrates the overlap of differentially expressed genes (DEGs) among these datasets. The genes are categorized into primary functional groups, including extracellular matrix (ECM) genes (e.g., COL5A2, COL5A1, COL4A2), immune-related genes (e.g., CCL5, CAVIN1, CALD1), migration and proliferation genes (e.g., VIM, TAGLN, FN1), epithelial tumor markers (e.g., EPCAM, KRT18, KRT19), smooth muscle-related genes (e.g., VIM, ACTA2, TAGLN), and macrophage-related genes (e.g., MRC2, MFGE8, LSP1). (C) Spatial scores of key gene characteristics in breast cancer samples are presented, with each row representing different breast cancer samples (10x-BRCA, 10X-BRCA2, GSM433610, GSM6177603) and each column corresponding to specific gene categories.

FIGURE 2

The analysis of the genetic functions of tumor boundary high-expression genes. (A) Gene Ontology (GO) enrichment analysis categorizes the border-related genes into three primary functional groups: biological processes (BP), cellular components (CC), and molecular functions (MF). (B) The analysis of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways reveals significant enrichment within specific pathways. (C) Examination of the Hallmark gene collection identifies five principal pathways.

FIGURE 3

Analysis of spatial transcription groups across various breast cancer subtypes. (A) Spatial clustering of breast cancer tissue samples, including ductal carcinoma in situ (DCIS), infiltrating ductal carcinoma (IDC), infiltrating lobular carcinoma (ILC), and triple-negative breast cancer (TNBC). Each point represents a spatial transcription group. (B) The spatial distribution of the C0 cluster. (C) Comparison of boundary gene expression between cluster C0 and other clusters. Each point represents a gene, indicating the proportion of the gene’s presence, with color denoting the average expression level. (D) DEGs among 12 clusters identified in the spatial transcription group data. Point size reflects the proportion of points corresponding to a given gene, while color intensity represents the level of gene expression. (E) Spatial localization of cancer-associated fibroblasts (CAFs) and tumor-associated macrophage (TAMs). The left panel for each subtype presented the CAF score distribution, the middle panel shows the TAM score distribution, and the right panel illustrates the combined positioning (CAF_TAM) of CAF and TAM scores, highlighting potential interactions within the tumor microenvironment.

FIGURE 4

Cell-cell colocalization analysis and interaction analysis. (A) Cell to cell colocalization analysis using SpaCET illustrates the spatial relationships among different cell types within the tumor microenvironment. Each point represents a pairwise co-localization of these cell types, with the Rho value (Spearman’s rank correlation coefficient) indicating the strength of the co-localization. Larger dots correspond to higher fraction products, indicating more frequent colocalization between cell types. The color scale represents Rho values, where red and blue indicate positive and negative correlations, respectively. Statistical significance of co-localizations was assessed using Spearman’s correlation. (B) The correlation between cell scores and reference cells is depicted, where each point represents a specific cell type. This representation shows the correlation between estimated cell scores (Y-axis) and reference profiles (X-axis). The size of each point corresponds to the score product, and the shaded area represents the confidence range of the model’s output. (C) Verification through multiple immunofluorescence staining demonstrates the co-localization of cancer-associated fibroblasts (CAF) and macrophages within the tumor microenvironment. The markers used include DAPI (nucleus, blue), ACTA2 (CAF, purple), FAP (fibroblast activation protein, green), CD68 (all macrophages, yellow), and CD163 (M2 macrophage marker, orange). The magnified area highlights the spatial interaction between CAF and M2 macrophages.

FIGURE 5

The development and verification of breast cancer prognosis is presented herein. (A) The Lasso (Least Absolute Shrinkage and Selection Operator) regression model is utilized for feature selection. The figure illustrates the relationship between the LAMBDA value and the number conversion, with the vertical dotted line indicating the optimal LAMBDA value determined through 10-fold cross-validation, corresponding to the minimum deviation. (B) The Lasso coefficient curve of the candidate prognostic genes is depicted, where each colored line represents a gene; the vertical dotted line denotes the selected LAMBDA value retained in the final model. (C) The genes and their coefficients used to build the model are showed. (D) Kaplan-Meier survival analysis is conducted based on the characteristics of the prognostic genes, comparing high-risk and low-risk groups using a median threshold value. (E) The clinical characteristics of patients within the risk groups are summarized. The heat map illustrates the overall survival period (OS), PAM50 molecular subtype, pathology, AJCC staging, age, and gender distribution among high-risk and low-risk groups, indicating a correlation between risk scores and clinical pathological characteristics.

FIGURE 6

The prognosis of risk characteristics across various breast cancer subtypes is analyzed. (A) Kaplan-Meier survival analysis of different PAM50 molecular subtypes of breast cancer, including Basal, HER2⁺, Luminal A (LumA), and Luminal B (LumB), reveals that patients are categorized into high-risk and low-risk groups based on risk scores. Each panel presents the hazard ratio (HR), 95% confidence interval (CI), and P-values from the risk assessment. Notable differences in survival rates are observed between the LumA and LumB subtypes. (B) Kaplan-Meier survival analysis for breast cancer subtypes, specifically IDC and ILC, indicates that high-risk patients in both subtypes exhibit significantly poorer outcomes. (C) A schematic representation illustrates the relationship between breast cancer pathology, PAM50 molecular subtypes, and risk groups. The classification of patients into high-risk and low-risk groups based on different pathologies and molecular subtypes suggests a correlation with their risk characteristics. (D) Kaplan-Meier survival analysis compares high-risk and low-risk groups for both PAM50 subtypes (left) and pathological subtypes (right). The survival curves demonstrate that the risk characteristics effectively stratify patients into distinct subtypes, with high-risk patients consistently exhibiting poor prognoses.

FIGURE 7

Chemotherapy sensitivity analysis for Malignant-boundary signature. (A) Comparison of the estimated effects of various chemotherapy drugs between the MBS-high group and the MBS-low group. The drugs analyzed include Oxaliplatin, Erlotinib, Bortezomib, Mitomycin C, 5-Fluorouracil, Vinorelbine, Gemcitabine and Pyrimethamine. We utilized the Wilcoxon test to assess statistical significance, displaying the P values in each panel. A lower IC50 value indicates higher drug sensitivity. (B) Comparison of IC50 of Dasatinib, Thapsigargin, and Sorafenib between MBS-high and MBS-low groups.

FIGURE 8

Tumor boundary characteristics and immune response in breast cancer. (A) Distribution of TIDE scores across all patients. Based on TIDE values, responders (blue) and non-responders (red) are clearly distinguished, indicating the predictive value of TIDE in assessing immune evasion in cancer treatment. (B) Proportion of responders and non-responders in low and high TIDE score groups. Chi-square analysis was used for statistical analysis. (C) Comparison of TIDE scores between the low and high MBS groups (* represents P < 0.05; ** represents P < 0.01; *** represents P < 0.001, **** represents P < 0.0001). (D) Comparison of expression profiles of various immune-related genes and pathways in the low and high signature groups. The box plots show the median and interquartile range (* represents P < 0.05; ** represents P < 0.01; *** represents P < 0.001, **** represents P < 0.0001).