We conducted a retrospective collection of 14 distinct breast cancer cohorts sourced from The Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), and Metabric (17). These cohorts comprised samples with comprehensive survival data, which were subsequently utilized for in-depth analysis. In total, our study encompassed 22,162 patients across the 14 cohorts for the purpose of prognostic assessment. The distribution of patients across cohorts was as follows: TCGA-BRCA (n = 1076), GSE202203 (n = 3206), GSE96058 (n = 3409), GSE20685 (n = 327), GSE86166 (n = 330), GSE1456 (n = 159), GSE21653 (n = 244), GSE7390 (n = 198), GSE11121 (n = 200), GSE6532 (n = 87), GSE88770 (n = 108), GSE48391 (n = 81), GSE131769 (n = 298) and Metabric (n = 1747).

To develop a LLPS signature specific to breast cancer, we adopted the methodology outlined by our previous research (18). Our approach entailed the utilization of ten diverse computational techniques: Random Survival Forest (RSF), Least Absolute Shrinkage and Selection Operator (LASSO), Gradient Boosting Machine (GBM), Survival Support Vector Machine (Survival-SVM), Supervised Principal Component (SuperPC), Ridge Regression, Partial Least Squares Cox Regression (plsRcox), CoxBoost, Stepwise Cox regression, and Elastic Net (Enet). Notably, RSF, LASSO, CoxBoost, and Stepwise Cox were selected for their capacity to diminish dimensionality and identify pertinent variables.

We then combined the machine learning algorithms into 108 unique combinations to construct a robust MDLS model. Each combination was trained on multi-omics data, with the average Concordance Index (C-index) used as the performance metric to identify the most predictive model. Through iterative cross-validation, we systematically evaluated the predictive accuracy of each algorithm combination across multiple breast cancer cohorts.

From the analysis, we identified the most consistent and predictive model. Four LLPS-related genes (POP1, TUBA1C, RACGAP1, and PLK1) were selected as key features based on their prognostic value, as determined by univariate Cox regression analysis. These genes formed the foundation of the final MDLS signature, which was optimized to predict patient outcomes in breast cancer.

To define the high- and low-MDLS groups, we used the surv_cutpoint function from the “survminer” R package. This function calculated the optimal cutoff value that maximally separated patients into high- and low-risk groups based on their survival data. Patients with MDLS risk scores above the cutoff were classified as the high-MDLS group, indicating a higher risk profile, while those with risk scores below the cutoff were classified into the low-MDLS group.

The performance of the MDLS signature was validated across 14 independent breast cancer cohorts, incorporating both bulk tumor and single-cell RNA sequencing data. In total, the cohorts involved over 9,723 breast cancer patients, ensuring a comprehensive evaluation. Additionally, we compared the MDLS signature with 69 published breast cancer signatures, and the MDLS demonstrated superior prognostic power across the cohorts.

To delineate genetic disparities between the two MDLS groups, we conducted an analysis of both genetic mutation levels and Copy Number Alterations (CNA) utilizing the TCGA-BRCA database.

The Tumor Mutation Burden (TMB) of both high- and low MDLS breast cancer patients was extracted from the raw mutation file. Utilizing the maftools landscape, we depicted the most frequently mutated genes (mutation rate > 5%). Additionally, patient-specific mutational signatures were identified using the deconstructSigs package (19). Notably, we emphasized four prominent mutational signatures (SBS3, SBS1, SB12, SBS11) within the TCGA-BRCA dataset that displayed heightened mutation frequencies. We identified the five most common regions of amplification and deletion and specifically highlighted the four predominant genes in chromosomal regions 8q24.21 and 5q11.2.

