Expression variations of lipid metabolism-associated genes in breast cancer
RNA-seq from 1100 breast cancer tissue samples and 112 non-tumor samples were downloaded from The Cancer Genome Atlas (TCGA). Given the criteria for FDR < 0.05 and | log2 (fold change) |> 2, we finally identified 80 up-regulated and 105 down-regulated lipid metabolism-associated genes (LMAGs) by “DESeq2” R package. The expression levels of these genes were shown in the heatmap (Fig. 1A), where blue represented the down-regulation of gene expression and red represented the up-regulation of gene expression. Additionally, the volcano map (Fig. 1B) was labeled with the names of differentially expressed LMAGs.
Development of a prognostic signature in the training cohort
A prognostic signature was constructed to apprise the prognostic value of the differentially expressed LMAGs in breast cancer. All enrollers were randomly divided into training and validation cohorts in a ratio of 6:4. We developed the model in the training cohort with 651 samples. Univariate Cox analysis was initially performed under the benchmark for P-value < 0.005. According to the forest map of the hazard ratio (Fig. 2A), three differentially expressed LMAGs (ALDH2, CYP21A2, and IL24) were screened out for subsequent analysis. Lasso regression revealed that these 3 genes were preserved when the tuning parameter (λ) was minimized to 0.0027 and log(λ) was − 5.93 (Fig. 2B).
After the preceding univariate and lasso analysis, ALDH2, CYP21A2, and IL24 were eventually included in the prognostic signature by multivariate cox regression. RiskScore was calculated as follows: RiskScore = h0(t) × exp [(-2.231 × expression level of ALDH2) + (-1.363 × expression level of CYP21A2) + (-2.353 × expression level of IL24)]. Based on the median RiskScore, patients in the training cohort were classified into low-risk (< median number) and high-risk (> median number) groups. The overall survival rates of patients in the high- and the low-risk group were calculated by Kaplan-Meier, indicating that higher RiskScore corresponded with poorer survival (P < 0.001, Fig. 2C). In the light of model diagnosis, the areas under the curve (AUC) of time-dependent ROC curves were 0.7, 0.7, and 0.65 for 1-, 3-, and 5-year survival, respectively (Fig. 2D). Patients in the high-risk group had high mortality and a shorter overall survival (OS) time than those in the low-risk group (Fig. 2E). A heatmap depicting the expression pattern of the genes incorporated in the risk model suggested that patients in the high-risk group tended to express a lower level of risk genes (Fig. 2F).
Internal validation of the risk model
Of the remaining breast cancer patients in the TCGA cohort, a total of 431 samples were employed as the validation set. According to the median RiskScore in the training cohort, 231 patients were encompassed in the low-risk group, while 200 patients were encompassed in the high-risk group. Kaplan-Meier analysis uncovered a significant difference in survival between the high- and low-risk groups (P < 0.01, Fig. 3A), with patients in a low-risk group having better survival outcomes. ROC curves verified that our risk model had fine predictive potency (AUC = 0.81 for 1-year, 0.7 for 3-year, and 0.64 for 5-year survival, Fig. 3B). Additionally, RiskScore, survival status distributions, and the risk genes expression pattern of the patients showed the same characteristics as the training cohort (Fig. 3C-D).
Clinical application value of the prognostic signature
Considering the significance of RiskScore in portending the prognosis of breast cancer patients, we further explored its value in clinical application. Univariable and multivariable Cox regression were performed to assess the potential independent prognostic role of RiskScore. Univariable Cox regression demonstrated that RiskScore was an independent predictor of poor survival in both the training and validation cohorts (HR = 1.534, 95%CI: 1.199–1.961 and HR: 1.510, 95%CI: 1.080–2.110, Fig. 4A). After adapting for other confounders, multivariable Cox regression also revealed that RiskScore (HR = 1.683, 95%CI: 1.308–2.167 and HR: 1.765, 95%CI: 1.232–2.527, Fig. 4B) played an independent role in predicting prognosis of breast cancer patients. A nomogram model was constructed based on the age, TNM stage, tumor stage, and RiskScore. The survival contingency of breast cancer patients for 1-, 3-, and 5 years were computed according to the total points (Fig. 4C). A C-index of 0.781 suggested that the nomogram had an agreeable predictive value. Calibration plot for the probability of survival also showed satisfactory compliance with the prediction of 1-, 3-, and 5-year OS (Fig. 4C). These results indicate that the risk model has a weighty hand in clinical prediction.
Immune microenvironment landscape between high- and low-risk group
The discrepancy in the immune microenvironment between high- and low-risk groups was evaluated by multiple immune analyses. Outcomes of the ESTIMATE algorithm demonstrated that patients in the high-risk group had lower Immune Score and ESTIMATE Score, as well as a higher level of Tumor purity (Fig. 5A). The abundance of immune cell infiltration in high- and low-risk groups predicted by the CIBERSORT analytical tool manifested that the high-risk group possessed a significantly lower abundance of immune cells, including naïve B cells, CD8 + T cells, naïve CD4 + T cells, resting memory CD4 + T cells, activated memory CD4 + T cells and monocytes. While the infiltration of T cells follicular helper, resting NK cells, macrophages (M0), and resting mast cells was prominently up-regulated (Fig. 5B). In addition, we investigated the expression of immune checkpoints in high- and low-risk groups, finding that CTLA4, PD1, LAG3, BTLA, and TNFRSF14 were significantly down-regulated in the high-risk group (Fig. 5C).
Weighted gene co-expression network analysis singled out the module closely related to immune infiltration
To figure out the predominant module that correlated with clinical traits or immune infiltration, we performed the weighted gene co-expression network analysis (WGCNA) on the TCGA cohort. After the soft threshold power β being set at 6 (scale-free R2 = 0.90, Fig. 6A), a total of 14 modules with unique colors were identified, each comprising genes with similar expression patterns recognized by hierarchical clustering (Fig. 6B). And genes contained in the gray module did not belong to any other modules. Module-trait relationships analysis was carried out to discern the association between co-expression modules and immune infiltration, as well as clinical traits (Fig. 6C). In accordance with the correlation heatmap, the MEbrown module had a strikingly positive correlation with ESTIMATE score (R = 0.81, P = 1e-211). Nevertheless, none of the modules had a strong link with clinical traits. Intriguingly, one of these genes in the risk model, ALDH2, was also located in the MEbrown module.
Potential prognostic value of ALDH2 in breast cancer
We explored the expression variants of ALDH2 in different TNM stages and tumor stages, noticing that ALDH2 was up-regulated in samples with lymphatic metastasis (P < 0.01, Fig. 7A). However, ALDH2 expression manifested no divergence in samples with T3-T4, M1 and Stage III-IV when compared with samples with T1-2, M0 and Stage I-II, respectively. Following the median expression level of ALDH2, patients were classified into two subgroups. Comparisons of the patients’ OS between the two subgroups were made by Kaplan-Meier curves, revealing that patients with high expression of ALDH2 had long-term OS (Fig. 7B).
Downregulation of ALDH2 in breast cancer impacted tumor malignancy
Western blotting was conducted to assess the knockdown efficiency after the transfection of siRNA targeting ALDH2 in MCF7 and MDA-MB-231 cells (Fig. 8A). Wound-healing assay exhibited a significant high migratory rate of ALDH2-silenced MDA-MB-231 cells, but not of MCF7 cells (Fig. 8B). We further performed cell invasion assay to compare the invasion ability of breast cancer cells between NC and si-ALDH2 groups. The quantity of cells invading to the lower chambers were more in MCF7 cells and MDA-MB-231 cells of si-ALDH2 group than that of NC group (Fig. 8C). These results indicate that downregulation of ALDH2 promotes the invasion and migration capacity of breast cancer cells.
Next, invasive breast ductal carcinoma tissue from 30 patients pasted on tissue microarray was analyzed by immunohistochemistry (IHC). Of the 30 patients, 16 patients were diagnosed with Nottingham grade II, and the 6 patients with grade III. The expression level of ALDH2 was compared between 16 patients with grade II and 6 patients with grade III, which showed no significant difference between the two groups (Supplementary Figure S1). The limited number of clinical samples may be one of main reasons for no differential expression of ALDH2 in breast cancer tissues with various pathological grades.
Potentially underlying mechanism of ALDH2 for modulating lipid metabolism
To investigate the predicted underlying mechanism of ALDH2 for modulating lipid metabolism, we performed the correlation analysis of ALDH2 in TCGA cohort. Under the screening condition of Pearson’s r > 0.2, 62 ALDH2-related genes were further probed by GO and KEGG enrichment analysis (Table 1). GO enrichment analysis showed that ALDH2-related genes involved in regulation of lipid metabolic process, lipid transport, alcohol metabolic process, lipid catabolic process and fatty acid metabolic process (Fig. 9A). KEGG pathways of ALDH2-related genes were mainly related to PPAR signaling pathway and AMPK signaling pathway (Fig. 9B). Most of the ALDH2-related genes enriched in PPAR signaling pathway, AMPK signaling pathway and fatty acid degradation were remarkably down-regulated (Fig. 9C). Network diagrams of the above three pathways showed the genes contained in each (Fig. 9D). Among them, ALDH2 was enriched in fatty acid degradation pathway. We also inspected the role of ALDH2-related genes played in PPAR signaling pathway, AMPK signaling pathway and fatty acid degradation pathway. Scatter plot demonstrated the correlation between ALDH2 and the genes were enriched in these pathways, whose Pearson value was greater than 0.3 (Fig. 9E).
Table 1
The co-expressed genes of ALDH2 (Pearson’s r > 0.4)
Gene | Pearson correlation coefficient | P-value |
G0S2 | 0.393086 | 5.83E-42 |
CEBPA | 0.387443 | 1.03E-40 |
GPD1 | 0.371073 | 3.07E-37 |
CES1 | 0.370515 | 4.00E-37 |
PPARG | 0.370499 | 4.03E-37 |
PLIN1 | 0.369284 | 7.17E-37 |
RBP4 | 0.364016 | 8.43E-36 |
CIDEC | 0.356032 | 3.24E-34 |
CIDEA | 0.353297 | 1.10E-33 |
HSD11B1 | 0.346122 | 2.60E-32 |
GPX3 | 0.341698 | 1.75E-31 |
CCL21 | 0.338768 | 6.08E-31 |
FABP4 | 0.334929 | 3.05E-30 |
ADH1B | 0.334714 | 3.34E-30 |
GPBAR1 | 0.328521 | 4.28E-29 |
FUT7 | 0.327943 | 5.42E-29 |
THRSP | 0.325781 | 1.30E-28 |
LIPE | 0.321347 | 7.66E-28 |
GPIHBP1 | 0.318527 | 2.33E-27 |
ADIPOQ | 0.316617 | 4.91E-27 |
CDO1 | 0.316036 | 6.16E-27 |
ADH1C | 0.31455 | 1.10E-26 |
PCK1 | 0.31381 | 1.46E-26 |
LEP | 0.311488 | 3.55E-26 |
ALDH1A2 | 0.306692 | 2.18E-25 |
CD19 | 0.291018 | 6.52E-23 |
ABCD2 | 0.290849 | 6.92E-23 |
FOXP3 | 0.286725 | 2.92E-22 |
ABCA9 | 0.285313 | 4.76E-22 |
HSD17B13 | 0.285249 | 4.86E-22 |
ABCA6 | 0.280071 | 2.84E-21 |
DPEP1 | 0.275443 | 1.33E-20 |
COL1A1 | 0.274019 | 2.13E-20 |
CXCL9 | 0.272996 | 2.98E-20 |
ABCA8 | 0.271679 | 4.58E-20 |
ABCB1 | 0.270824 | 6.05E-20 |
IL12B | 0.261906 | 1.04E-18 |
CRHBP | 0.255216 | 8.13E-18 |
AKR1C1 | 0.253674 | 1.30E-17 |
ZFP36 | 0.251817 | 2.26E-17 |
BMX | 0.242248 | 3.73E-16 |
CD36 | 0.239409 | 8.37E-16 |
P2RY12 | 0.237706 | 1.35E-15 |
CAV1 | 0.235943 | 2.22E-15 |
F7 | 0.235352 | 2.61E-15 |
LPL | 0.232621 | 5.55E-15 |
ENPP2 | 0.231905 | 6.75E-15 |
TACR1 | 0.226428 | 2.96E-14 |
ABCA10 | 0.221977 | 9.56E-14 |
KL | 0.221878 | 9.81E-14 |
AKR1C2 | 0.21886 | 2.14E-13 |
TIMP4 | 0.217554 | 2.99E-13 |
KLF4 | 0.214531 | 6.43E-13 |
FOS | 0.214395 | 6.66E-13 |
ACADL | 0.211172 | 1.49E-12 |