{'fraction_non_alphanumeric': 0.028015100337770712, 'fraction_numerical': 0.010878203854559905, 'mean_word_length': 5.183353808353808, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 4, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of a Lipid Metabolism-Associated Gene Signature Predicting Survival in Breast Cancer

Background Cancer metabolism and specifically lipid metabolism play an important role in breast cancer (BC) progression and metastasis. However, the role of lipid metabolism-associated genes (LMGs) in the prognosis of breast cancer remains unknown. Methods The expression profiles and clinical follow-up information of 1053 BC were downloaded from The Cancer Genome Atlas (TCGA), and metabolic genes were downloaded from the Gene Set Enrichment Analysis (GSEA) dataset. Univariate Cox regression and least absolute shrinkage and selection operator (LASSO) regression analyses were performed on the differentially expressed metabolism-related genes. Then, the formula of the metabolism-related risk model was composed, and the risk score of each patient was calculated. The breast cancer patients were divided into high-risk and low-risk groups with a cutoff of the median expression value of the risk score, and the prognostic analysis was also used to analyze the survival time between these two groups. Finally, we analyzed the expression, interaction and correlation among the lipid metabolism-associated genes risk model. Results The results from the prognostic analysis indicated that the survival was significantly poorer in the high-risk group than in the low-risk group in TCGA, and single-sample gene set enrichment analysis (ssGSEA) shows it is plausible that lipid metabolism is highly correlated with tumor immunity. Conclusion Lipid metabolism-associated genes may become a new prognostic indicator predicting the survival of BC patients. The prognostic genes (n = 16) may help provide new strategies for tumor therapy.


Background
Female breast cancer has surpassed lung cancer as the most commonly diagnosed cancer. Despite the dramatic improvement in breast cancer prognosis due to recent therapeutic advances, such as more effective adjuvant and neo-adjuvant chemotherapies, together with more radical and safer surgery, advances in early diagnosis and treatment over the past decades, breast cancer prognosis is still very poor, and the death rate is as high as 6.9%. 1 The exploration of potential biomarkers and regulatory mechanisms for early diagnosis and therapeutic targets of BC have important scientific significance and application values.
Cancer is basically a disease of abnormal cell growth and proliferation: in order to meet its needs, cells synthesize nucleic acids, proteins and lipids at an accelerated rate. 2 Unlimited cell proliferation is one of the most important characteristics of tumor cells. Increased cancer cell proliferation requires the rapid synthesis of lipids for the generation of biological membranes and lipids that provide cancer cells with energy during times of nutrient depletion. The synthesis of lipids in healthy organisms is tightly controlled and responds to the nutrient's status of the cell. However, many human cancers display aberrant lipid metabolism, 3 and reprogrammed metabolism is considered a hallmark of cancer. 4,5 Next-generation sequencing capabilities are revolutionizing approaches to all fields of medicine: recent advances in our understanding of tumor biology have uncovered a growing list of clinically relevant biomarkers. 6 However, the vast majority of studies have concentrated mainly on a single gene, and its predictive ability is insufficient compared with multiple biomarkerbased models. 7 In clinical practice, the more accurately a patient's OS stage can be predicted, the sooner the clinician can make the clinical decision. Lipid decomposition and anabolism of tumor cells, which are different from normal cells, have attracted more and more attention, and it is also a new hotspot of tumor targeted therapy.
In this study, we established a robust 16-gene prognostic signature for breast cancer by integrating TCGA datasets, which might complement classical clinical prognostic characteristics, and further aid clinicians in personalized treatment planning.

Data Collection
The RNA sequencing (RNA-seq) expression data and clinicopathological information of female breast cancer patients from 1053 breast cancer tissue samples and 111 normal tissue samples were downloaded from the TCGA BRCA dataset (https://portal.gdc.cancer.gov/). Probes were transformed into corresponding Entrez gene names referring to the annotation files. In addition, 193 genes in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway associated with metabolism were also extracted from the "c2.cp.kegg.v7.0.symbols" gene sets in GSEA platform (https://www.gsea-msigdb.org/gsea/downloads. jsp). Pathways include: glycerolipid metabolism, glycerophospholipid metabolism, ether lipid metabolism and sphingolipid metabolism. The genes in pathways were extracted for further analysis. The detailed flow-process diagram of this study is shown in Figure 1.

Construction of a Lipid-Metabolism Model
First, the expression level of lipid metabolism genes was extracted from the total gene expression list. If a gene appeared more than once in the same sample, the limma of Bioconductor R package was utilized for averaging operations. 8 Second, the limma was utilized to identify differentially expressed lipid-metabolism-associated genes (DE-LMGs) between breast cancer tissue samples and normal breast tissue samples. The false discovery rate (FDR) threshold was set at FDR < 0.05 for DE-LMGs calling. To establish the lipid metabolism-related risk model, univariate Cox regression analysis was performed on the lipid metabolism-related genes. A total of 18 prognostic related differential genes were obtained by the intersection of DE-LMGs and prognostic genes. Heatmaps of the 18 genes were plotted using the "heatmap" R package. The mutation rates of prognostic DE-LMGs were analyzed by the cBioPortal for Cancer Genomics online website (https://www.cbioportal.org). 9 The OncoPrint schematic was constructed in cBioPortal (TCGA provisional) to directly reflect the mutation of 817 BRCA patients. A protein-protein interaction (PPI) network of proteins encoded by all overlapping DE-LMGs with prognostic value was visualized using String (http://string-db.org). 10 To avoid overfitting, the least absolute shrinkage and selection operator (LASSO) was utilized to select variables with high prognostic value. 11 Next, 1000 LASSO iterations were performed for the prognostic model construction using the "glmnet" package in R and their regression coefficients were obtained. Finally, the formula of the risk score was composed as follows, and risk scores were computed: Risk score = ∑ni = ∑Coefi × xi, where xi represents the normalized expression level of target gene i and Coefi represents the regression coefficient. According to the median risk score in TCGA dataset, 1014 patients in the data set were divided into high-risk and low-risk groups after samples with a survival time of zero were removed. Kaplan-Meier plots were used to evaluate survival differences between the high-and lowrisk groups. Receiver operator characteristic curves (ROC) and area under the curves (AUC) were used to evaluate the availability of the prognosis model via the "survivalROC" R package. To analyze the distribution differences between different groups, PCA was performed using the "prcomp" function in the STATS package in R.

Univariate and Multivariate Cox Regression Analysis
Univariate Cox regression analysis was presented for assessment of the prognostic values of the risk score and clinical features (age, stage, T classification, N classification, M classification). Then, multivariate Cox regression analysis was used to determine which prognostic factors could independently predict the survival of patients.

Functional Enrichment and Pathway Analysis
BC patients were divided into the high-and low-risk groups based on the median risk score, Gene Ontology (GO) 12 and Kyoto Encyclopedia of Genes and Genomes (KEGG) 13 pathway enrichment analyses for all selected DEGs between the two risks. Cohorts were performed with the "clusterProfiler" package in BioConductor using |log2FC|≥1 and FDR < 0.05 as thresholds. Then, we determined the scores of 16 tumor-infiltrating immune cells and 13 immunerelated functions for samples by ssGSEA.
Bioconductor R package "GSVA" was used to compare ssGSEA enrichment scores for immune cells and immune-related pathways between the two groups (ie high-or low-risk groups). 14 Finally, correlations between the risk signature and the key immune regulators, PD-L1 and PD-L2 were evaluated.

Validation of the Protein Expression Levels of LMGs in the HPA
The Human Protein Atlas (HPA) is a database that aims to map all the human proteins in cells, tissues, and organs using an integration of various omics technologies (https:// www.proteinatlas.org/). We also used the HPA database to verify the protein expression levels of LMGs based on immunohistochemistry.

Statistical Analysis
We used R software (version 3.6.2) to perform all statistical analyses. The Student's t-test was used to compare gene expression levels between BC samples and noncancer samples. The OS for the two risk groups was Log rank test. The correlation of prognostic model risk score with the key immune regulators was tested by Pearson analysis. The P value was adjusted by BH method. Unless otherwise stated, p<0.05 was considered statistically significant.

Identification of Prognostic LMGs in the Breast TCGA Cohort
Among the 193 LMGs, 116 genes were significantly differentially expressed between BC samples and adjacent breast tissue samples (FDR < 0.05); 21 genes were statistically significantly related to OS in the univariate Cox regression model. Finally, a total of 18 DEGs were related to OS and were selected as hub genes for further analyses (Figure 2A). Additionally, a heatmap showing the expression profiles of the 18 genes was plotted ( Figure 2B), and the prognosis of 18 genes are shown in Figure 2C. The PPI network for removing free nodes that indicates a tight interplay of lipid metabolism genes is shown in Figure 2D. More genes connected to CHPT1 than the other genes, indicating that CHPT1 may have a significant role in BC through a PPI network. Besides, the results of Figure 2E also showed a correlation network between co-expression levels of lipid metabolism genes in another manner. Mutational information of these 18 genes showed that amplifications and deep deletions were the most frequent mutation types ( Figure 2F). Three LMGs had mutation rates ≥ 4%, and LYPLA1 had the highest mutation rate.

Construction of a Prognostic LMGs Signature
According to LASSO, the optimal tuning parameter λ was identified based on a 1-SE (standard error). As a result, 16 prognosis-related key LMGs were identified and integrated to construct a prognostic signature for BC ( Figure 3A and B). Figure 3C shows the risk score distribution of patients in TCGA database. As we can see from PCA mappings ( Figure 3D), patients formed two distinct clusters. As the risk score increased, the patients' survival time decreased and the death risk increased. Patients in the high-risk group were more likely to die earlier than those in the low-risk group ( Figure 3E). The results from the Kaplan-Meier plot indicated that the survival was significantly poorer in the high-risk group than in the low-risk group ( Figure 3F)  Samples have been chosen from well-established clinical data. Among the 867 patients, we explored risk stratification for patients with different clinicopathologic factors, including risk parameters, age, histological grade and clinical stage. These variables indicated significant differences in univariate analysis and age, the stage, N classification, M classification showed significant differences in multivariable analysis. The risk score was significantly related to overall survival (OS) in both univariate and multivariate Cox regression analyses ( Figure 4A and B and Table 1).

Functional Enrichment Analyses
We performed GO and KEGG functional enrichment analyses on risk-related DEGs to analyze the relationship between biological functions and risk score. The results indicated that the DEGs mainly focused on immunoglobulin complex, lymphocyte mediated immunity and immune response-activating ( Figure 5A). KEGG functional enrichment analysis suggested that the DEGs were mainly related to hematopoietic cell lineage, cytokine-cytokine receptor interaction ( Figure 5B). To further explore the relationship between the BC prognosis and immune status, we further used ssGSEA to quantify the infiltrating scores of immune cell and immunity-related functions. The enrichment score of aDCs, DCs, iDCs, B cells, CD8+T cells, T helper cells, NK cells, TIL were higher in the low-risk group than in the high-risk group (adjusted p value<0.05, Figure 5C). Meanwhile, low-risk group had a higher score of CCR, HLA, check-point, cytolytic activity, APC co-inhibition, inflammation promotion, para inflammation, T cell costimulation, T cell co-inhibition, and type II INF response (adjusted p value < 0.05, Figure 5D). In terms of immune checkpoints, considering the roles of the immune checkpoint protein PD-L1 (also known as CD274) and PD-L2 (also known as PDCD1LG2) in immune evasion, the relationship between these sites and the LMGs signature was comprehensively analyzed. The gene expression levels of PD-L1 and PD-L2 in the high-risk subgroup were significantly lower than those in the low-risk subgroup ( Figure 5E and F). Expression levels were also significantly negatively correlated with LMGs-based risk scores ( Figure 5G and H).

Expression Levels of Key Genes in the Clinical Samples
We compared the mRNA and protein expression levels of 16 genes in clinical database. The result showed that the mRNA expression levels in TCGA were presented in Figure

Discussion
Previous studies have developed a number of tumour molecular signatures for breast cancer. 6 The assessment of estrogen receptor, progesterone receptor, and human epidermal growth factor receptor 2 is mandatory for all invasive breast tumors, and has been widely used in clinical practice. 6 However, outcomes of studying patients based on a single gene lack of accuracy. The lack of effective and reliable prognostic biomarkers still remains

9509
a major problem in improving the clinical outcomes of breast cancer patients. The relationship between lipid metabolism and cancer is multifaceted. Recent studies have shown that reprogramming of cellular lipid metabolism directly leads to malignant transformation and progression. 15,16 It broadened our understanding of how lipid metabolism was relevant to cancer biology. To our knowledge, this study is the first to use LMGs to predict the prognosis of BC. In this study, we comprehensively analyzed 193 LMGs in BC samples and their relationships with prognosis. In functional analyses, Differentially expressed genes between high-and low-risk groups were significantly enriched in several immune cell types and immune-related pathways. We identified 116 DE-LMGs that were significantly differentially expressed between BC samples and normal samples, supporting the prevalence of lipid metabolism in the pathogenesis and progression of BC. Moreover, we found 18 of these genes were related to OS in a univariate Cox regression analysis. These results demonstrated the importance of studying LMGs in breast tumors. The combination of univariate analysis and LASSO Cox regression was conducted to screen genes to indicate either poor or good prognosis and to construct a robust gene signature, which has been widely used in studies. [17][18][19] In the result, we established a new lipidrelated 16-gene signature. According to the results of the ROC curve, the AUC of the risk score established by 16 LMGs was more than 0.7. This result showed that the risk score established by these 16 LMGs had a high value for predicting prognosis.
In recent years, the tumor microenvironment (TME) has been a research hotspot. During BC progression, tumor immune microenvironment remodeling with the change of the ratio of immune cells and release of multiple immune inhibitory and reactive cytokines is a critical feature. 20,21 The role of the tumor microenvironment (TME) in breast cancer immunomodulation is vitally important, 22 and a better understanding of the immune cell infiltrate in the breast cancer microenvironment is crucial for the development of more effective therapeutic approaches. 23 By analyzing the pathways and functions of genes enriched in differences between the high-and low-risk groups, we observed a strong correlation between gene signature and immune function. Given the above situation, we intended to unearth more information about the immunological characteristics of the LMGs signature. As excepted, ssGSEA analysis revealed a higher proportion

9511
of CD8+T cells and NK cells in the low-risk group. CD8 +T cells are the main kind of cytolytic lymphocytes in tumor microenvironment, while NK cells, a type of cytotoxic lymphocytes, are crucial constituents of the innate immune system and play a key role in immune surveillance. 22 This phenomenon may indicate that CD8 +T cells and NK cells have the ability to target and kill tumor cells, consistent with previous studies, 24-28 which may be the reason why the low-risk group has a better prognosis than the high-risk group. Tumor infiltrating lymphocytes (TILs) can recognize and kill malignant tumor cells and are important endogenous inhibitors of tumor growth. These important effector cells often lose their ability to kill tumor cells due to exhaustion. Lipid metabolism reprogramming in TIL can be triggered by the TME, and lipids and their derivatives can impact TIL function. 29 PD-1 and PD-L1 constitute an essential inhibitory mechanism, which causes T cell exhaustion. This is the main reason why PD-L1 has drawn increasing attention from researchers concerned. [30][31][32] Pd-L1 and PD-L2 are key regulators of immune response. 33 This study also verified that pd-l1 and pd-l2 had significant differences in the two risk subgroups and were negatively correlated with risk score. There have been many studies of lipid metabolism effects on tumor immunity. Recently, a potent sphingolipid metabolite regulated tumorigenesis, and responded to chemotherapy and immunotherapy by affecting the trafficking, differentiation or effector function of tumorinfiltrating immune cells (TIICs), and implicated in many processes that are important for BC, 34 but the potential mechanism of LMGs and tumor immune microenvironment in BC remains unclear, and requires further investigation.

Conclusions
In summary, we established a robust 16-gene prognostic signature for breast cancer by comprehensive mining transcriptional profiles of lipid metabolism-associated genes. Our study provides a new understanding of the prognostic value of tumor progression in BC and will benefit the prognosis assessment of BC patients. We also confirmed that there was a close relationship between lipid metabolism and breast cancer microenvironment, and cancer cells may influence tumor infiltrating immune cells for immune escape through lipid metabolism regulation. Targeting lipid metabolism in cancer cells could have therapeutic benefits. However, there are several limitations to our study. First, the present research is a retrospective study based on TCGA databases. Therefore, prospective clinical research should be performed to validate the application of this model. Secondly, the underlying mechanisms of lipid metabolism-associated genes on BC, and its relation to tumor immune status remain relatively enigmatic and warrant further investigation.

Ethical Statements
TCGA belongs to public databases. The patients involved in the database have obtained Medical Ethics Committee of The Affiliated Hospital of Qingdao University approval. Users can download relevant data for free for research and publishing relevant articles. Our study is based on opensource data, so there are no ethical issues.