This study investigates the prognostic value of matrix genes in breast invasive carcinoma (TCGA-BRCA) using data from the TCGA database. We downloaded the dataset consisting of 1,222 samples, integrated clinical data, and filtered for 1,109 patients with complete survival and TNM staging information.

Differential expression analysis was conducted using the limma package to identify differentially expressed genes (DEGs) linked to patient survival, categorizing patients into long-term (≥1 year) and short-term (<1 year) survival groups. Gene Ontology (GO) and KEGG pathway enrichment analyses were performed on these DEGs to understand their biological roles.

We developed a risk signature utilizing elastic net penalized logistic regression, optimizing the model to predict patient outcomes based on gene expression profiles. Each patient received a matrix risk score, enabling classification into high- and low-risk groups, followed by survival analysis using Kaplan–Meier curves.

To identify molecular subtypes, we implemented consensus clustering on MRDEG expression and validated results with independent datasets. We also assessed immune cell infiltration using the CIBERSORT algorithm. Pathway activity was analyzed with GSVA, focusing on hallmark pathways. Additionally, ligand-receptor interactions were examined to explore signaling dynamics in the tumor microenvironment. Statistical analyses were performed in R, with p<0.05 considered significant Figure 1 .

We downloaded the breast invasive carcinoma (BRCA) dataset, TCGA-BRCA (n=1,222), from the TCGA database (25) using the TCGAbiolinks package (26). The data type was selected as count and converted to FPKM format. In addition, we obtained clinical data corresponding to the matched samples of the TCGA BRCA dataset from the TCGA GDC ¹ official website, including age, survival status, follow-up time, and tumor stage. Our study excluded patients with no survival information and incomplete TNM staging information, and 1,109 patients were included in our subsequent analyses.

To examine the gene mutation status of TCGA-BRCA patients, we acquired ‘Masked Somatic Mutation’ data from the official website of TCGA. The obtained data served as the somatic mutation dataset for patients with breast cancer. We preprocessed the data using VarScan software and visualized the somatic mutations of patients using the maftools package (27). We obtained the body genes from the study of Naba et al. (28), which contains 1,062 matrix body genes, and the specific information is presented in Supplementary Table S1 and Table 1 .

In addition, we assessed the Molecular Signature Database (29) (MSigDB) ² . The 50 Hallmark gene sets were obtained from “h.all.v2023.1.Hs.symbols.gmt” on the database website, from “c2.cp.kegg.v 7.4.symbols.gmt” file to obtain the KEGG pathway gene set for subsequent Gene Set Enrichment Analysis (GSEA).

To further validate our approach, we acquired a set of scRNA-seq data, GSE161529, from the Tumor Immune Single-cell Hub 2 ³ database (30). Additionally, we obtained an independent validation set for breast cancer, UCSC (Caldas 2007), from the Xena platform ⁴ , and another independent validation set, GSE20685, from the Gene Expression Omnibus (GEO) database.

To elucidate the potential mechanism of gene action and related biological characteristics influencing the prognoses of patients with breast cancer, we first divided patients with an overall survival of >1 year into the long survival group according to their prognoses. Patients aged<1 year were classified into the short-term survival group. Genes in different groups were subjected to differential analysis using the limma package. Genes with an absolute value of log fold change (|logFC|) >0.25 and a p-value<0.05 were considered as differentially expressed genes (DEGs) associated with prognosis in patients with breast cancer. To obtain breast cancer-related matrix body-related DEGs (microsomal-related (MR) DEGs), we compared the DEGs obtained from the differential analysis of TCGA-BRCA datasets with matrix body-related genes (microsomal-related genes, MRGs) at the intersection and drew a Venn diagram. The findings from the differential analysis were visualized using the R package ggplot2 to generate a volcano plot and the R package pheatmap to create a heatmap.

GO analysis (31) is a standard method for large-scale functional enrichment studies, including biological processes, cell components, and molecular functions. KEGG (32) is a database containing information on genomes, biological pathways, diseases, and drugs. We used R-package clusterProfiler (33) to analyze the GO and pathway (KEGG) enrichment of differentially expressed matrix genes, and the screening criteria for entries were adj. p<0.05. An FDR value (q-value)<0.25 was considered statistically significant, and the correction method of p-value was Benjamini–Hochberg.

