Prognostic analyses of genes associated with anoikis in breast cancer

Data collection

RNA sequencing data of patients with breast cancer were downloaded from TCGA database (https://portal.gdc.cancer.gov/) using the R package TCGAbiolinks (Colaprico et al., 2016). After eliminating samples with incomplete clinical annotations, 1,083 cases of breast cancer projects, including 111 cases with matched adjacent tissues, were obtained (TCGA-BRCA), and the data format was level 3 HTSeq-fragments per kilobase per million. Additionally, corresponding clinical information was obtained from the UCSC Xena database (Goldman et al., 2020). The mRNA counts of TCGA-BRCA dataset were normalised by using the limma R package (Ritchie et al., 2015).

The breast cancer-related datasets GSE42568 (Clarke et al., 2013), GSE102484 (Cheng et al., 2017), and GSE20685 (Kao et al., 2011) were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2013) using the R package GEO query (Davis & Meltzer, 2007). The detection platform for the three datasets was GPL570 (HG-U133_Plus_2). GSE42568 included 104 breast cancer samples and 17 matching normal tissues, GSE102484 comprised 683 breast cancer samples, and GSE20685 comprised 327 breast cancer samples.

Eight gene sets of ARGs were retrieved from the Molecular Signatures Database (MSigDB) (Liberzon et al., 2015), which provides the most comprehensive information on human gene sets, resulting in 862 ARGs. Furthermore, 794 ARGs were retrieved from the GeneCards database (Stelzer et al., 2016). By taking the intersection of the two databases and omitting those genes with missing probes, we obtained 100 ARGs that are common genes in both databases (Table S1).

Prognostic model construction

The least absolute shrinkage and selection operator (LASSO) prognostic model of 100 ARGs for breast cancer was constructed using the TCGA-BRCA dataset, employing 10-fold cross-validation with seeds 2021 and a p-value of 0.05. We ran 1000 iterations to guard against overfitting. The outcomes of the LASSO analysis were visualised. All samples from the TCGA-BRCA dataset and the BRCA dataset were divided into high- and low-risk groups, respectively, according to the median risk score of the prognostic model, and the survival outcomes were displayed in spot graphs. Moreover, the ARG expression profiles for prognostic model construction were visualised using heatmaps.

Differential expression genes in the high- and low-risk groups

First, the batch effects of GSE42568, GSE20685, and GSE102484 were eliminated using the R package sva (Leek et al., 2012), and merged BRCA datasets were retrieved from the GEO database. Second, datasets from TCGA and GEO were standardised and then grouped into high- and low-risk groups according to the median risk score of the LASSO prognostic model. Differential gene expression analysis was performed between the two groups, and a volcano plot (ggplot2) and heatmap (R package pheatmap) were created to visualise up-regulated genes (log FC > 0.5 and p < 0.05) and down-regulated genes (log FC < −0.5 and p < 0.05).

Gene enrichment analysis and the PPI network

The clusterProfiler R package was used to perform Gene Ontology (GO) (Gene Ontology Consortium, 2015) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) (Kanehisa & Goto, 2000) pathway enrichment analyses (Yu et al., 2012) with p < 0.05 and a false discovery rate (FDR) value ( q-value) of < 0.2. The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database was accessed to provide information regarding known and predicted PPIs (Szklarczyk et al., 2019), From which information regarding ARGs was extracted from the STRING database, and a PPI network was constructed using Cytoscape (Shannon et al., 2003). Additionally, GeneMANIA was used for gene prediction of the selected ARGs, and the results were visualised (Franz et al., 2018).

Gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA)

GSEA is an analytical method that determines whether predefined gene sets demonstrate statistically significant and consistent differences between two biological states (Subramanian et al., 2005). In the current study, TCGA-BRCA genes were divided into phenotype-related high- and low-risk groups based on the risk score median of the LASSO prognostic model, and GSEA was performed using the clusterProfiler R package. The procedure was repeated 1,000 times for each analysis, and c2.cp.v7.2.symbols were obtained from MSigDB to serve as the reference gene collection. Statistical significance was set at FDR < 0.25 and p < 0.05.

We downloaded “h.all.v7.4.symbols.gmt” from the MSigDB database to perform GSVA based on gene expression levels in the TCGA-BRCA dataset and analyse variations of functional enrichment between the high- and low-risk groups. Statistical significance was set at p < 0.05.

