Development and Validation of a Prognostic Classifier Based on Lipid Metabolism-Related Genes for Breast Cancer.

Background: The changes of lipid metabolism have been implicated in the development of many tumors, but its role in breast invasive carcinoma (BRCA) remains to be fully established. Here, we attempted to ascertain the prognostic value of lipid metabolism-related genes in BRCA.
Methods: We obtained RNA expression data and clinical information for BRCA and normal samples from public databases and downloaded a lipid metabolism-related gene set. Ingenuity Pathway Analysis (IPA) was applied to identify the potential pathways and functions of Differentially Expressed Genes (DEGs) related to lipid metabolism. Subsequently, univariate and multivariate Cox regression analyses were utilized to construct the prognostic gene signature. Functional enrichment analysis of prognostic genes was achieved by the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). Kaplan-Meier analysis, Receiver Operating Characteristic (ROC) curves, clinical follow-up results were employed to assess the prognostic potency. Potential compounds targeting prognostic genes were screened by Connectivity Map (CMap) database and a prognostic gene-drug interaction network was constructed using Comparative Toxicogenomics Database (CTD). Furthermore, we separately validated the selected marker genes in BRCA samples and human breast cancer cell lines (MCF-7, MDA-MB-231).
Results: IPA and functional enrichment analysis demonstrated that the 162 lipid metabolism-related DEGs we obtained were involved in many lipid metabolism and BRCA pathological signatures. The prognostic classifier we constructed comprising SDC1 and SORBS1 can serve as an independent prognostic marker for BRCA. CMap filtered 37 potential compounds against prognostic genes, of which 16 compounds could target both two prognostic genes were identified by CTD. The functions of the two prognostic genes in breast cancer cells were verified by cell function experiments.
Conclusion: Within this study, we identified a novel prognostic classifier based on two lipid metabolism-related genes: SDC1 and SORBS1. This result highlighted a new perspective on the metabolic exploration of BRCA.

Introduction

According to 2020 GLOBOCAN statistics of A Cancer Journal for Clinicians, breast invasive cancer (BRCA) is the most common type of cancer worldwide, comprising approximately 11.7% of new cancer cases.1 Despite advances in screening, research, and treatment, BRCA incidence continues to increase.2,3 BRCA continues to be a major public health problem globally. The risk factors for BRCA are complex and numerous, such as postpartum childbirth, obesity, metabolic syndrome, life stress, and lack of exercise.4 However, the correlation between these risk factors and the biological behavior of cancer cells has not yet been studied in depth, resulting in low value for the accurate prediction of the prognosis.5,6 Among these risk factors, the role of lipid metabolism disorders in the occurrence and development of BRCA has started to be appreciated.7–9 Establishing a more accurate model for predicting the prognosis of BRCA patients based on biomarkers related to BRCA lipid-metabolism is expected to become a valuable new research direction.

Changes in lipid metabolism and signal transduction have been recognized as one of the signs of abnormal cell growth and cancer progression.10 In addition to being a component of the cell membrane and a source of energy storage, lipids can also activate cell proliferation, differentiation, and migration, as well as induce cell death.11,12 Abnormal lipid metabolism has been reported to be related to the occurrence and development of a variety of malignant tumors, such as prostate cancer,13 ovarian cancer,14 pancreatic cancer,15 liver cancer,16 and BRCA.17 Studies have found that changes in lipid metabolism exist during the early stages of BRCA development, and a lipid-rich environment promotes the proliferation and migration of BRCA cells.18,19 Some bioinformatics studies have found that tumor gene sets and inflammation-related gene sets are highly expressed in adipocytes adjacent to BRCA, while histocytology studies have found that the lipid metabolism of residual cancer cells in patients undergoing neoadjuvant therapy for BRCA differs from that before treatment and in normal tissues.20,21 These findings suggest that genes related to lipid metabolism may become targets for BRCA treatment. Changes in lipid metabolism may also affect the target organ tropism of BRCA metastasis. It was found that the BRCA cells that are more prone to brain metastasis have significant changes in lipid metabolism.22 Although current studies have determined that abnormal lipid metabolism is related to BRCA progression and resistance to treatment, there is currently no specific diagnosis and treatment plan incorporating this information due to the complexity of the effects of lipid metabolism-related genes on BRCA.23–25

The use of bioinformatics to analyze the targets of disease diagnosis and treatment has been widely used.26,27 In this study, we analyzed and obtained the biological processes of differential lipid metabolism genes with prognostic value in BRCA based on The Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Molecular Signatures Database (MSigDB) databases. We identified two genes, SDC1 and SORBS1, which we screened, analyzed, and verified as biomarker genes that may play important roles in BRCA. In addition, we not only verified the expression changes of these two genes in BRCA samples and breast cancer cell lines (MCF7, MDA-MB-231), but also verified the functions of these two genes in breast cancer cells by cell function experiments. Finally, we constructed a prognostic classifier and preliminarily verified it with datasets in the GEO database and clinical follow-up results we collected. The validation results showed that our prognostic classifier has a clinical value and can provide a new reference for the diagnosis and treatment of BRCA patients. In addition, we used Connectivity Map (CMAP) and the Comparative Toxicogenomics Database (CTD) to investigate compounds targeting prognostic genes, and to provide new insights for the synthesis of new drugs.

Materials and Methods

Cell Lines

MDA-MB-231 and MCF-7 cells were obtained from ATCC (Shanghai, China) and cultured in DMEM (Invitrogen, Carlsbad, CA, USA) at 37°C under 5% CO2.

Clinical Sample Acquisition

Tumor tissue and normal breast tissue adjacent to the tumor were collected from 30 BRCA patients. All the patients were introduced to surgery at the First Affiliated Hospital of Zhengzhou University, Henan, China, from December 2020 to January 2021 and did not receive any anti-cancer treatment before surgery. All patients were between 18 and 80 years old and agreed to accept the postoperative follow-up plan. Tissue specimens were collected within 30 min after surgery and quickly frozen in liquid nitrogen. Postoperative monitoring and treatment continued in accordance with relevant consensus guidelines. The degree of tumor differentiation was graded according to the WHO grading system. Overall survival (OS) was defined as the time between surgery and death. The study was approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University and conducted according to the principles expressed in the Declaration of Helsinki. Written informed consent was obtained from each patient prior to collection of any samples. Detailed information of each patient, including serial number, sampling time, pathology type, ER positivity rate, HER2 status, Ki-67 expression, lymph node stage, tumor size, age, and side of affected breast, is available in .

In order to verify the clinical applicability of our prognostic model, we collected the clinical data and pathological sections of 50 patients with BRCA who underwent surgery in our hospital from January 2015 to March 2015. The plan for clinical patient follow-up and collection of patient clinical data to verify the prognosis model was approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University. We followed up by phone or in-person to collect the survival information of the 50 patients. All patients were between 18 and 80 years old and agreed to accept the post-operative follow-up plan. The personal information of each patient is strictly confidential, and as such, written informed consent was obtained from each patient or its entrustment. The clinical conditions and follow-up results of the patients are shown in .