We performed risk signature selection on matrix genes differentially expressed between long- and short-term survival in the TCGA-BRCA dataset. We used Elastic Net penalized logistic regression to select highly correlated differentially expressed matrix genes based on the correlation between long and short survival. The elastic net penalized logistic regression was implemented using the glmfit function in the R package glmnet, where the parameter alpha was set to 0.5. Using α=0.5 in penalized logistic regression is to combine the advantages of Lasso (L1 penalty) and Ridge (L2 penalty) regression, allowing for both variable selection and handling of feature correlation issues. We select the shrinkage coefficient λ through cross-validation, specifically by finding a λ value that minimizes prediction error and ensuring that this value is within one standard error range, which helps prevent overfitting and improves the model’s predictive ability on new data. We initially normalized the expression profiles of the samples in the TCGA-BRCA dataset using the Z-scale. Subsequently, we used the createDataPartition function in the caret package to split the samples into training and test sets with 80% and 20% allocations, respectively. In this study, we developed an elastic net penalized logistic regression model using only the training set. The shrinkage coefficient (lambda) was selected as a value within a standard error range to minimize the cross-validation prediction error rate, and the model feature with a minor prediction error was selected as the final marker gene. To generate a gene-based matrix risk score for the samples, Firth’s correction was used to calculate the odds ratios using the logistic function in the logistic package. The matrix of each sample risk score is the sum of the product of the risk ratio and the expression value of each marker gene; that is Matrix Risk Score  = ∑ i = 1 n z i β i , where n is the length of the marker gene, z i is the expression of gene i, and β i is the log-odds ratios of gene i. Subsequently, according to the dataset, the patient matrix risk score was used to determine the best grouping through the surv_cutpoint function and divide patients into high- and low-risk groups. The Kaplan–Meier test was used to compare differences in overall survival among the different sample groups.

We used a consensus clustering algorithm based on Monte Carlo references (Monte Carlo reference-based consensus clustering, M3C) (34) based on MRDEG expression to identify matrix-associated molecular subtypes. M3C is a consensus-clustering algorithm that involves using Monte Carlo simulations to mitigate the overestimation of K and effectively reject the null hypothesis of K=1. Real data were compared to eliminate bias, and statistical tests for the presence of structures were used to correct for inherent bias in consensus clustering. The optimal number of clusters K has the largest relative cluster stability index, the proportion of Monte Carlo P value is<0.05, and the fuzzy clustering score (Proportional Ambiguous Clustering, PAC) is the smallest. To confirm the accuracy of consensus clustering, the results were validated using a validation set. Subsequently, the R package ggpubr was used to generate a box plot, with the sample cluster labels as groups. Group differences were assessed for statistical significance using the Wilcoxon rank-sum test, with a p-value<0.05 indicating statistical significance.

Gene mutation analysis of stromal molecular subtype populations in different breast cancers

Breast cancer data were downloaded from GDC, and all non-synonymous mutations were selected for downstream analysis. R package map tools were used to display the related gene mutations, the biological functions affected by the mutations, and the classification of potentially druggable genes in different groups of breast cancer stromal molecular phenotype characteristics.

To identify the underlying molecular mechanisms of different stromal molecular subtypes in patients with breast cancer, we first performed ESTIMATE (35) on the TCGA-BRCA dataset. We analyzed and calculated four tumor-related scores, namely the matrix score, immune score, tumor purity, and ESTIMATE score, the immune score and matrix score calculated based on the ESTIMATE algorithm can facilitate the quantification of immune and matrix components in tumors; in this algorithm, immune and matrix scores are calculated by analyzing the specific gene expression characteristics of immune and matrix cells to predict the infiltration of non-tumor cells. Subsequently, the CIBERSORT algorithm (36) was applied to assess the infiltration status of immune cells within integrated datasets of various tumor samples. Next, differences in immune cell infiltration among different tumor subgroups were examined using the Wilcoxon test. Statistical significance was set at p<0.05. CIBERSORT ⁵ involves using linear support vector regression and serves as an R/web tool for deconvoluting expression matrices of human immune cell subtypes. It is used to evaluate the infiltration status of immune cells in sequenced samples using a gene expression signature set of 22 known immune cell subtypes. In addition, we analyzed the differential expression of immune checkpoint genes across different matrix molecular subtypes.