For the preparation of the dataset for single-cell RNA sequencing analysis, we employed Seurat (v4.0) to process the data extracted from GSE161529 (20). This process involved filtering out genes with no expression and retaining those with nonzero expression levels. Normalization of the expression matrix was performed using Seurat’s “SCTransform” function. Dimensionality reduction of the dataset was achieved through (principal component analysis) PCA and UMAP reductions. To identify distinct cellular groupings, Seurat’s “FindNeighbors” and “FindClusters” functions were utilized. To ensure dataset integrity and reliability, the DoubletFinder package was employed to eliminate potential doublets (21). Cells failing to meet defined quality standards, such as exhibiting mitochondrial gene content exceeding 15% or containing fewer than 500 genes, were excluded. Through stringent quality control measures, a total of 47,784 cells were retained for subsequent analysis. Cell types were determined by manual annotation based on the presence of established marker genes.

In our investigation, we adopted the Single-Cell rEgulatory Network Inference (SCENIC) approach to construct gene regulatory networks (GRNs) from single-cell RNA sequencing data. SCENIC involves a three-step process: initially, it identifies co-expression modules between transcription factors (TFs) and their potential target genes. Subsequently, for each module, it identifies the direct target genes, prioritizing those enriched with the motif of the associated TF. A regulon is then defined, comprising a TF and its direct targets. Finally, the regulatory activity score (RAS) is computed for each cell by assessing the area under the recovery curve.

While the conventional SCENIC protocol encounters challenges with scalability for extensive datasets and is susceptible to variations in sequencing depth, we introduced a modification to enhance both scalability and robustness. This involved partitioning the data into metacells before applying SCENIC to these gene expression profiles (22). This adjustment significantly enhances data quality and reduces computational demands, representing a notable advancement in the application of SCENIC to single-cell RNA-seq data analysis.

We implement a comprehensive computational approach to elucidate the regulatory interplay between transcription factors (TFs) and their corresponding target genes, with a particular focus on TF clustering. Initially, the method entails filtering TF-target interaction data to isolate pairs that exceed a predefined significance threshold (>1), ensuring prioritization of regulatory interactions of utmost relevance. Subsequent analysis aims to identify key regulatory TFs by assessing the extent of their target gene regulation, highlighting them as central nodes within the regulatory network for in-depth investigation.

To visually represent the complex network of TF-target interactions, an undirected graph model is constructed. The spatial arrangement of this graph is refined using a force-directed algorithm to intuitively depict the network’s architecture, emphasizing the interplay between TFs and their targets. Additionally, to enhance comprehension of the network’s structure, the Leiden algorithm is applied for community detection. This process unveils the modular configuration of TFs based on their regulatory interconnections, assigning each TF to a distinct cluster. This facilitates a nuanced analysis of the regulatory landscape, enabling insights into the functional organization of TFs within the network.

Using the “CellChat” R package, CellChat objects were generated based on the UMI count matrices for each respective group (23). The “CellChatDB.human” database served as the reference for ligand-receptor interactions. Interpreting intercellular communication was executed using the default settings provided by the package. To assess and compare interaction counts and intensities, CellChat objects from each group were merged using the “mergeCellChat” function. Variations in the number and intensity of interactions among specific cell types across different groups were visualized using the “netVisual_diffInteraction” function. Additionally, changes in signaling pathways were identified with the “rankNet” function, and the distribution of signaling gene expression among the groups was illustrated using both the “netVisual_bubble” and “netVisual_aggregate” functions.

Furthermore, we applied the NicheNet package to analyze intercellular communication from the perspective of ligand activity and the expression patterns of specific downstream targets regulated by these key ligands (24). This approach enables a detailed understanding of the signaling processes underlying interactions between different cell types, leveraging information about ligand-target relationships to infer communication pathways within the cellular microenvironment.

In our endeavor to comprehensively and accurately assess immune cell infiltration levels, we conducted an analysis of adverse infiltrated immune cells across multiple algorithms, including MCPcounter, EPIC, xCell, CIBERSORT, quanTIseq, and TIMER, among patients stratified by the MDLS (25–31). Additionally, to depict the immune landscape and architecture within the tumor microenvironment (TME) with precision, we evaluated the ESTIMATE and TIDE indices. These measures provide crucial insights into the potential for immunotherapy and offer prognostic implications for breast cancer patients.