Data Sources

We acquired RNA expression data from 1171 samples (1072 breast cancer samples and 99 normal samples) from TCGA database. Of the 1072 BRCA samples, only 1069 BRCA samples recorded complete survival information, so candidate prognosis-related lipid metabolism differentially expressed genes (DEGs) were identified based on these samples. In addition, only 881 of the 1069 BRCA samples had complete clinical information, which was needed for the sample to be used in the construction of the prognostic classifier. The patients’ characteristics were presented in .

The GSE109169 dataset28 () was downloaded from the GEO database (25 BRCA samples, 25 normal samples) and was utilized to legitimize the expression pattern of candidate prognosis‐related lipid metabolism DEGs. The patients’ characteristics were presented in . Furthermore, we downloaded the GSE2068529 dataset from the GEO database, which was utilized as independent external validation set to authenticate the prognostic efficacy of the multigene signature. The GSE20685 dataset contains 327 BRCA samples with complete survival information and matched transcriptomic data (). The characteristics of the BRCA samples in the GSE20685 dataset were available in . Moreover, 1499 genes related to lipid metabolism were obtained by de-duplication from 13 gene sets related to lipid metabolism, which were downloaded from the Gene Ontology (GO) and the Molecular Signatures Database (MSigDB) websites. These 13 gene sets were GO_GALACTOLIPID_METABOLIC_PROCESS, GO_GLYCEROLIPID_METABOLIC_PROCESS, GO_GLYCEROPHOSPHOLIPID_METABOLIC_PROCESS, GO_GLYCOSPHINGOLIPID_METABOLIC_PROCESS, GO_LIPID_METABOLIC_PROCESS, GO_LIPOPROTEIN_METABOLIC_PROCESS, GO_MEMBRANE_LIPID_METABOLIC_PROCESS, GO_NEUTRAL_LIPID_METABOLIC_PROCESS, GO_PHOSPHOLIPID_METABOLIC_PROCESS, GO_REGULATION_OF_LIPID_METABOLIC_PROCESS, GO_REGULATION_OF_LIPOPROTEIN_METABOLIC_PROCESS, GO_REGULATION_OF_PHOSPHOLIPID_METABOLIC_PROCESS, GO_SPHINGOLIPID_METABOLIC_PROCESS. All lipid metabolism-related genes are integrated into .

BRCA-FPKM data were obtained from the TCGA database and log2 transformed for the gene expression files before use; the expression profiles acquired from the GEO database were log2 transformed. Moreover, we used the normalize Between Arrays function of the limma package to normalize the gene expression profiles of GSE109169 and GSE20685 datasets.

Translation level validation of identified prognostic genes was performed using the Human Protein Atlas30 database (), an online database containing immunohistochemical expression data for approximately 20 of the most common cancers.

Differential Analysis

Transcriptome RNA-seq data of 1171 cases (normal samples, 99 cases; BRCA samples, 1072 cases) were downloaded from TCGA database. R package limma31 was used to perform differentiation analysis of the gene expression, and DEGs were generated by the comparison between the BRCA samples vs the normal samples (). DEGs with |log2 fold change | > 1 and adjusted P < 0.05 were considered significant.

Ingenuity Pathway Analysis (IPA)

Overlapping genes () for DEGs and lipid metabolism-related genes were analyzed by the IPA software (version 1–19-00). RNA sequencing results were imported into IPA, and the enrichment status of canonical pathways and diseases and functions were assessed by P-values.32,33 Z scores were indicated by orange or blue to denote activation or inhibition of pathways. In IPA, terms with P < 0.05 were considered significantly enriched; furthermore, in IPA-Canonical pathway analysis, the activation and inhibition of the enriched pathway were assessed by Z-score, specifically, Z-score > 0 and P < 0.05 indicated that the pathway was active, while Z-score < 0 and P < 0.05, then the pathway was considered to be in the inhibited state.

Identification of Candidate Prognosis-Related Lipid Metabolism DEGs

Lipid metabolism DEGs were included in the univariate Cox regression analysis, and genes with a P < 0.05 were selected to be used in the multivariate Cox regression analysis. Similarly, genes with a P < 0.05 were considered as candidate prognosis-related lipid metabolism DEGs. Subsequently, GO and KEGG analyses were carried out based on the above genes using the R software clusterProfiler package. GO analysis consisted of three components: biological process (BP), molecular function (MF) and cellular component (CC). Entries with an adjusted (adj.) P < 0.05 were considered significantly enriched.

Knock-Down and Overexpression Experiments of SDC1 and SORBS1

The empty plasmids used for constructing controls of SDC1 gene overexpression and SORBS1 gene overexpression (oe-SDC1control and oe-SORBS1control), the non-targeting plasmid used for constructing controls of SDC1 gene knockdown and SORBS1 gene knockdown (si-SDC1 control and si-SORBS1 control), SDC1 knockdown plasmids, SORBS1 knockdown plasmids, SDC1 overexpression plasmids and SORBS1 overexpression plasmids were designed and constructed by Genechem (Shanghai, China). The sequences of the plasmids used in this study are provided in . According to the manufacturer’s instruction, cells of different groups were transfected with these plasmids using Lipofectamine 2000. 72 hours after transfections, the effects of knockdown or overexpression were examined by qRT-PCR and Western blot.

RNA Isolation, Reverse Transcription, and Real-Time PCR

Total RNA was isolated from 30 paired tissues and 8 groups of experimental cells using TriQuick Reagent (Solarbio, Shanghai, China) according to the manufacturer’s instructions. Then the concentration and purity of the RNA solution was quantified using a NanoDrop 2000 nucleic acid protein quantifier (Thermo Fisher Scientific, Waltham, MA, USA). qRT-PCR was performed as described previously.34 Briefly, the extracted RNA was reverse-transcribed to cDNA using the FastQuant RT Kit with gDNA eraser (TIANGEN, Beijing, China) prior to qRT-PCR. The reverse transcription reaction consisted of 4 µL of 5x reaction buffer, 1 µL of oligo (dt) primer, 1 µL of random hexamer primer, 1 µL of Servicebio RT Enzyme Mix, and 0.1 ng–5 µg total RNA, and finally the reaction was brought to 20 µL using RNase-free water. The qRT-PCR reaction consisted of 2 µL of reverse transcription product, 10 µL of 2X SYBRGreen qPCR Master Mix (High ROX, Servicebio, Wuhan, China), 0.4 µL each of forward and reverse primer, and 7.2-µL nuclease-free water. PCR was performed in a MiniAmp Thermal Cycler (A37834, Thermo Fisher Scientific, Waltham, MA, USA) under the following conditions: 95℃ for 3 s, followed by 40 cycles of 95℃ for 15s and 60℃ for 30s. The GAPDH gene served as an internal control. 2-ΔCt method was used to calculate the RNA levels of tumor samples, paired adjacent samples and 8 groups of experimental cells. Primer sequences used for qRT-PCR are shown in Table 1.

Table 1

Primer Sequences Used for qRT-PCR