Analysis of the prognostic value and receiver operating characteristic (ROC) curve for ARGs

Kaplan–Meier curves were used to estimate the overall survival of patients with breast cancer based on the ARG expression derived from the TCGA-BRCA dataset. The ROC curve was used to assess the relationship between ARGs and tumourigenesis, and genes with an area under the curve (AUC) of the ROC curve >0.6 were considered key genes.

Correlation analysis between clinical characteristics and prognosis

First, we performed univariate Cox regression to evaluate the prognostic value of clinical characteristics and key ARGs in patients, and a multivariate Cox regression model was constructed from factors with a p-value of <0.1 in univariate Cox regression. Subsequently, nomograms were used to predict the 2-, 3-, and 4-year survival probabilities individually based on a multivariate Cox regression model. The calibration curves were evaluated graphically by plotting the nomogram-predicted probabilities against the observed occurrences, and the 45° line represented the best predictive values. The nomograms and calibration curves were produced by the R package rms. Decision curve analysis (DCA) was performed using the R package ggDCA (Tataranni & Piccoli, 2019) to evaluate the predictive effect of the nomograms for survival probability.

Cell cultures and real-time qPCR

The human normal breast epithelial cell line MCF-10A and breast cancer cell line HCC1954, MDA-MB-231, and MCF-7 were procured from the American Type Culture Collection. MDA-MB-231 is oestrogen (ER)/progesterone (PR)/HER2 negative, HCC1954 is ER/PR negative and HER2 positive, and MCF-7 is ER/PR positive and HER2 negative. MCF-10A cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM)/F12 medium (GIBCO, USA) with 5% horse serum (HyClone, USA) and other supplements. MDA-MB-231and MCF-7 cells were cultured in DMEM (GIBCO, USA) containing 10% foetal bovine serum (Gibco, Billings, MT, USA), and HCC1954 cells were cultured in Rosewell Park Memorial Institute 1640 medium (GIBCO, USA) containing 10% foetal bovine serum (Gibco). All cell lines were incubated at 37 °C under 5% carbon dioxide according to the recommended protocols.

Total RNA was extracted from cells which had been cultured in dishes using TRIzol (Life Technologies, USA). Reverse transcription was then performed using an Evo M-MLV RT Mix Tracking Kit with genomic deoxyribonucleic acid (DNA) Clean (Accurate Biology, Changsha, China), and qPCR was performed using an HS SYBR Green Premix Pro Taq qPCR Tracking Kit (Accurate Biology, Changsha, China). All steps were performed according to the manufacturer’s instructions. The primer sequences were as follows:

The ΔΔCt method was used to analyse qPCR data. Normalised ΔCt values were obtained by subtracting the Ct values of the housekeeping reference gene (β-actin) from the target gene Ct values. ΔΔCt values were calculated by subtracting the testing group target gene ΔCt values from the average ΔCt of the control group. Finally, relative mRNA expression was obtained by 2^−ΔΔCt.

Immunohistochemical (IHC) analysis

IHC information on key genes in breast cancer and normal breast tissues was obtained from the Human Protein Atlas (HPA) database (Colwill & Gräslund, 2011).

Statistical analysis

Clinical data handling in the study was performed using the statistical software R (version 4.1.2; R Core Team, 2021). Continuous variables were summarised using the mean ± standard deviation. The Wilcoxon-rank sum test was used to compare the two groups, and Student’s t-test was used to compare the normally distributed data of the two groups. The Kruskal–Wallis test was used for comparing more than two groups, and categorical variables were compared using the chi-square or Fisher’s exact tests. Survival analyses and univariate and multivariate Cox regressions were performed using the survival package in R, and LASSO analyses were performed using the R package glmnet. Furthermore, survivalROC in the R package was used for creating the ROC curve. The significance of survival difference was tested with the log-rank test. In addition, two-sided Spearman’s rank correlation coefficients were employed. One-way analysis of variance was used for comparing multiple samples of experimental data. Differences with p-values <0.05 were considered significant.

LASSO prognostic models were established and 11 ARGs with possible prognostic values were screened