Furthermore, we quantified immune checkpoints, which serve as indicators of the immune state and offer preliminary predictions of patient responsiveness to immune checkpoint inhibitor (ICI) therapy. This comprehensive approach to evaluating the immune profile within the TME is essential for advancing personalized medicine and refining treatment strategies for breast cancer patients.

Our methodology for identifying therapeutic targets and drugs for high MDLS patients commenced by filtering out duplicate compounds from the Drug Repurposing Hub, resulting in a refined list of 6,125 compounds. The selection of therapeutic targets associated with breast cancer outcomes was established through Spearman correlation analysis. Specifically, we assessed the relationship between the MDLS and gene expression levels, selecting genes with a correlation coefficient greater than 0.3 and a P-value less than 0.05. Additionally, genes demonstrating a correlation coefficient below -0.3 and a P-value below 0.05 were identified as linked to poor prognosis. The significance of these genes was further evaluated by examining the relationship between CERES scores from the Cancer Cell Line Encyclopedia (CCLE) and risk scores (32).

To enhance predictions regarding drug responsiveness, we utilized data from the Cancer Therapeutics Response Portal (CTRP) and the PRISM project, both of which offer extensive drug screening and molecular data across diverse cancer cell lines. Differential expression analysis was conducted between bulk samples and cell lines. Subsequently, the pRRophetic package was employed to implement a ridge regression model for predicting drug response. This model, trained using expression data and drug response metrics from solid Cancer Cell Lines (CCLs), demonstrated excellent predictive accuracy, validated through 10-fold cross-validation (33).

Furthermore, to identify the most promising therapeutic drugs for breast cancer, Connectivity Map (CMap) analysis was performed. This entailed comparing gene expression profiles across different risk subgroups and submitting the top 300 genes (comprising 150 up-regulated and 150 down-regulated genes) to the CMap website. Interestingly, a negative CMap score indicated a higher therapeutic potential against breast cancer, suggesting an inverse relationship between the CMap score and a compound’s effectiveness as a potential treatment.

To assess gene expression in breast cancer specimens, RNA extraction was performed using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), followed by cDNA synthesis and quantitative reverse transcription PCR (qRT-PCR) using GoScript reverse transcriptase and Master Mix (Promega), following the manufacturer’s guidelines. The CFX96 Touch Real-Time PCR Detection System (BioRad, Hercules, CA, USA) was employed for data acquisition. Gene expression quantification was executed using the 2^-ΔΔCq method, with GAPDH utilized as the normalization control. Subsequently, patients were categorized based on their gene expression profiles, utilizing a predefined formula derived from the MDLS. This stratification played a crucial role in identifying patients with distinct risk profiles, thereby enabling tailored therapeutic interventions.

Tissue samples were obtained from 30 breast cancer patients undergoing surgery at Guizhou Provincial People’s Hospital ( Supplementary Table S1 ). These samples underwent Hematoxylin and Eosin (H&E) staining using established protocols. The diagnosis was independently confirmed by two pathologists.

For the immunohistochemistry (IHC) analysis, procedures for paraffin-embedded samples were followed, as outlined in previous studies (34, 35). Protein expression levels were assessed independently by two pathologists, adhering to standardized protocols and scoring systems consistent with methodologies from prior research (35).

Statistical analyses were performed to assess the differences between the high- and low-MDLS groups. Continuous variables were compared using the unpaired Student’s t-test or Mann-Whitney U test, depending on the distribution of the data. Chi-square test or Fisher’s exact test was employed to compare categorical variables. Pearson’s correlation test was used to evaluate correlations between gene expression levels and clinical variables.

Survival outcomes were analyzed using Kaplan-Meier analysis, and differences between groups were evaluated with the log-rank test. Cox proportional hazards regression models were used for both univariate and multivariate survival analyses to identify independent prognostic factors.

All statistical analyses were conducted using R (version 4.0.5), and the survminer R package. Statistical significance was defined as a p-value of less than 0.05, with differences considered significant at *p < 0.05, **p < 0.01, or ***p < 0.001.