Gene	Primer	Sequence (5’→3’)	Length	Tm	GC %
SDC1	Forward Primer	ACGGCTATTCCCACGTCTC	19	54.58	57
	Reverse Primer	TCTGGCAGGACTACAGCCTC	20	56.54	60
SORBS1	Forward Primer	CACAATCGAGAACAGCAAAAACG	23	56.05	43
	Reverse Primer	ACCCGCCTACTGTCATCCTTT	21	56.95	52
GAPDH	Forward Primer	CTGGGCTACACTGAGCACC	19	55.46	63
	Reverse Primer	AAGTGGTCGTTGAGGGCAATG	21	57.03	52

Western Blot Analysis

Western blot analysis was conducted using an SDS-PAGE electrophoresis system (Bio-Rad Laboratories, Hercules, CA, USA). Briefly, the total protein content was extracted from tissue or 8 groups of experimental cells using a RIPA buffer (P0013C, Beyotime, Shanghai, China). Protein samples were separated by SDS-PAGE and transferred onto nitrocellulose membranes (88,520, Thermo Fisher Scientific, Waltham, MA, USA), which were subsequently blocked 12 h at 4℃ with 5% skimmed milk containing TBST (Tris-buffered salt solution, containing 50 mmol/L Tris-HCl, 150 mmol/L NaCl, 0.1% v/ v Tween-20, pH 7.4) solution. Antibodies against SDC1 (ab128936; Abcam, Cambridge, MA, USA), SORBS1 (ab224129; Abcam, Cambridge, MA, USA), and GAPDH (T0004; Affinity Biologicals, Cincinnati, OH, USA) were used as primary antibodies. The samples were incubated with horseradish peroxidase–conjugated secondary antibodies at 37℃ for 1 h. The membrane was imaged using an Amersham Imager 600 (GE Healthcare UK Limited, Little Chalfont, UK).

Immunohistochemistry

In the present study, we performed immunohistochemical analysis on 30 BRCA tissue samples and their paired cancer-adjacent normal breast tissue samples. The tissues were collected and fixed in 4% paraformaldehyde, embedded in paraffin, and sectioned at 6-µm intervals. The sections were placed in EDTA (pH 9.0) for antigen repair, washed with PBS (pH 7.4), treated with 3% H2O2, blocked with goat serum, and then incubated 12 h at 4℃ with SORBS1 (1:500; ab224129, Abcam, Cambridge, MA, USA) and SDC1 (1:1600; 10593-1-AP, Proteintech, Rosemont, IL, USA). Subsequently, the secondary antibody goat anti-rabbit IgG (1:200; GB23303, Servicebio) was incubated for 1 h at 37℃, and the positive sites were labelled with a DAB (diaminobenzidine) color development solution (G1211, Servicebio, Wuhan, China). Finally, hematoxylin staining (G1004, Servicebio, Wuhan, China) was performed to visualize the nuclei. Images (200×) were captured with a microscope (XSP-C204, COIC, Chongqing, China), and three different visual fields were analyzed. Image-Pro Plus 6 processing software (Media Cybernetics, Inc., Rockville, MD, USA) was used to analyzed results and to calculate average optical density,35 where AOD=integrated OD/measurement area.36 The higher the AOD, the higher the expression of the protein.37,38 The person imaging the slides was uninformed during the experiment and analysis.

Invasion Assays

Cell invasion assay was done performed according to the manufacturer’s instructions. Briefly, the upper surface of a Transwell (8-μm pore size; Corning, Corning, NY, USA) was coated with 20μL diluted Matrigel (BD Biosciences, Franklin Lakes, NJ, USA). Equal numbers of the indicated cells were seeded in the upper chamber of the Transwell in serum-free medium, with FBS added to the lower chamber. After incubation for 48 hours, non-invading cells were removed with a cotton swab, and the invaded cells were stained with crystal violet. Average invaded cell number per field of view was obtained from 5 random fields.

Scratch Assays

Cells were plated in 24-well cell culture plates (1 × 105 cells/well) followed by overnight incubation. When culture was confluent (24 h post transfection), a p-200 pipette tip was used to score two vertical lines and one horizontal line (average width: 700–900 μm) simulating a “wound” by scratching the culture. All the wells were washed with sterile PBS twice and reference markings were drawn near the scratch area from the bottom side of the plate with a fine tip marker. Scratch images were captured within the marked area using an inverted microscope (Laxco, Mill Creek, WA, USA) with a 10X magnification objective piece.

Prediction Model

We randomly divided TCGA database of 881 BRCA samples containing complete clinical and survival information into a training set (n = 617) and a testing set (n = 264) based on a 7:3 ratio. The characteristics between the BRCA samples in the TCGA-training and -testing sets were not significantly different (), indicating that they could be used for subsequent analysis. In the training set, univariate and multivariate Cox regression analyses were applied to extract the optimal prognosis-related lipid metabolism DEGs based on selected candidate prognosis-related lipid metabolism DEGs. Ultimately, genes with a P < 0.05 were used to construct a prognostic gene signature. Prognosis-related lipid metabolism DEGs were visualized on the forest plots. Risk scores were simultaneously calculated for each patient in the training set, testing set, and external validation set of the following equations: risk score = (expression of Gene1 × β1Gene1) + (expression of Gene2 × β2Gene2) + (… expression of Genen × βnGenen).39,40 β was represented by the regression coefficient, which was generated by the step multivariate Cox proportional hazards regression model. Based on this score, patients were classified into high- and low-risk score groups on the basis of the median classification method.

Prognostic Analysis

Risk curves were plotted with the pheatmap package. The R survival package41 was employed to assess the relationship between lipid metabolism-related gene signature and OS (P < 0.05) and to plot Kaplan-Meier (K-M) curves. Subsequently, the survival ROC package42,43 was executed to create ROC curves and calculate AUC in order to estimate the predictive accuracy of the prognostic gene signature. The above analyses were performed simultaneously in the training set, testing set, and the external validation set (GSE20685 dataset). Univariate and multivariate independent prognostic analyses were developed by survival software packages to confirm whether risk scores could be applied as the independent prognostic indicator. The rms44 and survival packages of R were performed to construct a nomogram. Subsequently, calibration curves were plotted to evaluate the agreement between actual and predicted survival.45 Decision curve analysis (DCA) was carried out to calculate the net clinical benefit for each model. The optimal model was the one with the greatest net benefit calculated.46

Connectivity Map (CMap) Analysis

CMap is designed to reveal the relationships among genes, compounds, and biological conditions, which is a systematic, data-driven process.47 We resorted to CMap02 in order to screen potential compounds that might target prognostic genes. Compounds with |enrichment score| ≥ 0.5 and a P < 0.05 were selected as potential therapeutic drugs for BRCA. The compounds were then further filtered by CMap’s CLUE tool () for related Mechanism of Actions (MoA) and their inhibitors in order to investigate their joint intrinsic mechanism of action.

Predicting Prognostic Gene-Chemotherapy Drug Interaction Network by CTD