Seen in different subtypes of matrix molecules, we performed pathway enrichment analysis based on the 50 hallmark and C2 oncogenic pathways in patients with different subtypes. Pathway activity was assessed for each sample using the GSVA algorithm, and differentially active pathways were identified using a t-test.

We annotated the genes in the RNA-seq dataset as ligands and receptors using a curated database of human ligand-receptor pairs previously published by Ramilowski et al. (37). We retained only ligands corresponding to core matrix genes identified by Naba et al. (28) for subsequent analyses. The interaction score between a core matrix gene and its receptor was computed as the product of the expression values of the ligand (core matrix gene) and its cognate receptor in each sample. We identified the relative enrichment of ligand-receptor interaction scores among samples of different matrix subtypes using the Wilcoxon test and visualized the results using Circos.

All single-cell data analyses and integrations were performed using R software Seurat v 4.0.6. Two-cell quality control was implemented using the R Scrublet package. Cells with fewer than 300 genes, as revealed by single-cell sequencing, were deleted through quality control. Similarly, cells with more than 20% of the mitochondrial gene reads were deleted. The normalization and standardization of each sample data were realized through principal component analysis, and the inter-batch difference between samples was determined using the Harmony package. We used the t-distributed stochastic neighbor embedding algorithm to reduce dimensionality and visualize the single-cell data. The ECM scores of different cells were calculated using the AddModuleScore function.

For qPCR, total RNA was extracted using RNAiso Plus (TaKaRa, Japan), followed by reverse transcription using PrimeScript™ RT Master Mix (TaKaRa, Japan). qPCR was conducted using AceQ Universal SYBR qPCR Master Mix (Vazyme, China). The primer sequences are listed in Supplementary Table S2.

All data processing and statistical analyses were conducted using the R software ⁶ . To compare two groups of continuous variables, the independent Student’s t-test was used to assess the statistical significance of normally distributed variables, whereas the Mann–Whitney U test was used for non-normally distributed variables. The U-test (i.e., the Wilcoxon rank-sum test) was used to analyze the differences among non-normally distributed variables. The chi-square or Fisher’s exact test was used to compare and analyze the statistical significance of categorical variables between the two groups. The survival package in R was used for survival analysis, using Kaplan–Meier survival curves to illustrate the survival differences. The significance of the survival time difference between the two patient groups was evaluated using a log-rank test. Univariate and multivariate Cox analyses were performed using the survival package in the R software. All statistical p-values were two-sided, and p<0.05 is considered statistically significant.

Detailed R packages can be found in Supplementary Table S3 .

To further explore the underlying molecular mechanisms affecting the prognosis of patients with breast cancer, we conducted differential gene expression analysis on the complete TCGA-BRCA dataset to identify genes that were differentially expressed between patients in the long and short survival groups. Genes with an absolute value of log fold change (|logFC|) >0.25 and a p-value<0.05 were considered as DEGs associated with prognosis in patients with breast cancer, and 127 DEGs were identified ( Figure 2A ). Differential analysis revealed that 18.1% of the DEGs were stromal ( Figure 2B ). Of these genes, 3.9% were core matrix, and 14.2% were matrix-related. Subsequently, through Gene Ontology (GO) and KEGG analysis, biological processes and functions related to differentially expressed genes were identified. Among them, red represents biological processes, purple represents cellular components, blue represents molecular functions, and orange represents KEGG pathways. The p-values for all enriched functions are presented in the form of -log10(padj). In GO functional enrichment analysis, the analysis revealed the enrichment of biological processes associated with protein metabolism, including the negative regulation of endopeptidase activity, peptidase activity, and proteolysis. KEGG enrichment analysis indicated that the DEGs were associated with cytokines, including cytokine-cytokine receptor interaction and the chemokine signaling pathway ( Figure 2C ).

