We used the UCSC Xena Browser (https://xenabrowser.net/datapages/) (17) to download the TCGA-BRCA dataset in the transcriptome fragments per kilobase million (FPKM) format; this dataset contains survival data and corresponding clinicopathological data for 1280 patients. The single nucleotide variant (SNV) data, copy number variation (CNV) data, methylation data, mutation count, MSI sensor data and fragment genomic change frequency information of the corresponding samples were downloaded from the cBioPortal (18) database.

In addition, we downloaded the GSE88770 (19) and GSE20711 (20) transcriptome chip data and the clinical information of BRCA patients from the GEO (21) database. The GSE88770 dataset was generated by sequencing via the [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. A total of 117 human BRCA tumor samples were selected as the validation set, and quality control steps were applied to exclude samples with low quality (based on background noise and signal intensity). The GSE20711 dataset was also generated by sequencing via the [HG-U133_Plus_2] [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, which contains information on 88 human BRCA tumor samples. Normal tissue samples from two patients were included as controls to analyze gene expression differences.

In addition, we obtained the datasets from GSM5457205 (22), which included information on human BRCA tumor samples. The Illumina NovaSeq 6000 (Homo sapiens) platform was used to generate the sequencing data, and initial data processing included quality control steps, such as filtering out low-quality reads and performing adapter trimming. We further obtained data from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb) (23). The online database includes 742 genes associated with lipid metabolism, and functional enrichment analysis will be performed to explore the biological significance of these genes in the context of BRCA.

We used single-cell count data from the original UMI. In addition, Seurat v4.0 was used for preprocessing steps, quality control, normalization, dimensionality reduction clustering and clustering. Specific quality control criteria were established to ensure data integrity. First, every gene had to be expressed in a minimum of three cells, and each cell was required to expressed at least 200 genes. Second, the genes were selected based on the number of expressed genes in each sample, using the median ± 3*MAD (median absolute deviation) standard for filtering. Moreover, thresholds of 10% for the proportion of mitochondrial gene proportion and 1% for the proportion of hemoglobin gene proportion were set according to different sample types to exclude potential low-quality or dying cells. After applying quality control measures, we processed with data analysis using the Seurat package’s default parameters and standard operating procedures, which included: standardization of counts, normalization of data, identification of highly variable genes, dimensionality reduction (using PCA), and clustering of cells. The Harmony package was utilized to integrate data from multiple samples, thus correcting for batch effects. The annotation of single-cell groups was performed both manually and with the assistance of the scType package (23). Finally, the FindAllMarkers function used the Wilcoxon signed-rank test to compute differential gene expression between clusters, thereby identifying marker genes for each cell type.

On the basis of the lipid metabolism genes in the MSigDB database, we performed gene scoring on the BRCA tumor tissue samples in the TCGA-BRCA dataset, and used the gene set variation analysis (GSVA) method to perform unsupervised scoring on each sample by gene set. The samples were subsequently divided into two groups according to the score ranking according to the GSVA score: the samples with the top 20% score were the Score+ group, and the samples with the bottom 20% score were the Score- group. We analyzed the DEGs in the TCGA-BRCA dataset between lipid metabolism score groups in human BRCA tumor samples using the rank-sum test.

We used univariate Cox regression based on TCGA-BRCA data to screen genes associated with survival from DEGs grouped by lipid metabolism score. We integrated nine machine learning algorithms into 184 machine learning algorithm combinations, including LASSO, Ridge, elastic network (Enet), StepCox, survival support vector machine (survivalSVM), CoxBoost, supervised principal components (SuperPC), random survival forest (RSF), and generalized boosted regression modeling (GBM), to develop prognostic models.

The specific algorithm and parameters are as follows: (a) LASSO achieves feature selection through L1 regularization, with a regularization intensity λ set to 0.1, demonstrating excellent feature selection ability while effectively balancing model complexity and data fitting. (b) Ridge utilizes L2 regularization with a λ value of 0.1, which enhances robustness on the training set, ensuring model stability and adaptability to high-dimensional data while mitigating potential overfitting issues. (c) Enet integrates both L1 and L2 regularization techniques, making it suitable for scenarios where the number of features exceeds the number of samples. The α value is determined as 0.4 based on cross-validation results, enabling effective handling of multicollinearity problems by properly combining the advantages of LASSO and Ridge methods. (d) StepCox regression adopts the “direction=forward” method for feature selection based on the Cox proportional risk model. This progressive variable addition approach ensures focus on influential key genes in constructing a final model that can effectively identify important variables. (e) The survivalSVM possesses capabilities in dealing with survival time and event states using default parameter settings for model training purposes. With extensive application experience in survival analysis field, this algorithm generally captures complex patterns within survival data quite well. (f) CoxBoost employs Boosting concept to iteratively enhance the Cox model by utilizing a combination of 10 trees aiming at optimizing model fitting effect without excessive complexity introduction. (g) The SuperPC utilized default parameter settings to generate a linear combination of relevant features that captures the direction of greatest variation in the dataset. Optimal threshold evaluation was performed using cross-validation results, and application of the “pre-validation” function prevented issues with fitting multivariate Cox regression models to validation datasets. (h) RSF classification and regression, combining multiple decision trees with n_estimators set at 100 and default settings is widely accepted as it provides sufficient decision trees for increased prediction accuracy while ensuring model stability. (i) GBM improves model generalization by combining multiple weak learning models (usually decision trees) with a learning rate between 0.1-0.3 and tree depth set at 5 based on step-by-step experimentation and observation to ensure gradual growth without compromising stability.

After the construction and tuning of these models, we used two datasets, GSE88770 and GSE20711, as validation sets, and the average C-index as the main evaluation index of the model. The combination algorithm that achieved the highest C-index was chosen as the ultimate model. In addition, we collected literature on prognostic model construction from the TCGA-BRCA dataset from September 1, 2022 to September 1, 2023 (24–28).During this time, we conducted a literature survey to ensure that all the latest research findings were covered and corresponding data were used to compare the performance of the models.

First, we used the rank-sum test to compare clinical features between patients with high and low lipid metabolism scores, and the results are shown as bars. We then obtained data from the MSigDB (23) (https://www.gsea-msigdb.org/gsea/msigdb) database to download 50 tumor-associated pathways for gene set enrichment analysis (GSEA) (24, 29)to investigate the distinctions between groups with high and low scores in lipid metabolism. The transcriptome data and clinical characteristics of 32 tumors were downloaded from UCSC Xena, and the best machine learning combination model was subsequently used to analyze the prognostic effect of the model score in different tumors one by one.

We examined the variations in the tumor mutation burden (TMB), SNV data, and CNV data. We then selected the most important genes that make up the lipid metabolism score model for analysis, and we analyzed the differences in the transcript levels and methylation levels of these genes and the effects of individual genes on survival. For the selection criteria of the key genes in the lipid metabolism scoring model, first, we used univariate and multivariate Cox regression analyses to screen genes, and the selection criterion was a P value less than 0.05. Second, among the genes screened, those with a relatively high risk ratio (HR) and biological significance were prioritized to ensure that these genes not only were statistically significant, but also played an important role in lipid metabolism processes. We also used a combination of multiple machine learning algorithms (such as LASSO, Ridge, and CoxBoost) to construct prognostic models, and evaluated the predictive efficacy of each model. Finally, we selected the combination algorithm that performed best in the training set and validation set to screen 21 genes as the gene set of the LMPGS model.

We analyzed information on 26 immunogenic cell death-related genes and 47 immune checkpoint inhibitor (ICI)-related genes. We then compared the differences between the two classes of genes in the TCGA-BRCA transcriptome data to determine whether the immunotherapy response varied between groups with high and low lipid metabolism scores.

To assess the immune response, we used the tumor immune dysfunction and exclusion (TIDE) tool (http://tide.dfci.harvard.edu/login/) (30). In addition, we integrated the transcriptome-level data of each cell line from the cancer cell line encyclopedia (CCLE), the clinical trials research platform (CTRP), the PRISM database and the area under curve (AUC) data of 981 drugs (31). The difference in the AUC between different lipid metabolism score groups and the correlation between the AUC and lipid metabolism were analyzed to further identify therapeutic drugs with different scores.

To further analyze the drug of choice from the last set of experiments and the composition of lipid metabolism genes encoding proteins, we used the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) to obtain information on drug and lipid metabolism, mainly protein structure files. Then, batch processing and AutoDock Vina v.1.2.2 were used for molecular docking to select the receptor and ligand pairs with the minimum binding energy.

To further analyze the differentiation status between subsets, we used the classical monocle2 package (30) to perform a pseudotime analysis of T-cell subsets. The single-cell data were processed by constructing monocle objects, normalizing, and filtering low-quality cells and other processes. Highly discrete genes were selected to reduce the dimensionality of the data via the DDRTree method, and then the data were subjected to pseudotime analysis of different types of cells and important genes.

We utilized data from the human protein atlas (HPA) (32) to corroborate the variations in key gene expression between human BRCA tissues and normal tissues.

The ATCC provided the MCF10A, MDA-MB231, BT549, and SUM149PT cells. MCF10A, MDA-MB231, BT549, and SUM149PT cells were acquired from the Type Culture Collection of the Chinese Academy of Sciences in Shanghai, China. Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% bovine calf serum (HyClone) was used to culture MDA-MB231, and SUM149PT cells. Ten percent bovine calf serum (HyClone) was added to complete Roswell Park Memorial Institute 1640 (RPMI 1640) medium to culture the BT549 cells. DMEM/F12 supplemented with 10% bovine calf serum (HyClone) was used to culture the MCF10A cells. Every cell line was grown in a humidified incubator with 5% CO2 at 37°C.

Follow the manufacturer’s guidelines to extract total RNA from cells using TRIzol (Invitrogen). Using an Applied Biosystems TaqMan reverse transcription reagent kit, 1 mg of RNA per sample was utilized to synthesize cDNA. Using the SYBR Premix Ex Taq kit for real-time PCR (TaKaRa), qRT−PCR was conducted using an Applied Biosystems 7500 real-time PCR system. The primer sequences are shown in Supplementary Table S1 .

Samples of adjacent matched nontumor tissues and human BRCA tissues were acquired from the Affiliated Cancer Hospital of Harbin Medical University. IHC labeling was performed on nontumor and BRCA samples from 50 female subjects. The percentage of positive cells was used to assess tissue section staining quantitatively.

The cells were lysed on ice with RIPA lysis buffer. The proteins were separated via SDS−polyacrylamide gel electrophoresis, transferred to a PVDF membrane, blocked with 5% skim milk powder and washed with PBST. The membranes were subsequently incubated with the primary antibody at 4°C overnight. The membranes were incubated with secondary antibody for 1 hour after washing with PBST, and the protein bands were detected with enhanced chemiluminescence (ECL) luminescent solution.

The cells were plated in 96-well plates at 3000 cells per well. After 24 hours, the control group was treated with DMSO, and the experimental group was treated with different concentrations of niclosamide dissolved in DMSO. Cell viability was assessed via CCK-8 kit, and the light absorption value was measured at 450 nm with microplate reader.

Antibodies against the following proteins were used for the IHC and western blot experiments: ACSL1 (Proteintech, 13989-1-ap, IHC: 1:200, western blot: 1:1000), ACSF2 (Proteintech, 16140-1-ap, IHC: 1:200, western blot: 1:1000), CPNE3 (Proteintech, 11186-1-ap, IHC: 1:200, western blot: 1:1000), MTMR9 (ABclonal, A13124, IHC: 1:200, western blot: 1:1000), and β-actin (Santa Cruz, western blot: 1:1000). The following reagents were used in the CCK8 experiments: niclosamide (CAS No. 50-65-7, MedChemExpress, BAY2353) and Cell Counting Kit-8 (LABLEAD, CK001).

In this work, R software (https://www.www.r-project.org/, version 4.2.1) was used for all data calculations and statistical analyses. If not otherwise specified, correlation analyses were performed via Spearman’s correlation analysis via the Hmisc function of the base package of R software. The Wilcoxon rank sum test was used to compare MSI scores, altered portions of the genome, and mutation counts between the Score+ and Score- groups. The C-index, cross-validation and model ranking were used to determine the best machine learning mix. We used the Wilcoxon rank sum test to compare differences in gene transcription and methylation levels between the high-rated and low-rated groups. To control the error finding rate of multiple tests, we also adopted the Benjamini−Hochberg correction method. In addition, the differences of immune-related gene expression and immunoinfiltration indexes between the high-low rating groups were statistically analyzed via the Wilcoxon rank sum test. All of the statistical P values were bilateral, and differential gene screening was statistically significant if the corrected p value was < 0.05 (ns stands for p value>0.05; * stands for 0.01< p value <0.05; ** stands for 0.001< p value <0.01; *** stands for p value <0.001; **** stands for p value <0.0001).

Our overall experimental design is shown in Figure 1 . We scored human BRCA tumor tissue samples in the TCGA-BRCA dataset based on the lipid metabolism genes obtained from the MSigDB database. We selected the top 20% of the samples as the Score+ group and the bottom 20% of the samples as the Score- group. First, we performed GSVA separately for the Score+ group and the Score- group ( Figure 2A ), where a t value of GSVA greater than 1 was taken as the cut-off value. We found that the samples in the Score+ group were significantly enriched in the Adipogenesis, Heme metabolism, Bile acid metabolism, and Fatty acid metabolism gene sets. The Score- group was significantly enriched in the E2f target, Myc target v2, Myc target v1 and G2M checkpoint gene sets. On the basis of this grouping, we compared whether the microsatellite instability (MSI) score ( Figure 2D ), fraction of the genome altered ( Figure 2C ) and mutation count ( Figure 2B ) differed between the two groups. We found that the MSI score, fraction of the genome altered and mutation count in the Score- group were noticeably greater than those in the Score+ group.

We selected the top 20% of the samples as the Score+ group and the bottom 20% of the samples as the Score- group. The Score+ group included 217 samples, and the Score- group also included 217 samples. We used the rank-sum test to analyze differences between the high- and low-lipid metabolism groups. To screen for significant DEGs, we set the thresholds to |logFC| > 0.5 and adjusted P value < 0.05; ultimately, we identified 259 genes related to lipid metabolism that were differentially expressed between the Score+ and Score- (high- and low-score) groups ( Figures 3A, B ).

To construct a prognostic model focused on lipid metabolism-related genes, we initially employed univariate Cox regression analysis to pinpoint 26 genes associated with prognosis from the pool of 259 DEGs (p value < 0.05). Then, we utilized a combination of machine learning algorithms. We carried out parameter tuning for each basic algorithm, and different combinations were generated under different parameter settings. As seen from the specific combinations listed on the left side of Figure 4 , not only were there multiple parameter combinations for each individual algorithm, but pin-to-pair combinations and even multiple combinations between algorithms were also explored. Through these parameter adjustments, different ratio feature selection strategies, cross-validation strategies, and the use of integration methods, 184 model combinations were generated from 9 basic machine learning algorithms. Next, we evaluated the performance of 184 combinations of the machine learning algorithms using two datasets, GSE88770 and GSE20711, as validation sets. The C-index was used as the main model evaluation metric. The results showed that the machine learning combination of the elastic network (alpha=0.4) combined with LASSO had the highest C-index, so we selected this model as the best model ( Figure 4 ). Notably, TCGA-BRCA has a larger C index than the other two validation sets do, which may be related to the larger data volume and low noise of the data, and the greater amount of clinical information contained in it may provide additional context for the model, further enhancing its predictive power.

To assess the influence of this score on the overall survival (OS) of BRCA patients, we utilized the optimal cutoff value to categorize patients into two groups according to the LMPGS. The group with a high score for LMPGS had a notably lower OS rate than did the group with a low score in the training set ( Figure 5A , p value < 0.0001). The same conclusion was drawn for the GSE88770 and GSE20711 validation sets ( Figure 5B , p value = 0.12; Figure 5C , p value =0.014). To further validate the precision of the prognostic model, time−ROC curves and Kaplan−Meier (KM) survival curves were plotted for each of the above three datasets. We found that the AUC was greater than 0.9 in the TCGA-BRCA cohort, greater than 0.6 in the GSE88770 test cohort for the first 5 years, and greater than 0.6 in the GSE20711 test cohort for most of the first 10 years (Figures 5D–F).

To compare the accuracy of our lipid metabolism-related prognostic gene score model and other score models in predicting overall prognosis in BRCA patients, we first performed a multivariate analysis including the lipid metabolism-related prognostic gene score and other clinicopathological features using the TCGA-BRCA dataset. The lipid metabolism-related prognostic gene score, age, and tumor stage were identified as prognostic factors for BRCA patients ( Figure 6A ). A nomogram model was constructed ( Figure 6B ), and we employed a multivariate Cox regression model to assess the prognostic lipid metabolism-related gene score, age, and tumor stage. We also conducted a literature survey and searched for articles from September 1, 2022, to September 1, 2023, with keywords such as TCGA, prognosis and BRCA, and a total of 5 relevant articles were collected. We combined the AUC scores of these models and the number of genes included. In the comparison of the 1-, 3-, and 5-year survival data, our model had the highest AUC of the multiple models, and at the same time, the number of genes included in our model was relatively large ( Figure 6C ).

To explore the clinical distinctions between groups characterized by high and low scores for prognostic genes linked to lipid metabolism, we initially conducted a comparative analysis of age ( Figure 7A ), tumor stage ( Figure 7B ), and tumor status ( Figure 7C ). We found a significant difference in tumor stage (p = 0.004) between the high- and low-score groups founded on LMPGS model, whereas no notable disparities were observed in age or tumor status between these groups. We next assessed the effect of the lipid metabolism-related prognostic gene score on OS in patients with 31 additional tumor types. The results revealed that the lipid metabolism-related prognostic gene score was associated with OS in adrenocortical carcinoma, bladder urothelial carcinoma, cholangiocarcinoma, colon cancer, and other tumors, suggesting that the lipid metabolism-related prognostic gene score on OS has broad applicability ( Figure 7D ).

To further understand the differences between patients with high and low scores, we performed a transcriptomics analysis. First, we compared genomic variations to compare the TMB and mutation, deletion, and amplification profiles between the two groups ( Figure 8A ). Next, we focused on determining the differences between the two groups in terms of the three most important genes that composed the prognostic model of genes. In terms of methylation levels, the methylation level of ACSL1 was significantly greater in the low-score group (p value<0.01), while the methylation level of MTMR9 was greater in the high-score group (p value<0.01), and the methylation level of ACSF2 was not significantly different between the two groups ( Figure 8B ). At the transcriptome level, ACSL1, ACSF2 and MTMR9 were highly expressed in the high-score group ( Figure 8C ). In addition, the transcriptome expression levels of these three genes were also significantly correlated with the OS of BRCA patients. Compared with those of the low-expression group, the survival rates of the high-expression groups of ACSF2 ( Figure 8D , p value=0.023) and ACSL1 ( Figure 8F , p value=0.0017) were significantly better. Patients in the low-MTMR9 expression group had a considerably better survival rates than those in the high-MTMR9 expression group ( Figure 8E , p value=0.0015).

To assess the differences in immune cell infiltration between the two score groups, we used the ESTIMATE method (Figure 9C). The findings indicated that the group with a low lipid metabolism-related prognostic gene score had higher stromal ( Figure 9A ) and immune ( Figure 9B ) scores, while the high lipid metabolism-related prognostic gene score group of had greater tumor purity ( Figure 9D ). These findings suggest that the low-score group may exhibit a greater degree of immune cell infiltration than the high-score group does. Moreover, we further analyzed the differences in immunogenic cell death-related genes between the lipid metabolism-related prognostic gene score groups and found that ANXA1, IFNE, LRP1 and other genes were highly expressed among the LMPGS in the low-score group. CALR, EIF2AK1, EIF2AK1 and other genes were highly expressed in the LMPGS high -score group, and the majority of genes linked to immunogenic cell death differed between the two groups ( Figure 9E ). Next, we examined the differences in ICI-related genes between the prognostic gene groups associated with lipid metabolism. TNFRSF9, CD200, PDCD1, IDO1 and other genes were highly expressed in the low lipid metabolism-related prognostic gene score group, and the expression of approximately half of the ICI-related genes differed between the two groups ( Figure 9F ).

The TIDE scores of the two groups varied, indicating that the two groups have different responses to immunotherapy ( Figure 10A ). In addition, the percentage of patients who responded to immunotherapy predicted by the TIDE database varied significantly ( Figure 10B ). To examine in more detail how groups with high and low lipid metabolism-related prognostic gene scores respond differently to chemotherapy, we used the first-line BRCA chemotherapy drug Plinabulin for validation. The findings indicated that the group scoring high exhibited greater AUC values, indicating good drug absorption ( Figure 10C ). Subsequently, we analyzed the drug data from the CCLE and PRISM databases and found that the sensitivity to 7 drugs (CERANIB-2, YK-4-279, NICLOSAMIDE, MONENSIN, CUDC-907, NOCODAZONE, and CETRIMONIUM) was significantly correlated with the lipid metabolism-related prognostic gene score ( Figure 10D , p value <0.05). Moreover, there were notable differences in the dose−response curves between the two score groups (p value <0.05, Figure 10E ).

To further examine the possibility of interactions between these seven possible medications and the proteins encoded by the three most significant genes (ACSF2, ACSL1, and MTMR9) in lipid metabolism-related prognostic gene set, we acquired structural information on the proteins produced by these three genes as well as the molecular structures of the pharmacological ligands from the PubChem database. We subsequently performed molecular docking via AutoDock Vina v.1.2.2 to identify receptor−ligand pairs with binding free energies less than -5 kcal/mol. The results showed that only ACSF2, MTMR9 and NICLOSAMIDE could bind to each other (Figures 11A–F). These findings may help elucidate the interactions between drugs and target proteins and provide important clues for further research and drug design.

We performed quality control of the GSM545720 BRCA dataset according to quality control standards and then performed dimensionality reduction and clustering. We subsequently set a resolution value of 0.8 as the clustering criterion and collected the signature markers of each cluster from published single-cell studies. Finally, we identified eight subpopulations: epithelial cells, endotheliocytes, fibroblasts, NK cells, T cells, plasma cells, CD4 cells, and macrophages ( Figure 12A ).

We subsequently performed differential expression analysis via the FindAllMarkers function using the Wilcoxon method and found that among the 22 genes related to lipid metabolism and prognosis, CPNE3 was upregulated in epithelial cells ( Figures 12B, D ). To further analyze the changes in CPNE3 expression during tumor progression, we performed a pseudotime series analysis of epithelial cells. In the initial stage, the expression of CPNE3 gradually decreased over time ( Figure 12C ). These results contribute to our understanding of the single-cell population and variations in gene expression associated with lipid metabolism and prognosis at the cellular level (Figure 12E).

To validate the disparities in the expression of crucial genes between human BRCA tissues and normal tissues, we obtained immunohistochemical profiles of key genes from the Human Protein Atlas database. We detected marked differences in the immunohistochemical staining of two key proteins, ACSF2 and ACSL1, between tumor tissue samples and normal tissue samples. As shown in the figure, the expression levels of ACSL1 and ACSF2 were higher in tumor tissues than in normal tissues ( Figures 13A–D ).

According to the qRT−PCR results, the mRNA expression of MTMR9 and CPNE3 was greater in triple-negative breast cancer (TNBC) cell lines (MDA-MB-231 and BT549) and ER-, PR-, and HER2+ BRCA cell line (SUM149PT). However, ACSL1 and ACSF2 mRNA expression was greater in normal cell lines (nontumorigenic breast cell line: MCF10A) ( Figure 14A ). The MCF10A cell line is considered the baseline for evaluating changes in expression levels. IHC demonstrated the upregulation of ACSF2, MTMR9, and CPNE3 in BRCA tissues, whereas ACSL1 expression was not significantly different between tumor tissues and adjacent tissues ( Figures 14B, C ). Furthermore, a correlation analysis between the positive areas of each protein in tissue samples and progression-free survival (PFS) time was performed among 42 patients whose complete follow-up data were available. The results revealed that the expression levels of MTMR9 and CPNE3 were negatively correlated with the survival of BRCA patients. MTMR9, CPNE3 and Ki67 expression levels were positively correlated, and the difference was statistically significant ( Figures 14D, E ). However, no significant correlations were found between ACSL1 or ACSF2 expression and survival.