A network of interactions between prognostic genes and chemotherapeutic agents was constructed using the Comparative Toxicogenomics Database (CTD) to obtain chemotherapeutic agents that could reduce or increase prognostic gene expression levels. Briefly, chemotherapeutic drugs for SORBS1 and SDC1 were searched for in the CTD, and the prognostic gene-drug interaction networks were visualized by using Cytoscape 3.8.2 ().

Statistical Analysis

The Jvenn website () was used to produce a Venn diagram for intersection analysis. A K–M curve with Log rank test was used to assess the OS differences between different groups. An ANOVA test was performed to detect the association of the risk score with clinical characteristics, and also, to reveal the differences in the levels of risk scores within different subtypes of clinical characteristics. The statistical analyses were conducted using the R software. A P < 0.05 indicated statistically significant differences.

Model Validation by Detailed Data from 2015.01 to 2015.03

After collecting clinicopathological and follow-up data and conducting immunohistochemistry staining, correlations between SDC1 and SORBS1 expression and survival prognosis were analyzed.

Results

Identification of Differentially Expressed Lipid Metabolism-Related Genes

RNA sequencing data of 1072 BRCA samples and 99 normal samples were extracted from TCGA database. A differential analysis based on the R package limma was performed with normal samples as controls in order to screen for genes that were aberrantly expressed in BRCA, and genes satisfying |log2 FC| > 1 and adjusted P < 0.05 were identified as DEGs. We identified a total of 1732 BRCA-related DEGs, of which 644 genes were upregulated and 1088 were downregulated in the BRCA group compared with the normal group (Figure 1A; ). Subsequently, to recognize lipid metabolism-related DEGs, we performed intersection analysis based on a list of 1732 DEGs and 1499 lipid metabolism-related genes. The results were shown in Figure 1B, and a total of 162 overlapping genes were identified, which were defined as lipid metabolism-related DEGs (). Further, the heatmap demonstrated the expression pattern of these lipid metabolism-related DEGs between 1072 BRCA and 99 normal samples in the TCGA database (Figure 1C), including 34 up-regulated genes and 128 down-regulated genes (BRCA vs normal). Canonical pathway analysis of these lipid metabolism-related DEGs in IPA revealed that several pathways related to lipid metabolism (eg, fatty acid β-oxidation I, fatty acid α-oxidation, and fatty acid activation), carcinogenesis (eg, AMPK signaling, JAK/STAT signaling, and ERK/MAPK signaling), and immune response (eg, MIF regulation of innate immunity, IL-7 signaling pathway, and glioblastoma multiforme signaling) were inhibited. Furthermore, we demonstrated that these genes were correlated with the negative regulation of estrogen-dependent BRCA signaling, HER2 signaling in BRCA, and BRCA regulation by stathmin1 (; ). The diseases and functions analysis confirmed that lipid metabolism and small molecule biochemistry were the most abundant pathways (; ).

Figure 1

Identification of metabolism-related genes and construction of a prognostic classifier. (A) A volcano plot of all DEGs is shown combined with |log2FC| and an adjusted p-value. Red represents 644 upregulated DEGs. Green represents 1088 downregulated DEGs. (B) Extraction of metabolism-related genes from the DEGs. (C) Heatmap of lipid metabolism-related DEGs between tumor and matched adjacent normal tissue. Different colors represent the expression trend of lipid metabolism-related DEGs in the two groups, red represents high expression of genes and blue represents low expression of genes.

Functional Enrichment Analysis of Candidate Prognosis-Related Lipid Metabolism DEGs

In the present analysis, we included 1069 BRCA samples that contained complete survival information. Univariate Cox analysis screened 76 genes with P < 0.05 from the identified 162 lipid metabolism-related DEGs (Table 2). Thirteen candidate prognosis-related lipid metabolism genes were obtained from the subsequent multivariate Cox analysis. Through examination with a forest plot, genes with HR > 1 were found to include LEPR, FGF2, GPAM, SDC1, CCDC3, PCK2, PLTP and HSD17B13, which may be risk factors for BRCA prognosis; while ABCD2, FABP4, SORBS1, FAM126A and GP1HBP1 with HR < 1 could be termed as protective factors for BRCA prognosis (Figure 2A).We then validated the expression patterns of the 13 genes mentioned above using the GSE109169 dataset. The results showed that the expression trends of these genes were consistent with those in the TCGA database. Except for SDC1, the expression of the other 12 genes was downregulated in the BRCA samples (Figure 2B and C). The underlying features of these genes were ascertained by GO and KEGG pathway analysis (). Undoubtedly, several GO terms related to lipid metabolism were found in biological processes, such as “positive regulation of lipid metabolic process,” “regulation of lipid biosynthetic process,” “regulation of lipid metabolic process,” “positive regulation of lipid biosynthetic process,” “glycerolipid metabolic process,” and “triglyceride metabolic process” (). Furthermore, “lipid droplet” and “external side of plasma membrane” were the most remarkably enriched GO items of CC (). As implied by the KEGG pathway analysis (), candidate prognosis-related lipid metabolism genes were primarily enriched along two pathways-The proliferator-activated receptor (PPAR) signaling pathway and adipocytokine signaling pathway ().

Table 2

Univariate COX Analysis of Different Lipid Metabolism Genes (P < 0.05)

ID	HR	HR.95L	HR.95H	p-value
LEPR	1.080015	1.0446	1.11663	6.04E-06
ADM	1.023029	1.012557	1.033609	1.44E-05
ADH1C	1.05337	1.025811	1.081671	0.000121
CCDC3	1.010499	1.005017	1.01601	0.000168
PLTP	1.005733	1.002641	1.008835	0.000274
PCK1	1.02409	1.010228	1.038143	0.000619
FGF2	1.037334	1.015572	1.059562	0.000703
G0S2	1.001482	1.000617	1.002347	0.000783
MGLL	1.01302	1.005283	1.020816	0.000943
HCAR2	1.042764	1.016626	1.069574	0.001225
ACACB	1.015828	1.00608	1.025671	0.001412
CD36	1.002847	1.001089	1.004609	0.001494
KLF4	1.008258	1.003137	1.013405	0.001546
ADH1B	1.001686	1.000613	1.002761	0.00207
DGAT2	1.007378	1.00267	1.012108	0.002099
PTGIS	1.025369	1.009068	1.041933	0.002183
LGALS12	1.011648	1.004028	1.019325	0.00268
CAV1	1.002493	1.000862	1.004126	0.002729
LRP1	1.007875	1.002627	1.01315	0.003231
TEK	1.054722	1.017956	1.092815	0.00325
ENPP2	1.015092	1.00498	1.025306	0.003363
RETSAT	1.005345	1.001763	1.00894	0.003423
ANXA1	1.003369	1.001071	1.005673	0.004047
GPAM	1.005001	1.001526	1.008489	0.004761
PLIN1	1.001213	1.00037	1.002057	0.004782
HSD17B13	1.090967	1.026844	1.159093	0.004846
ADIPOQ	1.001466	1.000426	1.002506	0.005711
ACSL1	1.00275	1.000788	1.004715	0.005984
BMX	1.168945	1.045744	1.30666	0.006012
CYP2U1	1.154779	1.042081	1.279664	0.00602
LIPE	1.003918	1.001119	1.006724	0.006048
AADAC	1.074246	1.020586	1.130726	0.006155
LPL	1.001335	1.000377	1.002294	0.006318
VAV3	0.989454	0.981893	0.997073	0.006749
FABP4	1.000299	1.000082	1.000517	0.007006
ANGPTL4	1.011239	1.003002	1.019543	0.007398
PDE3B	1.038544	1.010131	1.067755	0.007536
MLXIPL	1.049762	1.012711	1.088168	0.008075
CIDEA	1.006311	1.001635	1.011008	0.008109
EGR1	1.000754	1.00019	1.001317	0.00872
QKI	1.056527	1.013967	1.100873	0.008763
PDGFRA	1.017727	1.004373	1.031259	0.00912
CEBPA	1.012191	1.002986	1.02148	0.009333