A flow chart of the study is provided in the Fig. 1. BRCA datasets of GEO were obtained from GSE42568, GSE20685, and GSE102484 after eliminating batch effects. Box plots and principal component analysis plots displayed the datasets before and after eliminating batch effects (Figs. S1A–S1D). The results revealed that the batch effect had been effectively eliminated. We performed LASSO analysis of ARGs based on the TCGA-BRCA dataset (Fig. 2A) to evaluate the prognostic value of 100 ARGs from the intersection of MSigDB and the GeneCards database (Table S1), and the variable trajectories graph was used to visualise the results (Fig. 2B). The risk factor graph demonstrated that those in low-risk groups had longer survival compared with those in high-risk groups (Fig. 2C) and 12 ARGs (BIRC3, CD24, FGFR1, IVL, KRT15, L1CAM, MIA, NDRG1, NOS2, PTPN3, SERPINA1, and TP63) of potential prognostic value were screened (Fig. 2C) (Table 1). Next, we explored the 12 ARG profiles between two risk groups in the TCGA-BRCA and BRCA datasets, respectively, and visualised the results with pheatmaps (Figs. 2F and 2G). The box plot further showed a comparison of each ARG between the two risk groups, and there were 11 significant differentially expressed ARGs (BIRC3, CD24, IVL, KRT15, L1CAM, MIA, NDRG1, NOS2, PTPN3, SERPINA1, and TP63) with the same expression trends in both datasets (Figs. 2H and 2I).

Subsequently, we found 878 differentially expressed genes between the two risk groups based on the TCGA-BRCA dataset following thresholds of |log FC| > 0.5 and P.adj < 0.05. Of these, 244 genes were up-regulated, and 634 were down-regulated in the high-risk group, and the volcano plot displayed the results (Fig. 2D). In addition, 196 genes were obtained from the BRCA dataset under the same thresholds comprising 57 up-regulated genes and 139 down-regulated genes (Fig. 2E).

Results of the gene enrichment analysis and PPI network based on screened ARGs

We performed GO and KEGG enrichment analysis, including cellular components, biological processes (BPs), molecular functions (MFs), and signalling pathways, with the criteria of p < 0.05 and FDR (q-value) < 0.25 (Table 2) to elucidate the biological functions of the selected 11 ARGs. GO analysis demonstrated that 11 ARGs are predominantly enriched in BP-related processes (such as keratinocyte differentiation, response to hypoxia, and negative proteolysis regulation) and MF-related processes (such as endopeptidase inhibitor activity, peptidase inhibitor activity, and flavin mononucleotide [FMN] binding) (Fig. S2A). KEGG pathway analysis showed that the 11 ARGs were primarily associated with small-cell lung cancer and toxoplasmosis (Fig. S2A). Figure S2B illustrates the potential links between these items. Furthermore, the chord diagram illustrated each ARG’s distribution in the GO items and KEGG pathways (Fig. S2C). CD24, NOS2, and NDRG1 were associated with the process of response to hypoxia. KRT15, TP63, and IVL were mainly enriched in the process of keratinocyte differentiation. BIRC3 and NOS2 were closely associated with small-cell lung cancer pathways. SERPINA1 and BIRC3 were associated with peptidase inhibitor activity.

PPI networks were constructed using the STRING database with an interaction threshold of 0.15 to ascertain the interactions between the 11ARGs in breast cancer. Visualisations were performed using the Cytoscape software (Fig. S2D). The results indicated potential interaction between ARGs. Furthermore, genes with similar functions to ARGs were predicted, and a protein network was constructed using GeneMANIA, with ephrin type-B receptor 2, RAN binding protein 9, and activated leukocyte cell adhesion molecule being closely associated with the 11 ARGs (Fig. S2E).

Signalling pathways and hallmark gene sets associated with risk groupings

We subsequently performed GSEA and GSVA between high- and low-risk groups based on the TCGA-BRCA dataset to explore potential relationships between ARGs and gene sets related to breast cancer. The results of GSEA indicated several significantly enriched BPs (Table 3), including pre-Notch expression and processing (Fig. 3B), signalling by Notch (Fig. 3C), T cell factor-dependent signalling in response to Wnt (Fig. 3D), and signalling by Wnt (Fig. 3E), as illustrated in Fig. 3A. Other BPs, including transcription of androgen receptor-regulated genes kallikrein-related peptidase (KLK) 2 and KLK3, DNA damage telomere stress-induced senescence, and processes associated with cell cycle regulation, are presented in Table 3. We performed GSVA on all genes in TCGA-BRCA to further explore the differentially expressed hallmark gene sets between the two risk groups. Table 4 reveals that 38 hallmarks exhibited differences, which were visualised using heat maps (Fig. 3F), and the gene sets with p-values <0.001 were compared between two risk groups and visualised using box plot diagrams (Fig. 3G). The results implied that gene sets of Myc targets, E2F targets, hypoxia, p53, tumour necrosis factor (TNF)-α/nuclear factor kappa B (NF- κB), interleukin (IL)-6/Janus kinase (JAK)/signal transducer and activator of transcription 3 (STAT3), and mammalian target of rapamycin complex 1 (mTORC1) signalling pathway had marked differences in their enrichment between the two groups.

