To assess the potential connection between cholesterol and the risk of breast cancer, genetic data on cholesterol (met-a-307, sample Size 7,813, number of SNPs 2545,608)and breast cancer(ieu-a-1132, ER+ Breast cancer (Oncoarray), sample size 833691, number of SNPs 10680275) were searched and obtained from the IEU Open GWAS project(https://gwas.mrcieu.ac.uk/). The data then were briefly collated and subjected to a two-sample Mendelian randomization (2-SMR) analysis. Mendelian randomization-Egger (MR-Egger) method analyses were the main way performed along with the inverse variance-weighted (IVW) method analysis, Weighted-median method analysis, Weighted mode method analysis and Simple mode method analysis (12).

The basic information on breast cancer RNA sequencing transcriptome data, CNV files, somatic mutation data, and clinicopathologic data were acquired from the publicly available TCGA database (http://xena.ucsc.edu/). Microarray dataset GSE58812 was downloaded from the GEO database (https://www.gov/geo/). A total of 1324 breast cancer samples were analyzed in this study. 1097 patients with a survival time greater than 30 days and 120 normal tissue samples were selected from the TCGA-BRCA cohort. The GSE58812 cohort contains 107 samples of breast cancer patients. The 73 Cholesterol homeostasis genes (CHGs) and 36 Angiogenesis genes (AAGs) were retrieved from the MSigDB team (Hallmark Gene set) as indicated in Table S1 .

9 CHGs were obtained with univariate Cox regression (UniCox) analysis. Consensus clustering was used to identify different cholesterol homeostasis-related patterns by the k-means algorithms with 1000 repetitions (13). The distinction in clinical characteristics between the C1 and C2 groups was assessed using a Chi-square test. Differences in the biological function of these patterns were investigated using Genetic Set Variable Analysis (GSVA) (14). OS time and OS state of various modes were compared using the Kaplan-Meier method (15). Additionally, we explored the association between molecular patterns of cholesterol homeostasis genes, clinical features, and survival differences.

The “Estimation” R package was used to present the proportion of immune cells and stromal cells in BC by analyzing gene expression, which can further calculate the tumor purity (16). Abundance of 23 specific immune cell subtypes was measured in tumors with the CIBERSORT algorithm to reveal the infiltration of immune cells (17). We predicted the sensitivity of immunotherapy by comparing the expression levels of several immune checkpoints among different subgroups. Moreover, the degree of immune cell infiltration in tumor and normal samples was determined by single sample Gene Set Enrichment Analysis (ssGSEA analysis) (18).

Using the “limma” package, we acquired DEGs for breast cancer in the TCGA dataset. DEGs should comply with the | log2 fold change (FC) | ≥ 0.5, p< 0.05. Pearson correlation analysis was used to obtain genes that were related to Cholesterol homeostasis, with |cor|≥ 0.6.

A CAG_Score was established to quantitatively evaluate the state of cholesterol homeostasis for individual BC patients. Firstly, we performed uniCox analysis and multi-factor Cox analysis (mulCox) for CHG-related genes to search for which has significant prognostic value. Then, we integrated OS time, OS, and gene expression data with the “glmnet” package and developed the CAG_Score by the Lasso Regression Algorithm (19).

The median CAG_score was adapted to classify breast cancer patients into low-risk and high-risk groups.

A CAGs-related nomograph was established to describe the clinical features and risk score of BC patients, as well as the clinical prediction of 3-year,4-year, and 5-year survival status. Calibration curves were generated to identify the accuracy of the predictive effect.

The IC50 of commonly used clinical drugs was numerically analyzed by the “pRRophic” package in order to compare the chemotherapy effects of different risk groups (20). Total RNA of breast cancer cells (MDA-MB-231, MCF-7, SKBR-3) and normal breast cells (MCF-10A) were prepared by TRIzol reagent (Thermo Fisher Scientific, Waltham, USA). cDNA was synthesized with TOROIVD qRT Master Mix kit (TOROIVD, shanghai, China) according to the manufacturer’s instructions. The qRT-PCR was performed using the TOROGreen qPCR Master Mix kit (TOROIVD, shanghai, China) on the ABI 7500 real-time fluorescence quantitative PCR system (Thermo Fisher Scientific). All sequences of primers used are shown in Table S2 .

The TCGA cohort was divided into a training cohort (n=824) and a test cohort (n=206) randomly. We defined BC prognostic risk markers as characteristic mRNA, lncRNA, and miRNA in the TCGA cohort. Screening for characteristic mRNAs, miRNAs, and lncRNAs based on high-risk score and low-risk score, for each type of data, the top 100 most relevant features were retained as BC-specific risk markers according to the P-value. Then, we performed lasso regression for further feature filtering and reduced the number of markers to 20 for each type of data. With 20 makers per molecular layer, we created a risk predictor of each single molecule layer with three machine learning models, such as Light GBM, Logistic regression, and Random forest (21). Finally, based on 60 BC-specific markers from three data types, we developed a LightGBM model (RiskLight) to distinguish breast cancer patients with different prognostic risks associated with dysregulated cholesterol homeostasis.

In the statistical analysis, p<0.05 was considered statistically significant. The t-test is used for the analysis of normally distributed data, while the Wilcoxon rank sum test is used for the analysis of abnormally distributed data. In addition, Pearson correlation analysis or Spearman analysis was used to describe the relationship between two numerical variables. The above algorithms are all implemented in R Software (version 4.1.2).

The flow chart of the research design is shown in Figure 1 . To evaluate the role of cholesterol in the occurrence of breast cancer, the Mendelian randomisation-Egger (MR-Egger) method was used first in the main MR analysis, as the detailed results are presented in Figure 2A (MR Egger, OR: 0.54, 95% CI: 0.35-0.84, p<0.006). This means that cholesterol levels may be a risk factor for breast cancer. We obtained 73 genes for cholesterol homeostasis from the MSigDB database and verified the expression levels of 73 CHGs in tumor specimens and normal control in the TCGA-BC cohort ( Figure 2B ). 63 CHGs had differential expression ( Figure 2C ). Correlations between 73 CHGs were analyzed with the String website ( Table S3 ). Protein interaction network (PPI) was constructed by Cytoscape software to explore the interactions between CHGs (22). And we identified SCD, PPARG, CTNNB1, FDPS, LDLR, ACSS2, FDFT1, FADS2, HMGCR, SREBF2, ACTG1 and HMGCS1 as the vital genes of cholesterol homeostasis ( Figure S1 ). We calculated the CNV mutation rate of CHGs, Figure 2D shows the results. In addition, we determined the incidence of SNV of 73 CHGs in BC, and 142 out of 981 BC samples (14.46%) showed mutations, which indicated that the mutation rate of 73 CHGs was less than 1% ( Figure S2 ).

Generation of a subset of genes related to cholesterol homeostasis regulation in BC to reveal the relationship between cholesterol homeostasis regulation and tumorigenesis. 1097 BC patients of TCGA-BC were included in this study, and uniCox analysis revealed 9 CHGs with prognostic significance ( Figure 3A ). To determine the relationship between CHGs expression patterns and BC subtypes, consensus cluster analysis was used to classify BC patients according to prognostic genes. When the clustering variable was 2, BC patients were well divided into the C1 group (n=510) and the C2 group (n=587) ( Figure S3 ). PCA analysis showed significant subpopulation differentiation in samples ( Figure 3B ). KM analysis revealed that cluster C2 showed a worse prognostic status ( Figure 3C ).The clinical features distinguishing the C1 and C2 groups are presented in Table S4 . In addition, the relationship between gene expression and clinical features of the two clusters was shown in Figure 3D . The heatmap indicated that the expression level of CHGs had a significant correlation with the clinical characteristics, and the genetic characteristics of the C2 subcluster were associated with distant tumor metastasis. The biological functions and signaling pathways of tumor cells were compared by the GSVA algorithm, and the findings showed that the C2 subcluster performed obvious immune pathway characteristics, lipid metabolism, and sterol metabolism-related pathways were down-regulated, and cancer metastasis-related pathways were significantly different as well ( Figure 3E ). This suggests that dysregulation of cholesterol metabolism is closely associated with tumor immunity and the development of tumors.

Investigating the infiltration extent of 23 human immune cells in both clusters by the CIBERSORT algorithm ( Figure 4A ), we found that the content of Macrophage M0, Macrophage M1, activated Dendritic cell and T cell include activated CD4 positive memory T cell, helper T cell, gamma and delta T cell were significantly higher in group C2, whereas Plasma cell, macrophage M2, resting dendritic cell mast cells behaved in an opposite way. Inter-individual differences in 23 immune cells were assessed by the “ssGSEA” algorithm and the number was generally higher in the C2 group ( Figure 4B ). The TME scores exhibited that patients in cluster C2 had a higher abundance of immune and matrix components ( Figure 4C ). In addition, PD-1, PD-L1, and CTLA-4 were shown a similar increase in cohort C2, which represents the critical expression status of the immune checkpoints (ICP) ( Figure 4D ). Meanwhile, the correlation analysis between CHGs and immune cells displayed that FBXO6, SEMA3B, GSTM2, and CXCL16 were correlated with immune cell abundance ( Figure 4E ).

The Pearson correlation algorithm was applied to analyze CHGs, resulting in 510 highly correlated DEGs ( Figure 5A ). Functional enrichment analysis of these DEGs was then performed to demonstrate the potential biological activity of cholesterol homeostasis genes. KEGG and GO analysis revealed an enrichment of cancer and metastasis-related pathways as well as blood vessel development and sterol metabolism, which suggested that cholesterol homeostasis is closely related to angiogenesis ( Figures 5B, C ). To reveal the association between cholesterol homeostasis and angiogenesis, we obtained 36 angiogenic genes (AAGs) from MsiGDB and explored the correlation between CHGs and AAGs. The results were as expected, especially when ANTXR2, GPX8, and AVPR1A were strongly associated with angiogenesis ( Figure 5D ).

Subsequently, we examined the expression of AAGs in groups C1 and C2 ( Figure 5E ), as well as in tumor and normal tissues ( Figure S4 ); however, a significant discrepancy exists. The GSVA algorithm was used to evaluate the cholesterol homeostasis score (CHG score) and angiogenesis score (AAG score) of TCGA BC samples based on 73 CHGs and 36 AAGs. And cholesterol homeostasis scores were positively correlated with angiogenesis scores in the TCGA-BC cohort ( Figure 6A ). Moreover, the cholesterol homeostasis score and angiogenesis score were compared between the C1 and C2 groups. We found that patients in the C2 group had a worse prognosis with higher cholesterol homeostasis score and angiogenesis score ( Figure S5 ). The correlation between vascular stability and cholesterol homeostasis score was also validated in addition. The abundance of genes related to vascular stability (CDH5, CLDN5, TIE1, JAM2, TEK) indicated that the group with a lower cholesterol score had higher vascular stability ( Figures 6B-F ), while low vascular stability often promotes cancer growth (23–27). All the findings were verified in the GSE58812 cohort ( Figures 6G-L ).

Considering that cholesterol homeostasis is closely connected with angiogenesis, we developed a prognostic CAG_score based on genes related to cholesterol homeostasis. The BC patients were randomly assigned to the training cohort (n=731) or the test cohort (n=366). We performed UniCox analysis of 786 cholesterol-related genes, and 49 DEGs with prognostic significance (logFC>0.5, P<0.05). Subsequently, LASSO and multi-Cox analyses were performed on 49 DEGs to establish the most suitable prediction model. We set the Lambda value to 0.00298971135072249 and finally obtained 7 genes ( Figures 7A, B ).

In the scoring model established by CAGs, we found that higher scores were associated with a worse survival rate and higher mortality rate ( Figures 7C, D ). The model genes also showed a trend with increasing CAG_score ( Figure 7E ). To evaluate the robustness of the CAG_score, we compared the CAG_score from the test to the whole cohort, and the results showed an excellent performance of the CAG_ score in assessing the prognosis of BC patients ( Figures 7F , S6 ). Figure S7 shows the distribution of CHGs and AAGs in the two CAG_score clusters. We found significant differences in gene expression in both groups.

Through the analysis of clinical indicators, we established a nomogram to predict 3, 4, and 5-year OS in BC patients ( Figure 7G ). The calibration curve shows that the method has a high forecasting accuracy ( Figure 7H ). Meanwhile, the R package “Rms” was conducted to integrate data on survival time, survival status, and 6 characteristics, and a nomogram was built using the Cox method to assess the prognostic value of these characteristics in a sample of 1030. The overall C-index of the model was: 0.783437290915762, 95%CI(0.740754644512089-0.826119937319435), p value=1.00165584281093e-38.

As mentioned above, CAG_score was positively correlated with the abundance of Macrophage M0, Macrophage M2, Plasma cell, and activated Dendritic cell, while CD8+ T cell, T cell gamma delta, activated or dormant CD4+ memory T cell, B memory cell, regulatory T cell, activated NK cell, macrophage M1 were negatively correlated with CAG_score ( Figure 8A ). In addition, there was a direct correlation between the CAG_ score and the TME score ( Figure 8B ). We explored the relevance between prognostic marker genes and 23 immune cells. We concluded that T cells and Macrophages are closely associated with the selected genes ( Figure 8C ). Furthermore, We evaluated the expression of ICPs in groups of different prognostic features. Figure 8D shows that the expression of 24 ICPs was inconsistent in both risk subgroups. The low-risk group showed a higher level of ICPs expression.

It is a meaningful research direction to select and guide the appropriate immunotherapy regimen for the patient (28). To examine the role of CAG_scores in clinical diagnosis, we evaluated the IC50 for 138 Common drugs in TCGA-BC patients. The results showed that BC patients with higher CAG_ scores were more sensitive to the AKT inhibitors VIII and Imatinib, while patients with low CAG_ scores responded better to Crizotinib, Saracatinib, Erlotinib, Dasatinib, Rapamycin, Roscovitine and Shikonin ( Figure 9A ).

We detected the RNA expression of the CAG_score genes in breast cancer cell lines. Our results indicate that all genes were highly expressed in MDA-MB-231, MCF-7 and SKBR-3 cell lines ( Figures 10A-H ), which was consistent with our prediction.

We have demonstrated the significance of metabolic regulation associated with cholesterol homeostasis as an immune micro-environmental factor in our study. To further identify molecular signatures associated with prognostic risk at the multi-omics level, we conducted an analysis of associations between three molecular layers (mRNA, miRNA, lncRNA) and high-low risk for each type of data, the top 100 most relevant features were retained as BC-specific risk markers according to the P-value ( Figure S8A ). We used the Light GBM framework to integrate multi-omics features to develop high- and low-risk prediction models as a way to emulate the tumor micro-environment in which cholesterol homeostasis is dysregulated. As a result, the three risk predictors based on the single molecular layer performed well in predicting high and low risk in the test cohort (AUC=0.8491 for the mRNA model, AUC=0.7939 for the lncRNA model, and AUC=0.7970 for the miRNA model) ( Figures 11A, C, E ). We also compared the superiority of random forest and logistic regression models with the Lightgbm model ( Figures 11B, D, F ). The results show that all three algorithms exhibit consistent results. Finally, we integrated 3 risk predictors, based on the LightGBM algorithm combined with multi-omics data to develop an integrated model (Risklight) for predicting cholesterol homeostasis-related risk patterns. Risklight is superior to all risk predictors based on single molecular layers (AUC=0.89) ( Figures 11G, H ).