From the DrLLPS database (36), we collected LLPS genes to conduct differential expression analyses in TCGA-BRCA between tumor and normal tissues. A ten-fold cross-validation process was applied to construct prediction models using 108 algorithm combinations and calculated the mean C-index of each algorithm in the training cohort (TCGA-BRCA) and 8 external cohorts. The combination of StepCox[forward] and survival-SVM with the highest mean C-index (0.649) was selected as the final model ( Figure 1A ). We assessed the prognostic value of these LLPS genes based on univariate Cox regression and calculated the Hazard ratio (HR) for these genes in the nine enrolled cohort ( Figure 1B ). As shown in Figure 1C , four positively correlated risk factor genes (POP1, TUBA1C, RACGAP1, PLK1) were screened to build the machine learning derived LLPS signature (MDLS). The findings revealed that our model successfully differentiated between high- and low-MDLS patients, suggesting that MDLS can offer valuable reference information for predicting the survival of breast cancer patients ( Supplementary Figure S1 ).

Univariate and multivariate Cox analysis showed that MDLS was an independent risk factor compared with other clinical indicators ( Supplementary Figure S2A ). The nomogram consisting of MDLS, stage and age was resorted to accurately predict the OS of BC patients at distinctive times, and the prediction effect of MDLS is better ( Supplementary Figures S2B–F ). The kernel-smoothing hazard plot show that high-MDLS patients had a poorer outcome and higher recurrence frequency than low-MDLS patients ( Supplementary Figure S2G ). To evaluate the stability of the predictive model of the MDLS, 68 published signatures in BC were manually collected and assessed in 10 independent cohorts. We demonstrated that only the MDLS had consistent statistical significance across all cohorts ( Figure 2A ). We compared the predictive power of MDLS with those 69 features across 10 cohorts using the C-index ( Figure 2B ). Our model showed significantly better accuracy than the others in almost all cohorts (ranking first in six cohorts, second in one cohort, fourth in one cohort, fifth in one cohort, and ninth in one cohort), revealing the stability of the MDLS.

To explore genomic heterogeneity between high- and low-MDLS groups, we analyzed gene mutation and copy number variation across the groups ( Figure 3A ). We observed high tumor mutation burden (TMB) in the high-MDLS patients ( Figures 3A, C ). Combining the 10 oncogenic signaling pathways in the TCGA database, we found that a total of 12 classical tumor suppressor genes and 6 oncogenes were mutated more frequently in the high-MDLS group ( Figure 3B ). Moreover, mutational signature analysis shows that the frequencies of SBS2, SBS13, SBS7b, SBS7d were significantly higher in low-MDLS ( Figures 3A, D ). The number of classical mutations (SBS2 and SBS13) in APOBEC3 was examined, the results showed the same trend ( Figure 3D ). Next, we delved deeper into the CNA scenery of the two groups. Compared to the low-MDLS group, the high-MDLS group owned evidently higher deletion or amplification in the chromosome arm levels, like the amplification of 3q26.32, 4q13.3, 8q24.21, 10p15.1, 12p13.33, and the deletion of 5q11.2, 5q21.3, 14q24.1, 15q13.1, 19q13.32 ( Figures 3A, E ). At the gene level, there was significant gene amplification on chromosome 8q24.21 (PVT1, MYC, CCDC26, GSDMC) and significant gene deletion on chromosome 5q11.2 (GPBP1, RAB3C, DDX4, ITGA1) in the high-MDLS group ( Figure 3A ). In conclusion, tumor suppressor gene deletion and oncogene amplification in the high-MDLS group may be one of the reasons for the poor prognosis.

We selected 8 patients with breast cancer for further evaluation of MDLS, including 4 tumor tissues and 4 normal tissues ( Supplementary Figures S3A, B ). The cells were divided into 17 clusters and 7 cell types ( Figures 4A, B ). Statistical analysis was conducted to determine the overall number and proportion of these seven types of cells in the bodies of these eight tumor patients ( Supplementary Figures S3C, D ). Representative markers for each cell type are shown ( Figure 4C ; Supplementary Figure S3E ). Single-cell sequencing revealed transcriptome differences for each cell type between tumors and normal tissues ( Figure 4D ). The results show that macrophages, T cells, and epithelial cells are notably enriched in tumor tissues, while other cells are highly represented in normal tissues.

