A cholesterogenic gene signature for predicting the prognosis of young breast cancer patients

TCGA (BRCA-US) and GEO (GSE131769) data collection

BC transcriptome data (HTSeq-FPKM) and clinical materials were downloaded from the TCGA (https://cancergenome.nih.gov) and GEO database (https://www.ncbi.nlm.nih.gov/geo). The BC patients from the TCGA were from the United State National Cancer Institute. While the patients from GSE131769 (GEO database) were collected from Yonsei University College of Medicine, Korea. The samples from these two datasets did not overlap. The selection criterion was age ≤45 years. A total of 183 tumour samples and 30 normal samples in the TCGA were included as the training set. The GSE131769 dataset, with data from the GPL6947 Illumina HumanHT−12 V3.0 expression bead chip, was used for clinical external validation. After excluding patients older than 45 years, 130 patients with expression array data and survival data were enrolled.

Cholesterogenic gene preparation and differentially expressed genes (DEGs)

Cholesterogenic genes were downloaded from the Molecular Signatures Database (MSigDB) (Liberzon et al., 2011). Cholesterogenic gene sets were considered if the description contained “cholesterol”. A total of 30 cholesterol-related gene sets were selected (Table S1).

The transcriptional expression of these cholesterogenic genes was incorporated into the TCGA samples using Perl (Perl is a highly capable, feature-rich programming language. The Perl package was uploaded to the sorce code). Then, we obtained the DEGs by comparing the expression of the RNA sequences of 183 tumour samples with those of 30 normal samples using R language, where the “limma” and “ggplot2” packages were applied for statistical analysis (“limma” package was used for differential analysis; “ggplot2” package was used for drawing volcano plot). In order to include all the possible genes in the preliminary screening, we just use P < 0.1 as the cut-off value.

Functional and pathway enrichment analysis

Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.8 (https://david.ncifcrf.gov/summary.jsp) was used to analyse the Gene Ontology (GO) functional enrichment and Kyoto Encyclopedia (KEGG) pathway enrichment of the DEGs (Jiao et al., 2012). The description of GO functional enrichment consisted of biological processes (BPs), molecular functions (MFs), and cellular components (CCs). R language was employed for visualization using the “ClusterProfiler”, “enrichplot” and “ggplot2” packages. The cut-off values were set at P < 0.05 and false discovery rate (FDR) < 0.05.

Construction of a five-gene signature in the TCGA dataset

First, we used the “survival” package in R to carry out univariate Cox regression to identify prognosis-related cholesterogenic genes. Then, by combining survival status and survival time of the patients, multi COX regression was performed to further screen out survival-related genes with corresponding hazard ratio value and correlation coefficient. “Survival” package in R was applied. Third, the genes closely correlated with survival were used to construct a risk scoring model. After filtering out 4 patients with incomplete follow-up data, we stratified 179 patients into high-risk and low-risk groups according to the median risk score and compared the overall survival between the two groups with Kaplan–Meier analysis using the “survminer” package. Fourth, univariate and multivariate analyses were implemented to assess the prognostic value of the model based on age, clinical stage, surgical method, and risk score. Fifth, the receiver operating characteristic (ROC) curve was drawn to evaluate the predictive value of the signature using the “survival ROC” package.

Subgroup analysis in TCGA and external validation in GSE131769

Subgroup analysis was carried out according to the clinical features of the patients in the TCGA database. The patients were divided into two groups by the dichotomy method according to age (≤40 years and >40 years), surgical method (breast-conserving surgery (BCS) and mastectomy), T stage (T1-2 and T3-4) and N stage (N0-1 and N2-3) to verify the predictive ability of the signature.

BC patients aged 45 years or younger in GSE131769 with microarray data and corresponding clinical information were used for further external validation.

Expression of the five screened cholesterogenic genes in normal tissues and tumour tissues

The expression of GRAMD1C, NFKBIA, INHBA, CD24, and ACSS2 in the TCGA training set was subjected to log2 transformation. We compared the expression of the above five cholesterogenic genes between normal tissues and tumour tissues. The P value was set at 0.05.

Cell culture and transient transfection

The human BC cell lines MDA-MB-231 and MCF7 were obtained from the American Type Culture Collection (Manassas, VA, USA). The cells were cultured in Dulbecco’s modified Eagle medium (DMEM, Gibco, USA) with 10% foetal bovine serum (FBS, Gibco, USA) and incubated at 37 °C in a humidified atmosphere with 5% CO₂.

Cells were digested with trypsin and seeded in 24-well plates at a density of 7 × 10⁴ cells/500 µL per well. Transient transfection of the plasmids was performed using Lipofectamine 3000 reagent (Invitrogen, Waltham, MA, USA) according to the manufacturer’s instructions.

Cell counting kit-8 (CCK-8) analysis, colony formation assay and transwell migration assay

For CCK-8 analysis, the cells (2 ×10³) were seeded into 96-well plates. At the indicated time point, 100 µL of CCK-8 dye solution was added to each well and incubated for 1 h, and the absorbance was measured at a wavelength of 450 nm. Cell viability was calculated by the following formula: cell viability = (optical density (OD) value of the experimental group/OD value of the control group) × 100%. For the colony formation assay, the indicated cells (0. 2 ×10³) were seeded into six-well plates and maintained at 37 °C. After 10 days, colonies were fixed with 10% formaldehyde for 15 min and stained with 1.0% crystal violet for 30 s. Colony morphologies were captured by a light microscope (Olympus).

Matrigel (BD Biosciences, Bedford, MA, USA) was prebalanced at 4 °C and mixed with 8 times the volume of precooled blank culture medium. An 80 µL Matrigel-medium mixture was coated on the upper polycarbonate 24-well culture chamber. Cells (1 ×10⁴) were plated on the top side coated with Matrigel and incubated at 37 °C for 22 h, followed by removal of the cells inside the upper chamber using cotton swabs. Invaded cells remaining on the lower membrane surface were fixed in 1% paraformaldehyde, stained with crystal violet, and counted (ten random fields per well, 100 × magnification).

Differentially expressed cholesterogenic genes and enrichment analysis

We downloaded TCGA-BRCA RNA-Seq data (HTSeq-FPKM) and clinical features from the TCGA database, which contained 1098 tumour samples and 123 normal samples. We selected patients aged 45 years or younger. Finally, 183 tumours and 30 normal tissues were included.

Then, cholesterol-related genes were extracted from MSigDB, with a total of 345 cholesterogenic genes (Table S2). We compared the expression of cholesterogenic genes in tumour tissues and normal tissues to obtain the DEGs. Considering only primary filtering, P = 0.1 was set as the criterion. A total of 97 upregulated genes and 124 downregulated genes were screened out (Fig. S1).

To explore the underlying mechanisms of the DEGs, we conducted GO functional and KEGG pathway enrichment analyses by using DAVID 6.8. In BP analysis, the lipid localization ranked the first. The membrane raf, membrane microdomain and membrane region were enriched in CC analysis. MF analysis showed that the DEGs were mainly enriched in lipoprotein particle binding and protein-lipid complex binding, which play equally important roles. Moreover, cholesterol metabolism was the most significant pathway in the KEGG analysis (Fig. 1).

Construction of the cholesterogenic gene signature in the TCGA-BRCA dataset and survival analysis

Univariate Cox regression was conducted to filter the survival-related cholesterogenic genes. GRAMD1C, NFKBIA, INHBA, CD24, MYLIP, ACSS2, and NR1H3 were screened out with the criteria of P < 0.05. Multi COX regression was performed to further select prognosis-related genes (GRAMD1C, NFKBIA, INHBA, CD24, and ACSS2). We constructed a five-gene signature in the TCGA cohort, with risk score formula = ∑ⁱ (coefi*expri), (with i meant the number of genes, coefi meant the coefficiency of the genes, and expri meant the genes in the signature). As a result, Risk scores were calculated according to the following formula: risk score = − 1.169 × GRAMD1C − 0.992 × NFKBIA + 0.432 × INHBA + 0.261 × CD24 − 0.839 × ACSS2. Four patients with follow up <3 months were excluded, the remaining 179 patients were divided into the high-risk group (n = 89) and the low-risk group (n = 90) according to the median risk score. Kaplan–Meier survival analysis showed that the low-risk group had a significantly better prognosis than the high-risk group (P < 0.001) (Figs. 2A, 2C, 2D, 2E). The ROC curve showed that the area under the curve (AUC) value was 0.810 in TCGA dataset, which implied that the five-gene signature had excellent power in predicting the survival of young BC patients (Fig. 2B).

We combined age, clinical stage, surgical method (BCS vs. mastectomy), T stage, N stage and risk score in univariate and multivariate factor analyses. We found that the risk score was an independent risk factor for prognosis in young BC patients (Fig. 3).

Subgroup analysis in TCGA and external validation in GSE131769

We performed subgroup analysis by stratifying the young BC patients in TCGA-BRCA according to clinical features (age, T stage, N stage, and surgical method). The five-cholesterogenic-gene signature could be applied to differentiate prognosis in different subgroups of patients, including the age ≤40 years and age>40 years subgroups, the T1-2 subgroup, the N0-1 and N2-3 subgroups, and the mastectomy subgroup (P < 0.05) (Figs. 4A, 4B, 4C, 4E, 4G, 4H). We noticed a different survival trend in the BCS and T3-4 subgroups, although the P value was above 0.05 (Figs. 4D, 4F). This result is likely because of the small sample sizes of these subgroups.

To further confirm the model, we performed external validation by using GSE131769, which contains 130 patients aged 45 years or younger. By calculating the risk score for each person, we found that the cholesterogenic gene signature also worked well in the external validation set. The low-risk group (n = 65) had a much more favourable prognosis than the high-risk group (n = 65, P = 0.033) (Fig. 5).

Expression of the five cholesterogenic genes in normal tissues and tumour tissues

The expression of CD24 and INHBA was significantly high in tumour tissue (P = 0.0102 and P < 0.001). The expression of GRAMD1C, NFKBIA, and ACSS2 was significantly higher in normal tissues than in tumour tissues (P < 0.001, P = 0.0096, and P < 0.001, respectively), which indicated that CD24 and INHBA were risk genes and that GRAMD1C, NFKBIA, and ACSS2 were protective genes (Fig. S2).