Abbreviations: COX analysis, Cox proportional hazards model; HR, hazard ratio.

Figure 2

Functional enrichment analysis of candidate prognosis-related lipid metabolism DEGs (A) Forest map of lipid metabolism-related DEGs associated with prognosis. (B) Expression of prognostic lipid metabolism DEGs in TCGA, where the horizontal axis represents different gene, the vertical axis represents the gene expression distribution, where different colors represent different groups. Asterisks represent levels of significance, ****p < 0.0001. (C) Expression of prognostic lipid metabolism DEGs in GSE109169, where the horizontal axis represents different gene, the vertical axis represents the gene expression distribution, where different colors represent different groups. Asterisks represent levels of significance ***p < 0.001, ****p < 0.0001.

Filtering Candidate Prognosis-Related Lipid Metabolism DEGs for Prognostic Prediction

Two prognosis-related genes, SDC1 (HR = 1.2638, P = 0.03076) and SORBS1 (HR = 0.7628, P = 0.03576), were retained based on the univariate Cox regression analysis in the training set (Figure 3A). Subsequently, we proceeded to risk model gene selection utilizing multivariate Cox regression (Figure 3B). The results indicated that SDC1 (P = 0.037) and SORBS1 (P = 0.043) were the optimal variables for the construction of prognostic signature. According to the hazard ratio (HR), SDC1 might be a promoter of BRCA (HR = 1.244, 95% CI: 1.014–1.523) and SORBS1 was a promising tumor suppressor gene (HR = 0.768, 95% CI: 0.595–0.991). Subsequently, based on the TCGA-training set, we revealed the association of the two identified prognostic genes with OS in BRCA patients by K-M curves. The relatively short OS of high expression of SDC1-BRCA compared to low expression of SDC1-BRCA (Figure 3C; P = 0.048); and the expression of SORBS1 was positively associated with OS in BRCA patients, with high expression of SORBS1 implying longer OS (Figure 3D; P = 0.0096), which evidence corroborated the inferences made based on HR values. Next, the expression values of prognostic genes and multivariate Cox regression coefficients () were utilized to generate risk scores for individual samples in the training set. Moreover, IPA exposed a sophisticated network of interactions involving the two prognostic genes ().

Figure 3

Risk score independent of the prognostic analysis. (A) Univariate Cox analysis showing the hazard ratio of each candidate prognosis-related lipid metabolism DEG in predicting overall survival in BRCA from the training set. (B) Multivariate Cox regression analysis of 2 prognostic genes. (C) Survival analysis of SDC1. (D) Survival analysis of SORBS1. (E) The distribution of risk score and survival status of the training set-BRCA patients in the high- and low-risk groups. (F) Survival analysis of the high- and low-risk groups in the training set. (G) ROC curve showing the moderate accuracy of the constructed prognostic model of BRCA from the training set. (H) The distribution of ROC analysis of the two-gene signature in the testing set. (I) The distribution of Kaplan-Meier survival analysis of the two-gene signature in the testing set. (J) The distribution of risk score and survival status analysis of the two-gene signature in the testing set. (K) The distribution of ROC analysis of the two-gene signature in the GSE20685 dataset. (L) The distribution of Kaplan-Meier survival analysis of the two-gene signature in the GSE20685 dataset. (M) The distribution of risk score and survival status analysis of the two-gene signature in the GSE20685 dataset.

The 2-gene based prognostic signature was assessed in the training set. Depending on the training set-median risk score, the 617 samples were divided into high (n = 308) and low (n = 309) risk groups, and the risk score curves and survival status of all samples were presented in Figure 3E. K-M survival analysis showed that the risk score based on the 2-gene prognostic signature could significantly discriminate the clinical outcomes of BRCA patients, with a low-risk score implying a better prognosis (P = 0.003; Figure 3F). Subsequently, ROC curve analysis showed that the risk score predicted the AUC of 1, 3, and 5-year OS for training set-BRCA patients to be 0.694, 0.681, and 0.618, respectively (Figure 3G). Further, we evaluated the general applicability of the risk score system in the testing set and external validation set (GSE20685 dataset). The performance of risk score in the TCGA-testing set (n = 264; Figure 3H) and the independent external validation set (GSE20685 dataset, n = 327; Figure 3I) was similar to that in the TCGA-training set. High-risk scores were significantly associated with poor prognosis (all P < 0.05). The AUC of the ROC curve at 1, 3, and 5 years was 0.635, 0.628, and 0.611 in the testing set, and the AUCs in the GSE20685 dataset were 0.819, 0.672, and 0.650, respectively. This evidence suggested that the SDC1 and SORBS1-based prognostic signature possessed robust prognostic predictive power. Additionally, heat maps were constructed to visualize the expression of prognostic genes in each cohort (), with SDC1 being relatively highly expressed in the high-risk group and SORBS1 being relatively overexpressed in the low-risk group. Moreover, immunohistochemical results, pathological results, and follow-up results of 50 patients () collected from January 2015 to March 2015 were divided into high-risk and low-risk groups according to this prognostic model, and the difference in the survival curves between the two groups was found to be statistically significant by survival analysis (P = 2.516e-02) ().

The Obtained Two-Gene Signature Was Linked to Pathological Features of BRCA

We further investigated whether the genetic signature was implicated in the pathological features of BRCA. Specifically, in the AJCC pathologic T-stage subgroup (Figure 4A), risk score levels were significantly higher in the T2/T4 subtype than in the T1 subtype; interestingly, risk score levels were markedly lower in the T3 subtype compared with the T2 subtype; compared with the T3 subtype, risk score levels were remarkably higher in the T4 subtype; however, risk score levels were comparable between the T1 and T3 and T2 and T4 subtypes. In the AJCC pathologic N-stage subgroup (Figure 4B), risk score levels were increased considerably in the N1 and N2 subtypes relative to the N0 subtype; however, risk score levels were not statistically different between the N0 and N3, N1 and N2/N3, and N2 and N3 subtypes. In the AJCC pathologic stage subgroup (Figure 4C), the risk score levels were proportional to the stage level, with the higher the stage level, the higher the risk score level; however, there was no significant difference in the risk score levels among stage II, stage III, and stage IV. We also assessed the relationship between different types of BRCA and risk scores (Figure 4D). Of the 4 types of BRCA, HER2-enrich BRCA had the highest risk score, followed by luminal B, basal-like, and luminal A BRCA. Such results suggested a substantial interaction between lipid metabolism genes and clinical molecular features.