Next, MDLS was integrated into the single-cell analysis for scoring ( Figure 4E ). The cells were segregated into two groups based on the media MDLS scores of the epithelial cells ( Figure 4F ). The potential pathways of MDLS were enriched and visualized by differential expression analysis and GSEA ( Supplementary Figures S3F, G ). A tumor microenvironment consists of tumor epithelial cells, tumor stromal cells, and extracellular matrix. Most cancers are caused by a destroyed epithelial cell population, causing tumor cells to grow rapidly. Take the epithelial cells for example, high-MDLS group was notably enriched in proteasome, focal adhesion, Ribosome, spliceosome. While the low-MDLS group was predominantly associated with reactive oxygen, oxidative phosphorylation ( Supplementary Figure S3G ). Further observation on copy-number alteration by the copyKat algorithm was employed to distinguish between normal cells and tumor cells ( Figure 4G ). Ultimately, a higher MDLS score was observed in tumor-aneuploid cells compared to tumor-diploid cells, indicating the significance of MDLS in breast cancer progression (Figure 4H).

We used the SCENIC pipeline to analyze single-cell RNA-seq data with cis-regulatory sequence information to comprehensively construct GRNs. The UMAP analysis visualized the cell data into seven cell clusters and performed RAS scores on the MDLS score ( Figures 5A, B ). In summary, the gene expression information was converted to the activity score of transcription factors known as RAS. PCA and variance analyses revealed the transcription factors associated with the formation of MDLS and cell types ( Figures 5C, D ).

We identified seven key regulatory factors for cell types and scored each regulator specificity according to Jensen-Shannon divergence. The regulator with higher regulon specificity score (RSS) might be correlated specifically with that cell type ( Figure 5E ). We selected the regulator with the highest RSS value for each cell type to further examine its functional properties ( Figure 5F ; Supplementary Figure S4A ). Take the epithelial cells for example, PBX1, SPDEF and ATF3 were found to be the most specific regulatory factors in the RSS sequencing of epithelial cells ( Figure 5F ).

Each cell type has its own form and function, and the characteristics of the cell type need to be maintained by the coordinated interaction of transcription factors and their corresponding target genes. According to the Leiden algorithm, we compared the RAS scores of each regulon pair across the atlas to characterize the combinatorial patterns of the MDLS. A highly modular diagram shows the formation of 12 modules ( Figure 5G ; Supplementary Figure S4B ). The module B and C play a key role in MDLS progression ( Figure 5H ; Supplementary Figure S4B ). We focus on the key transcription factors that drive epithelial cell transcriptome changes in MDLS. Using GSEA analysis to identify multiple pathway variants in epithelial cells, MAPK/MTORC1 signaling and G2M checkpoint were activated in the high-MDLS group ( Figures 5I, J ). The results further confirmed that transcription factors are involved in these pathways and the progression of MDLS ( Figure 5K ). The detailed regulatory network of transcription factors that influence MDLS progression were further demonstrated ( Figure 5L ).

To emphasize the complex interactions between cells interactions in BC progression, we applied the CellChat analysis to assess communication between seven different types based on MDLS. The number and strength of interactions were assessed, revealing that the group with high-MDLS had higher cell-to-cell communication ( Figure 6A ). In the high-MDLS group, endothelial cells, epithelial cells, and fibroblasts displayed a massive amount of interaction. Nevertheless, macrophages and T cells have weaker interactions with other cells ( Figure 6B ). We further investigated 51 signaling pathways of cell communications, therein, 48 signaling pathways were primarily activated in high-MDLS cells, including laminin, collagen, and cxcl, except for FN1, MK, and CCL ( Figure 6C ). The interaction between endothelial and epithelial cells in the high-MDLS group was stronger, while the communication with T cells was weaker ( Figure 6D ). Using epithelial cells as an example, the high-MDLS group had higher pathway specificity compared to the low-MDLS group, such as MIF, VISFATIN, LAMININ and THBS pathway ( Figure 6E ).