Recent advancements have underscored that individual matrix molecules rarely operate independently but as integral components within a dynamic three-dimensional supramolecular network comprising structurally and functionally integrated matrix constituents (41). We performed a correlation analysis of differentially expressed matrix genes in patients with breast cancer to determine whether these genes are also regulated in the disease ( Figure 2D ). Unsupervised clustering revealed two significant stromal body gene clusters related to somatic genes.

Volcano map display of differential analysis, in which red is the gene with up-regulated expression, and blue is the gene with down-regulated expression; (B) The proportion of differential matrix genes, red is the proportion of core matrix genes, and blue is the gene proportion of other matrix bodies; (C) Functional enrichment analysis of differentially expressed matrix genes; (D) Correlation heat map of differentially expressed matrix genes. GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

To explore the potential biological functions of different gene clusters, we performed GO and KEGG functional enrichment analyses on these two gene clusters. Gene Cluster 1 was primarily enriched in salivary secretion ( Figure 3A ), whereas gene Cluster 2 was primarily enriched in viral protein interactions with cytokines and cytokine receptors, chemokine signaling pathways, and cytokine-cytokine receptor interactions ( Figure 3B ). Concerning functional enrichment analysis with GO, gene Cluster 1 was found to be primarily enriched in the negative regulation of proteolysis and peptidase activity ( Figure 3A ). Gene Cluster 2 was primarily enriched in the chemokine-mediated signaling pathway, response to chemokines, and cellular response to chemokines ( Figure 3B ).

The enrichment result display of gene cluster1; (B) The enrichment result display of gene cluster2. GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

We identified 15 matrix-related marker genes using the elastic net penalized logistic regression method. To generate a matrix gene-based risk score for the samples, we used Firth’s correction to calculate the odds ratios( Figure 4A ) and the logistic function in the logistic package. Subsequently, patients were divided into high-and low-risk groups based on their matrix risk score. The Kaplan–Meier test was used to compare the differences in overall survival among the different sample groups. Survival analysis showed that the constructed stromal body risk signature could be used to accurately distinguish and predict patient prognosis ( Figures 4B, C ). Furthermore, significant differences were found between patients with different clinical characteristics. For example, in patients with breast cancer who died, it was significantly higher ( Figure 4D ); in older patients, it was also higher than that in younger patients ( Figure 4E ) and significantly lower in patients with T1 stage disease ( Figure 4F ).

To verify the effectiveness of our model, we applied our matrix body risk model to the GSE 20685 dataset, and the UCSC results of survival analysis on the Caldas 2007 dataset showed that our model could be used to significantly distinguish patients with breast cancer with different prognoses in the independent validation set ( Figures 5A, C ). The Receiver Operating Characteristic (ROC) curve analysis demonstrated that our model has a certain prognostic predictive ability and may have some clinical reference value; the Area Under Curve (AUC) of the GSE20685 dataset was 0.617 ( Figure 5B ), and that of the UCSC Caldas 2007 dataset was 0.571 ( Figure 5D ).

To explore which cell types express the marker genes we identified for constructing our stromal risk signature based on cell annotation information from a single-cell dataset (GSE 161529), we calculated positive ratios (positive score) and the enrichment degree of stromal risk genes with a negative ratio (negative score). CD4 Tconv, endothelial cells, epithelial cells, fibroblasts, malignant cells, Mono/Macro, pericytes, and plasma cells differed significantly between tumor and normal cells in the negative score ( Figure 6A ), and CD4 Tconv, endothelial cells, epithelial cells, fibroblasts, malignant cells, mono/macrophages, NK cells, pericytes, and plasma cells were significantly different between tumor and normal cells in the positive score ( Figure 6B ).