Figure 4

Scatter dot plot shows the association between risk score of TCGA breast cancer samples and clinical characteristics. (A) Risk scores in the different pathologic T stages of BC. (B) Risk scores in the different pathologic N stages of BC. (C) Risk scores in the different pathologic stages of BC. (D) Risk scores of different BC subtypes. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Furthermore, we performed stratified survival analysis in the TCGA-training set, TCGA-test set, and independent external validation set, respectively, to explore whether the prognostic value of the risk score was applicable to other clinical factors. All patients were divided into designated subgroups based on clinical characteristics in the corresponding datasets. BRCA samples in the TCGA-training set and -test set were analyzed stratified according to age, gender, pathologic T-stage, pathologic N-stage, pathologic M-stage, and pathologic tumor stage. In the TCGA-training set, the 2-gene signature was useful in >60 years, female, T2, T3, N2, N3, M0, and stage III subgroups, with clinically and statistically significant prognostic value (). In the TCGA-test set, the prognostic model was able to significantly differentiate the clinical outcomes of patients in the <60 years, female, N3, M1, and stage IV subgroups (). BRCA samples in the independent external validation set were analyzed stratified according to age, gender, metastatic events, regional recurrence or not, adjuvant chemotherapy, pathologic N stage, pathologic M stage, pathologic tumor stage, and subtype. In the independent external validation set, the 2-gene signature had clinically and statistically significant prognostic value in the <60, female, no regional recurrence, N0, and type I subgroups ().

The Independent Utility of the Prognostic Classifier

First, univariate Cox regression analysis indicated that age, tumor stage, AJCC pathologic stage, AJCC pathologic TNM stage, and risk score were correlated with OS (Figure 5A). Subsequent multivariate Cox regression analysis suggested that age, AJCC pathologic stage, and risk score could serve as independent prognostic factors (Figure 5B). A nomogram was derived based on the Cox regression analysis results that employed four independent prognostic factors to predict the prognosis of patients with BRCA in TCGA database (Figure 5C). The concordance index of the nomogram was 0.739. The calibration curves suggested that the nomogram (combined model) was possibly under- or overestimating the mortality (Figure 5D). Decision curve analysis (DCA) evidenced that the combined model exhibited the best net returns for the 5-year OS as opposed to 3-year OS (Figures 5E and F). Regrettably, we were unable to capture the net returns of the combined model on 1-year OS, possibly owing to the inadequacy of the 1-year OS data. Collectively, these results indicated that the nomogram constructed using a combined model may be the optimal nomogram in predicting long-term survival (5 years) for patients with BRCA as compared to the nomogram obtained using individual prognostic factors, which may assist in clinical management.

Figure 5

The Independent Utility of the Prognostic Classifier. (A) Univariate Cox regression analysis of clinical features and the prognostic signature. (B) Multivariate Cox regression analysis of clinical features and the prognostic signature. (C) Nomogram to predict survival probability at 1, 3, and 5 years. (D) The calibration curves for the nomogram. (E) The 3-year decision curve for the nomogram and other clinical traits. (F) The 5-year decision curve for the nomogram and other clinical traits.

Validation of Prognostic Genes Based on Clinical Samples

We performed qRT-PCR analysis to assess the expression levels of the two genes used to construct our prognostic model. Consistent with the results of the biochemical analysis, SORBS1 was significantly downregulated in the tumor samples as compared to the adjacent healthy samples (Figure 6A), while SDC1 was markedly overexpressed in the tumor samples (Figure 6B). Furthermore, we performed Western blot analysis and found that protein expression of SORBS1 was downregulated in BRCA (Figures 6C), and the protein expression level of SDC1 was upregulated in BRCA (Figure 6D). The results of immunohistochemistry showed positive staining for SORBS1 and SDC1 proteins in a brown or tan color, which revealed the localization of SORBS1 and SDC1 to be in the cytoplasm/membrane of BRCA tissue and adjacent normal breast tissue, respectively (Figure 6E and F). The AOD of SORBS1 was higher in the adjacent healthy group as compared to the BRCA group (Figure 6G; P < 0.05). The AOD of SDC1 protein was significantly higher in the BRCA group than in the adjacent healthy group (Figure 6H; P < 0.05). Besides, we explored the protein expression levels of these prognostic genes in BRCA using the HPA database. Consistently, protein levels of SDC1 were not expressed in normal breast tissue, while moderate expression levels of this gene were observed in breast cancer tissue (); meanwhile, moderate protein expression levels of SORBS1 were observed in normal breast tissue, while low protein expression levels of these genes were observed in breast cancer tissue (). In conclusion, the present results suggested that the transcriptional and translational expression levels of SDC1 were overexpressed in BRCA patients while SORBS1 were overexpressed in normal tissues.

Figure 6

Expression levels of SORBS1 and SDC1 proteins in clinical samples. (A) Expression levels of SORBS1 measured by qRT-PCR analysis. (B) Expression levels of SDC1 measured by qRT-PCR analysis. (C) Representative image of Western blot analysis of SORBS1 protein. N, Adjacent healthy samples; T, BRCA samples. (D) Representative image of Western blot analysis of SDC1 protein. (E) Immunohistochemical staining (SORBS1 protein) of BRCA and adjacent healthy groups; magnification, ×200. (F) Immunohistochemical staining (SDC1 protein) of BRCA and adjacent healthy groups; magnification, ×200. (G and H) The AOD of SORBS1 and SDC1 protein in BRCA and adjacent healthy groups. Mean ± SEM; n = 30; ***P < 0.001 vs adjacent healthy; ****P < 0.0001 vs adjacent healthy.

CMap Analysis Identified Novel Candidate Compounds Targeting the Prognostic Genes

To identify potential compounds capable of targeting prognostic lipid metabolism-related genes (SDC1 and SORBS1), we screened a total of 76 compounds that met the |enrichment score| > 0.5, P < 0.05 by the CMap02 database (). Unfortunately, only 37 of the 76 compounds could be traced in the CLUE tool. The 37 compounds that were able to inhibit prognostic gene expression as described above are illustrated in . Subsequent CMap MoA analysis revealed 44 mechanisms of action common to the above compounds. Amitriptyline was involved in four MoAs, including norepinephrine inhibitor, norepinephrine reuptake inhibitor, serotonin receptor agonist, and serotonin reuptake inhibitor, while the MoAs involving acetylcholine release stimulant, acetylcholinesterase inhibitor, butyrylcholinesterase inhibitor, and potassium channel antagonist were all modulated by the compound tacrine. Furthermore, a total of 11 compounds shared the following six mechanisms: 1) adrenergic receptor agonist, 2) bacterial cell wall synthesis inhibitor, 3) glucocorticoid receptor agonist, 4) histone deacetylase (HDAC) inhibitor, 5) progesterone receptor agonist, and 6) sterol demethylase inhibitor.