Nichenetr analysis was utilized to assess the activity of ligands regulating epithelial cell incoming and outgoing in different MDLS groups. Further analyses were focused on the differences in the activity of the ligand-receptor pairs. A Circos diagram illustrates the interaction of ligand-receptor pairs in different cell types ( Figure 6F ). The high interaction of TGM2-SDC4 indicates that fibroblasts and endothelial cells are the primary transmitter cells influencing changes in the epithelial cell pathway ( Figure 6G ). Figure 6H provides a detailed roadmap of TGM2 ligand reaching the target receptor SDC4 through other receptors or transcription factors.

We applied six algorithms to evaluate the immune infiltration. High-MDLS group had more immune infiltration, such as CD4⁺ memory T cells, CD8⁺ T cells and M1 macrophages ( Figure 7A ). The expression of ICIs is also assessed as a key indicator of immunotherapy responsiveness. The expression of ICIs was higher in the high-MDLS group, such as TNFRSF14, PD-1, PD-L1, LAG3, IDO1 ( Figure 7B ). IHC was performed to support the above results using the representative cell markers and clinical ICIs ( Figure 7C ).

Next, we used ESTIMATE analysis to predict tumor purity and the presence of infiltrating stroma/immune cells in tumor tissue. The results showed that low-MDLS patients had higher ESTIMATE scores and Stromal scores compared to the higher score patients but had lower tumor purity ( Figure 8A ). It suggests that the low-MDLS patients are more likely to receive immunotherapy. Meanwhile, TIDE value, dysfunction and exclusion value is higher in low-MDLS patients ( Figure 8B ). Notably, patients with low-MDLS combined with high-TIDE had a higher survival rate than patients with other types ( Figure 8C ). The results showed that the anti-tumor immune activity of low-MDLS patients was higher than that of high-MDLS patients ( Figure 8D ). Immune checkpoints have long been used in the immunotherapy of cancer, so we evaluated the ability of the MDLS to predict the immune checkpoint blocking response. MDLS in the IMvigor210 (anti-PD-L1) and GSE78220 (anti-PD-1) cohorts was further assessed. In IMvigor210 ( Figures 8E–H ) and GSE78220 ( Figures 8I–L ), patients with low-MDLS exhibited better survival rates and clinical benefits compared to those with high-MDLS. In summary, patients with low-MDLS may derive greater benefits from ICIs treatment.

Chemotherapy remains a cornerstone treatment for cancer. In this study, we devised a targeted approach for breast cancer patients with high-MDLS levels, leveraging sensitivity data collected from multiple datasets. Initially, we employed Spearman’s correlation analysis to identify key therapeutic targets. The analysis revealed a positive correlation between MDLS and the abundance of six potential targets (PSMA7, PRMT5, SLC25A13, INCENP, TREM1, FOXM1). Importantly, these targets also exhibited notably negative correlations with their respective CERES scores, indicating their potential as therapeutic targets for patients with high-MDLS levels ( Figure 9A ). Therein, five of six targets were found to be closely linked to various drug action pathways, underscoring their significance as critical therapeutic targets for this patient subgroup ( Figure 9B ).

Subsequently, we obtained 12 chemical compounds from the CTPR and PRISM datasets, including. An analysis comparing the AUC values of these compounds between the two MDLS groups revealed higher AUC values in patients with low-MDLS, indicating a less favorable response to chemotherapy in this demographic ( Figures 9C, D ). A comprehensive multiple-perspective analysis was then conducted to select the most effective therapeutic drugs from these 12 candidates. This analysis included detailed evaluations of each compound’s clinical status, experimental evidence, mRNA expression levels, and CMap scores. Ultimately, methotrexate was identified as the most suitable therapeutic drugs for patients with high-MDLS, based on their CMap score ( Figure 9E ).