We used a consensus clustering algorithm based on a Monte Carlo reference (Monte Carlo reference-based consensus clustering, M3C) to identify the matrix-associated molecular subtypes based on the expression of matrix body-associated DEGs (MRDEGs). Five matrix-associated molecular subtypes were identified ( Figure 7A ). Survival analysis revealed that the five distinct stromal-associated molecular subtypes had significantly different survival rates, with patients in Cluster 3 having the worst prognosis ( Figure 7B ). To further explore the underlying molecular mechanism, Cluster 5 samples were significantly enriched in ECM interactions ( Figure 7C ), ECM proteoglycans ( Figure 7D ), and KEGG_ECM ( Figure 7E ) pathways, indicating that Cluster 5 samples had higher matrix body-related activity. To validate the feasibility of our clustering results, we performed the same clustering on the GSE20685 dataset and found that the samples clustered into four categories, with significant differences in survival. The lack of significance in cluster 4 may be due to the biological characteristics of the samples in this category, the small sample size, or the heterogeneity of clinical features. This indicates that the prognosis of patients in cluster 4 is relatively uniform, with no obvious survival differences. Additionally, the gene expression or related molecular pathways in this category may not have had a significant impact on patient prognosis, thus failing to achieve statistical significance in the survival analysis ( Figure 7F ).

Furthermore, we analyzed the differences in immunity between the different sample clusters. The ESTIMATE analysis indicated that Cluster 2 samples had the highest matrix ( Figure 8A ), immune ( Figure 8B ), and ESTIMATE ( Figure 8C ) scores and exhibited the lowest tumor purity ( Figure 8D ). Conversely, Cluster 1 showed the highest tumor purity. The CIBERSORT analysis revealed significant differences in the infiltration of various immune cells among patients with different molecular subtypes. Notably, CD8+ T cells and activated NK cells showed higher enrichment in Cluster 2 samples but lower enrichment in those of Cluster 1; however, T cells CD4 memory resting cells are more enriched in Cluster-5 samples and less enriched in Cluster-2 samples. Macrophages M0 cells are more enriched in Cluster-1 samples and less enriched in Cluster-5 samples. Macrophages M2 cells are more enriched in Cluster-3 samples and less enriched in Cluster-2 samples ( Figure 8E ).

Furthermore, significant differences were observed in the expression of immune checkpoint genes among patients with the five molecular subtypes ( Figures 9A, B ). These results indicate that stromal body genes may affect the prognosis of patients with breast cancer by regulating their immune response and infiltration.

The mutation characteristics of the above stromal-associated breast cancer subgroups were analyzed using the R package maftools. The Cluster 1 subtype primarily had mutations in TP53, TTN, and GATA 3 ( Figure 10A ); the Cluster 2 subtype primarily had TP53, TTN, and PIK3CA mutations ( Figure 10B ); the Cluster 3 subtype primarily had PIK3CA, TP53, and KMT2C gene mutations ( Figure 10C ); the Cluster 4 subtype primarily exhibited PIK3CA, GATA3, and TP53 gene mutations ( Figure 10D ); the Cluster 5 subtype primarily had PIK3CA, CDH1, and TP53 gene mutations ( Figure 10E ).

Furthermore, we analyzed the mutations of the three patient subtypes to explore the gene druggability and the interaction between drugs and genes (from Drug Gene Interaction database, DGIdb database) and found that the genes to predict that the drug might act on Cluster 1, 2, 3, 4, and 5 subgroups are DRUGGABLE GENOME (FCGBP, HMCN1, MUC16, MUC17, and OBSCN) ( Figure 11A ), DRUGGABLE GENOME (CDH1, DST, FAT3, HMCN1, and MUC16) ( Figure 11B ), DRUGGABLE GENOME (CDH1, HMCN1, MAP2K4, MAP3K1, and MUC16) ( Figure 11C ), CLINICALLY ACTIONABLE (ARID1A, CBFB, CDH1, GATA3, and KMT2C) ( Figure 11D ), and DRUGGABLE GENOME (ABCA13, CDH1, HMCN1, MAP3K1, and MUC16) ( Figure 11E ), respectively, indicating that these mutated genes can be used for subsequent studies on the development of drug targets.