Prognostic Gene-Drug Interaction Network Analysis

Interactions between prognostic genes and drugs used for cancer treatment were probed through the CTD database and visualized via Cytoscape. As depicted in , there were 50 and 42 drugs that affected the expression of SDC1 and SORBS1, respectively, of which 16 drugs regulated the expression of both prognostic genes. For instance, abrine reduced the SDC1 expression levels, while upregulating the SORBS1 expression ().

SDC1 Could Promote the Migration and Invasion of Breast Cancer Cells

To verify transduction efficacy, RT-PCR and Western blot were used to detect the expression of SDC1 following 24 h of transduction (Figure 7A and B). To explore the effect of SDC1 on breast cell invasion and migration, transwell invasion assays and wound healing assays were performed. Transwell invasion assays demonstrated that si-SDC1 could reduce the invasion ability of MCF-7 and MDA-MB-231 cells and oe-SDC1 could enhance the invasion ability of MCF-7 and MDA-MB-231 cells in vitro (Figure 7C and D). The wound healing assays demonstrated that the migratory abilities of MCF-7 and MDA-MB-231 cells were significantly suppressed following transfection of si-SDC1 (Figure 7E and F). Similarly, compared with the corresponding negative control groups, the speed of wound closure was significantly faster in the MCF-7 and MDA-MB-231 cells transfected with SDC1 plasmid (Figure 7E and F). These results suggest that SDC1 can increase the invasiveness of breast cancer cells.

Figure 7

Effects of SDC1 knockdown or overexpression on BC cells’ invasion and migration. (A) Quantitative PCR analysis of SDC1 expression levels in si-SDC1 control, si-SDC1, oe-SDC1 control and oe-SDC1 breast cancer cell lines (MCF-7 and MDA-MB-231). (B) Western blot of SDC1 expression levels in si-SDC1 control, si-SDC1, oe-SDC1 control and oe-SDC1 breast cancer cell lines (MCF-7 and MDA-MB-231). (C) Effects of SDC1 knockdown or overexpressed on invasion of MDA-MB-231 breast cancer cells. (D) Effects of SDC1 knockdown or overexpressed on invasion of MCF-7 breast cancer cells. (E) Effects of SDC1 knockdown or overexpressed on migration of MDA-MB-231 breast cancer cells showed by scratch assays. (F) Effects of SDC1 knockdown or overexpressed on migration of MCF-7 breast cancer cells showed by scratch assays. ***:P < 0.001.

SORBS1 Could Inhibit the Migration and Invasion of Breast Cancer Cells

To verify transduction efficacy, qRT-PCR and Western blot were used to detect the expression of SOBRS1 following 24h of transduction (Figure 8A and B). In the transwell invasion assays, the invasion ability of MCF-7 and MDA-MB-231 cells in the si-SORBS1 group was significantly increased compared with si-SORBS1 control group, and the invasion ability of MCF-7 and MDA-MB-231 cells in the oe-SORBS1 control groups was significantly increased compared with oe-SORBS1 groups (Figure 8C and D). The wound healing assays demonstrated that the migratory abilities of MCF-7 and MDA-MB-231 cells were significantly promoted following transfection of si-SORBS1 (Figure 8E and F). Similarly, the speed of wound closure was significantly slower in the MCF-7 and MDA-MB-231 cells transfected with SORBS1 plasmid comparing with the negative control groups (Figure 8E and F). Taken together, SORBS1 inhibits the migration and invasion of breast cancer cells, and SORBS1 depletion conversely promotes the migration and invasion of breast cancer cells.

Figure 8

Effects of SORBS1 knockdown or overexpression on BC cells’ invasion and migration. (A) Quantitative PCR analysis of SORBS1 expression levels in si-SORBS1 control, si-SORBS1, oe-SORBS1 control and oe-SORBS1 breast cancer cell lines (MCF-7 and MDA-MB-231). (B) Western blot of SORBS1 expression levels in si-SORBS1 control, si-SORBS1, oe-SORBS1 control and oe-SORBS1 breast cancer cell lines (MCF-7 and MDA-MB-231). (C) Effects of SORBS1 knockdown or overexpressed on invasion of MDA-MB-231 breast cancer cells. (D) Effects of SORBS1 knockdown or overexpressed on invasion of MCF-7 breast cancer cells. (E) Effects of SORBS1 knockdown or overexpressed on migration of MDA-MB-231 breast cancer cells showed by scratch assays. (F) Effects of SORBS1 knockdown or overexpressed on migration of MCF-7 breast cancer cells showed by scratch assays. **:P < 0.01; ***:P < 0.001.

Discussion

The high annual incidence rate of BRCA has a serious impact on human health and the social economy.48 Detection of BRCA at an early stage is a critical step for successful treatment and improvement of the prognosis.49 The transformation of metabolic pattern is one of the important characteristics of tumor cells. A number of studies have found that tumor cells will undergo lipid metabolism reprogramming to meet the needs of rapid proliferation.50,51 Because the tumor microenvironment is rich in adipokines, abnormal lipid metabolism plays a more important role in breast cancer than other malignant tumors. In recent years, great achievements have been made in abnormal lipid metabolism in tumors.12,19 Lipid molecules with the diagnostic potential for BRCA are constantly being discovered.12,52 Studies have found that the abnormal lipid metabolism of BRCA cells is closely related to their resistance to HER2 inhibitors and CDK4/6 inhibitors.53 In addition, some studies have reported changes of phospholipids in the plasma, serum,54,55 and urine56,57 of BRCA patients. However, lipid metabolism involves both exogenous and endogenous processes, and its impact on BRCA risk and prognosis needs to be further determined. Thus, lipid metabolism-related genes in BRCA may be a breakthrough point that can further improve the prognosis of breast cancer patients.

In this study, a risk assessment classifier closely related to the prognosis of BRCA was constructed through the analysis of differentially expressed lipid metabolism-related genes with prognostic value in BRCA and revealed the correlation between BRCA and lipid metabolism. We not only used the GEO database, but also the cell function experiments and the clinical data of tissue samples to verify that this prognostic classifier can offer good clinical application value.

In addition, we analyzed the correlation between the prognostic classifier and pathological conditions and clinical characteristics and found that the risk of TNBC was significantly higher. TNBC is histologically defined as lacking estrogen and progesterone receptors (ER/PR) and human epidermal growth factor receptor 2 (HER2) amplification, it is a subtype of BRCA with the worst prognosis and highest risk because of the lack of effective therapeutic targets.58,59 However, no targeted drugs for TNBC have been discovered so far.60,61 At present, the main systemic treatment for TNBC is still chemotherapy.62 Antibody drug conjugates (ADCs), immune-oncology (I/O) therapies, and polyadenosine diphosphate-ribose polymerase (PARP) inhibitors have joined the armamentarium against specific types of TNBC but have limited durable efficacy in TNBC patients.63 Our findings may provide information for the development of drugs that target lipid metabolism genes, thus, making it possible to find new pathways in the treatment of TNBC.