Analysis of the prognostic significance and ROC curve of the ARGs

Kaplan–Meier survival curves were plotted for 11 ARGs based on the TCGA-BRCA dataset (Figs. 4A–4I) to further explore the roles of ARGs in the prognosis of patients with breast cancer. We found that nine ARGs impacted the prognosis of patients with breast cancer (p < 0.05). Specifically, BIRC3 (hazard ratio [HR] = 0.67, p = 0.015), KRT15 (HR = 0.61, p = 0.002), MIA (HR = 0.68, p = 0.021), SERPINA1 (HR = 0.58, p = 0.001), and TP63 (HR = 0.64, P = 0.007) were positively correlated with overall survival, while CD24 (HR = 1.69, p = 0.002), IVL (HR = 1.71, p = 0.001), NDRG1 (HR = 1.41, p = 0.036), and NOS2 (HR = 1.46, p = 0.022) were associated with a poor prognosis. We also performed ROC curves to evaluate the correlation between ARGs and tumourigenesis using TCGA-BRCA dataset. As a result, the AUC of the ROC curve was calculated, and five key ARGs (CD24, KRT15, MIA, NDRG1, and TP63) with AUC > 0.6 were demonstrated (Figs. 4J–4N).

Correlation analysis between the clinical characteristics and prognosis

Clinical information of patients from TCGA-BRCA dataset was extracted to assess the LASSO prognostic model (Table 5). We initially performed univariate Cox regression to evaluate the prognostic value of crucial ARGs (CD24, KRT15, MIA, NDRG1, and TP63) in patients and then formulated a multivariate Cox regression model comprising of those factors with a p-value < 0.1 in the univariate Cox regression (Table 6). Forest plots revealed that two genes (CD24 and NDRG1) were associated with increased risk with HRs > 1, while three other genes (KRT15, MIA, and TP63) were protective genes with HRs < 1 (Fig. 5A). Subsequently, we constructed nomograms for ARGs to evaluate the predictive power of the Cox regression models (Fig. 5B). Furthermore, we performed calibration curves for 2-, 3-, and 4-year survival predictive power of Cox regression model nomograms, respectively (Figs. 5C–5E). The results indicated that the blue line corresponding to the 3-year survival was closest to the ideal diagonal 45° line (grey line), implying that the Cox regression models yielded the most reliable prediction for the 3-year survival. Finally, DCA was used to assess the clinical applicability of the LASSO-Cox prognostic model for the 2-, 3-, and 4-year survival (Figs. 5F–5H), and the blue line in the graphs represents the model’s predictive power. The range on the x-axis of the blue line was higher than that of the red line (all positive) and grey line (all negative) for the 3-year survival, indicating the best clinical applicability in terms of the 3-year survival.

mRNA levels of ARGs in MCF-10A, MDA-MB-231, HCC1954, and MCF-7 cell lines

Gene mRNA levels of key ARGs in normal mammary and breast cancer cell lines were detected using qPCR. We ascertained that CD24 expression was up-regulated in breast cancer cells compared with normal breast epithelial cells (Fig. 6A). Moreover, KRT15, MIA, and TP63 levels significantly decreased in breast cancer cells (Figs. 6B, 6C and 6E). However, mRNA levels of NDRG1 varied across the three breast cancer lines, being up-regulated in MDA-MB-231 cells and down-regulated in HCC1954 and MCF-7 cells (Fig. 6D).

CD24 and NDRG1 protein levels in breast cancer tissues

The expression levels of CD24 and NDRG1 that had an HR of >1 in the multivariable Cox regression model were detected using IHC based on the HPA database. It was observed, in comparison to normal breast tissues, there was a significant increase in the CD24 and NDRG1 expression in breast cancer (Figs. S3A–S3D).