Subsequently, we calculated the scores of marker genes related to the matrix (CCL5, CLEC3A, CST1, CST7, CTSW, CXCL9, EMID1, HPX, IGFALS, LAMB3, PLAT, S100B, SCUBE2, SEMA3B, and SERPINA1) using the ssGSEA algorithm. Comparing the Figure 12A score of different typing scores of samples with the sample typing information revealed that the Cluster 1 score was the lowest; we defined it as the ECM-low group. Cluster 5 score was the highest, and we defined it as the ECM-high group. The hallmark ( Figure 12B ), and C2 ( Figure 12C ) enrichment pathway analyses for different patient groups, the colors in the heatmap indicate the relative expression levels: red for high expression, blue for low expression. The clustering on left shows the hierarchical relationship of samples based on their gene expression profiles. It revealed that the ECM-high group samples were primarily enriched in APOPTOSIS, HALLMARK IL2 STAT5 SIGNALING, HALLMARK TNFA SIGNALING VIA NFKB, HALLMARK KRAS SIGNALING UP, HALLMARK EPITHELIAL MESENCHYMAL TRANSITIO, HALLMARK COAGULATION, HALLMARK INTERFERON ALPHA RESPONSE, HALLMARK INFLAMMATORY RESPONSE, HALLMARK ESTROGEN RESPONSE EARLY, HALLMARK COMPLEMENT, HALLMARK INTERFERON GAMMA RESPONSE, and HALLMARK ALLOGRAFT REJECTION, and the ECM-low group samples are primarily enriched in LI_CISPLATIN_RESISTANCE_DN, LI_CISPLATIN_RESISTANCE_UP, KANG_CISPLATIN_RESISTANCE_DN, and BRACHAT_RESPONSE_TO_CISPLATIN. Apoptosis, IL2 stat5 signaling, Tnfa signaling via NFKB, Kras signaling up, Epithelial-mesenchymal transition, Coagulation, Interferon alpha response, Inflammatory response, Estrogen response early, Complement, Interferon-gamma response, Allograft rejection, and the ECM-low group samples are mainly enriched in Cisplatin resistance dn, Cisplatin resistance up, Kang cisplatin resistance dn, and Brachat response to cisplatin.

Similarly, we assessed the expression of matrix-related marker genes (CCL5, CLEC3A, CST1, CST7, CTSW, CXCL9, EMID1, HPX, IGFALS, LAMB3, PLAT, S100B, SCUBE2, SEMA3B, and SERPINA1) in various breast cancer subtypes ( Figure 13 ). The results showed that different markers, such as CTSW and S100B, were specifically overexpressed in Luminal A breast cancer cells. CST1, EMI D1, and HPX SCUBE2 were specifically overexpressed in Luminal B breast cancer cells. CCL5, CLEC3A, CTSW, CXCL9, IGFALS, and SEMA3B were specifically overexpressed in HER2-positive breast cancer cells. CST7, PLAT, and SERPINA1 were specifically overexpressed in Basal-like breast cancer cells.

We used ESTIMATE to assess tumor purity between different groups (ECM-high vs. ECM-low), revealing a notable difference in the tumor purity between them (p<0.05, Figure 14A ). Furthermore, the ECM-high group exhibited lower tumor purity owing to its higher matrix and immune scores (p<0.05; Figure 14B ). Moreover, we developed a reference matrix for CIBERSORTx using the cell types identified in the single-cell dataset (GSE161529). Deconvolution methods were used to calculate the scores of TCGA-BRCA samples for different cell types.

The findings revealed that samples from the ECM-low group exhibited higher malignant tumor scores ( Figure 14C ). Conversely, samples from the ECM-high group had a significant enrichment of immune cells, particularly B cells, macrophages, monocytes, and CD4 T cells ( Figure 14D ). Furthermore, the ECM-high group samples demonstrated significant enrichment of stromal cells, specifically endothelial cells and fibroblasts, whereas the ECM-low group samples were notably enriched in epithelial cells ( Figure 14E ).

ECM components interact directly with cell surface receptors, regulating the activity of numerous signaling pathways, including those related to epithelial-mesenchymal transition (EMT) and ECM production. We conducted a ligand-receptor interaction analysis to elucidate the potential direct effects of these maturation forms on cell signaling. The results showed that the matrix body genes CCL19, CCL5, CXCL9, LTB, MMP12, PLAT, SEMA3B, and SERPINA1 interacted with many receptors ( Figure 15A ). Interacting receptors are crucial in cancer development, participating in the IL-6/JAK/STAT 3 signaling pathway ( Figure 15B ).