In order to further explore the mechanism of the lipid metabolism DEGs related to BRCA survival, we performed network analysis using IPA and KEGG and GO enrichment analysis on the 13 genes selected. We found that PPAR signaling pathways and adipokine signaling pathways had the highest correlation with these genes, and IPA uncovered a significant enrichment of lipid metabolism pathways and estrogen-dependent BRCA signaling. Moreover, the effect of PPAR signaling pathway and estrogen-dependent signaling in BRCA was previously confirmed and the related drugs are in clinical application,64–67 so these results verified our research’s directions and methods are scientific from the other side. The top GO terms involved in lipid metabolism were “lipid droplet” and “external side of the plasma membrane.” Adipokine signaling pathways are the main pathways of lipid metabolism, and their relevance to tumors has been confirmed in many studies.68–71 They play an independent and combined role in the activation of intracellular signaling networks related to the malignant phenotype of BRCA cells.72 The GO term “external side of the plasma membrane” has also been reported to be associated with testicular cancer.73 Lipid droplets74 are the natural immune center that integrates cell metabolism and host defense and are gradually being recognized as a prominent feature of many cancers.75,76 Studies have reported that lipid droplet and the PPAR signaling pathway are closely related to the metastasis of TNBC.77–79 Furthermore, lipid droplets have been found to have increased expression in drug-resistant TNBC cells.21 In addition, based on the current progress of immunotherapy for BRCA and the increasing recommendation of immunotherapy in the treatment regimen of TNBC,80 although our functional enrichment analysis did not indicate that the lipid metabolism-related DEGs were related to immune cells, we still calculated the correlation between SDC1 and SORBS1 and immune cells based on TIMER database, but no positive result was obtained. This means that SDC1 and SORBS1 do not affect the prognosis of BRCA through immune function pathway, which also suggests that it may be a new direction for us to explore the diagnosis and treatment of BRCA from the factors related to lipid metabolism. These findings have a certain correlation with the poor prognosis of TNBC in the prognostic classifier we established in this study and support the development of therapies for TNBC targeting lipid metabolism.

The genes we used for prognostic classifier constructing were SORBS1 and SDC1. SORBS1 encodes a CasitasB.1ineagelymphoma (CBL)–associated protein that functions in the signaling and stimulation of insulin, and it also been confirmed to be related to COPD.81–83 Studies have found that SORBS1 can participate in immune-related gene signal transduction, and its increased expression is closely related to poor prognosis in cervical cancer,84 gastric cancer,85 prostate cancer,86 and colorectal cancer.87 However, the exact role and diagnostic significance of SORBS1 in BRCA are not yet clear. Syndecan-1 (SDC1, also known as CD138) is a key regulator of fatty acid synthesis and catalyzes the formation of mono-unsaturated fatty acids (MUFA) from saturated fatty acids (SFA).88–91 SDC1 has been implicated in various cancers, including BRCA.92–95 The inactivation of SDC1 can lead to the consumption of unsaturated fatty acids, thereby inhibiting cancer cells.91,96,97 However, the imbalance in the ratio of saturated to unsaturated fat can lead to a cellular stress response and death in both normal and transformed cells.91,98 Therefore, although SDC1 has potential as a biomarker, its complex relationship with the body limits its clinical application.99–101 In this study, we constructed a prognostic classifier based on SORBS1, which is currently less studied, and SDC1, which has a complex effect on the body. In addition, we examined the correlation between invasiveness and expression of SDC1 and SORBS1 in MCF-7 and MDA-MB-231 breast cancer cell lines by cell function experiments. This is the first prognostic classifier that combines the two genes in order to jointly assess their prognostic risk in BRCA and this classifier can change the limited application status of SDC1 and SORBS1.

Through the CMap02 database, we screened potential drugs that target genes related to lipid metabolism as a means to control the progression of BRCA. We also used the CTD database to explore the interaction between selected prognostic genes and cancer treatment drugs. We found that vorinostat, abrine, and NSC 689534 are very promising drugs. Vorinostat is a simple molecule histone deacetylase (HDAC) inhibitor which has both oral and vena formulation. The main indication currently approved by the FDA is the treatment of skin T-cell lymphoma and its combination with pembrolizumab (a monoclonal anti-PD1 antibody) can increases the response rate of head and neck malignancies to immunotherapy.102 Whereas, pembrolizumab combined with standard chemotherapy plays an important role in the treatment of advanced TNBC, and this may mean that we can try adding vorinostat to this treatment to increase the tumor response rate to pembrolizumab. In addition, the main metabolic pathway of vorinostat is via hydroxamic acid glucuronidation and oxidative cleavage of the aliphatic methylene chain.67 This confirms our findings that this drug may be effective in treating BRCA patients carrying the prognosis-related lipid metabolism genes we screened. Abrine is an alkaloid, and the main component of the traditional Chinese medicine cochinchinensis.103 In addition, abrine has a structure that can form hydrogen bonds with peroxisome PPARs, thereby affecting the PPAR pathway.104 Our analysis found that abrine reduces the expression level of SDC1 while increasing the expression of SORBS1, which indicated that it may have a dual protective effect on BRCA patients. The inhibitory effect of abrine on BRCA cells has been confirmed in several in vitro experiments.105,106 It is our hope that our research results can prompt experts in the field of chemistry and pharmacy to pay attention to this compound with therapeutic potential. NSC 689534 (2-pyridinecarbaldehyde N, N-bis (2-pyridinylmethyl) thiosemicarbazone) is a member of the heterocyclic thiosemicarbazone family of compounds, and its biological activity is significantly different between the unbound state and the metal-bound state.107,108 In recent years, there have been reports about the effects of these drugs on human cancer, but most of them are pharmacological studies. It is found that the chelating ability of thiosemicarbazone NSC 689534 can induce oxidative/ER stress and inhibit tumor growth in vitro and in vivo.108 Our findings suggest that NSC 689534 is a potential drug for the treatment of BRCA which is worthy of further preclinical investigation.

Our analysis provides a reference for the prognosis of BRCA patients, a basis for the diagnosis and treatment of patients, and a new research direction for BRCA research. This study still has some certain limitations, these results are preliminary that mainly based on the secondary mining and analysis of previously published datasets. Although PCR experiments, Western blot experiments, immunohistochemistry experiments and cell function experiments were performed in this study to initially validate the association between the two genes and breast cancer, we still need to conduct prospective studies in order to further verify the clinical applicability of this prognostic classifier.

Conclusions

In summary, we have constructed a prognostic classifier of BRCA based on two lipid metabolism-related genes: SDC1 and SORBS1. This classifier has a value in predicting the disease-free survival rate of BRCA patients and identifying high-risk patients. Our results strengthen the underestimated role of abnormal tumor lipid metabolism in the prognosis of BRCA. The translational application of this classifier will guide clinicians to make a more informed decision regarding adjuvant systemic treatments and choosing follow-up plan.