Prognostic analyses of genes associated with anoikis in breast cancer

Abstract

Breast cancer (BRCA) is the most diagnosed cancer worldwide and is responsible for the highest cancer-associated mortality among women. It is evident that anoikis resistance contributes to tumour cell metastasis, and this is the primary cause of treatment failure for BRCA. However, anoikis-related gene (ARG) expression profiles and their prognostic value in BRCA remain unclear. In this study, a prognostic model of ARGs based on The Cancer Genome Atlas (TCGA) database was established using a least absolute shrinkage and selection operator analysis to evaluate the prognostic value of ARGs in BRCA. The risk factor graph demonstrated that the low-risk group had longer survival than the high-risk group, implying that the prognostic model had a good performance. We identified 11 ARGs that exhibited differential expression between the two risk groups in TCGA and Gene Expression Omnibus databases. Through Gene Ontology and Kyoto Encyclopaedia of Genes and Genomes enrichment analyses, we revealed that the screened ARGs were associated with tumour progression and metastasis. In addition, a protein-protein interaction network showed potential interactions among these ARGs. Furthermore, gene set enrichment analysis suggested that the Notch and Wnt signalling pathways were overexpressed in the high-risk group, and gene set variation analysis revealed that 38 hallmark genes differed between the two groups. Moreover, Kaplan-Meier survival curves and receiver operating characteristic curves were used to identify five ARGs (CD24, KRT15, MIA, NDRG1, TP63), and quantitative polymerase chain reaction was employed to assess the differential expression of these ARGs. Univariate and multivariate Cox regression analyses were then performed for the key ARGs, with the best prediction of 3 year survival. In conclusion, ARGs might play a crucial role in tumour progression and serve as indicators of prognosis in BRCA.

Keywords: Anoikis; Breast cancer; CD24; NDRG1; Prognosis.

Figure 1. Flowchart of the study.

The TCGA-BRCA and BRCA datasets were obtained from the TCGA and GEO databases, respectively. LASSO analysis was constructed using 100 ARGs based on the TCGA-BRCA dataset, and we screened 12 ARGs with possible prognostic values. Subsequently, all samples from TCGA-BRCA and BRCA datasets were divided into high- and low-risk groups based on the prognostic model’s median risk score. Each ARG was compared between the two risk groups, and 11 significant differentially expressed ARGs with the same expression trends in both datasets were identified. Next, GO and KEGG enrichment analyses and PPI network construction were performed using the 11 ARGs. We further performed GSEA and GSVA based on the TCGA-BRCA dataset and obtained five key ARGs by plotting KM survival curves and ROC curves (AUC > 0.6) for each ARG based on the TCGA-BRCA dataset. Finally, univariate and multivariate Cox regression analyses were performed for key ARGs, and qPCR was used to detect gene expression in different cell lines. Results of ARG protein expression with an HR of >1 were downloaded from the HPA database. TCGA, The Cancer Genome Atlas; BRCA, breast invasive carcinoma; GEO, Gene Expression Omnibus; LASSO, least absolute shrinkage and selection operator; ARGs, anoikis-related genes; GO, Gene Ontology; KEGG, Kyoto Encyclopaedia of Genes and Genomes; PPI, protein-protein interaction; GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; KM, Kaplan–Meier curve; ROC, receiver operating characteristic curve; AUC, area under the curve; qPCR, quantitative polymerase chain reaction; HR, hazard ratio; HPA, human protein atlas.

Figure 2. Prognostic model establishment and DEG screening between the high- and low-risk groups in breast cancer.

(A) A LASSO prognostic model was constructed using 100 ARGs. (B–C) Variable trajectories graph and risk factor graph of the prognostic model. (D–E) Volcanic maps of the DEGs between the two risk groups based on the TCGA-BRCA and BRCA datasets, respectively. (F–G) Heat maps of ARG expression between the two risk groups in TCGA-BRCA and BRCA datasets, respectively. (H–I) Box plots comparing each ARG between the two risk groups in TCGA-BRCA and BRCA datasets. ns, p ≥ 0.05, *, p < 0.05, **, p < 0.01, ***, p < 0.001, and p < 0.05 were considered statistically significant. DEGs, differentially expressed genes; LASSO, least absolute shrinkage and selection operator; TCGA, The Cancer Genome Atlas; BRCA, breast invasive carcinoma; ARGs, anoikis-related genes.

Figure 3. GSEA and GSVA.

(A) Four biological processes of GSEA for TCGA-BRCA dataset. Enrichment analysis of the pre-Notch expression and processing (B), signalling by Notch (C), TCF dependent signalling in response to Wnt (D), signalling by Wnt (E) in TCGA-BRCA dataset. (F) Heat maps of GSVA for TCGA-BRCA dataset. (G) Box plots of significantly differently enriched hallmark genes between the two risk groups. ***, p < 0.001 and p < 0.05 were considered statistically significant. GSEA, gene set enrichment analysis; TCGA, The Cancer Genome Atlas; BRCA, breast invasive carcinoma; TCF, T cell factor; GSVA, gene set variation analysis.

Figure 4. KM survival curves and ROC curves.

KM curves of BIRC3 (A), CD24 (B), IVL (C), KRT15 (D), MIA (E), NDRG1 (F), NOS2 (G), SERPINA1 (H), and TP63 (I). ROC curves of CD24 (J), KRT15 (K), MIA (L), NDRG1 (M), and TP63 (N). KM, Kaplan–Meier; BIRC3, baculoviral IAP repeat-containing protein 3; CD24, cluster of differentiation 24; IVL, involucrin; MIA, melanoma-inhibiting activity; KRT15, keratin 15; NDRG1, N-Myc downstream regulated 1; NOS2, nitric oxide synthase 2; SERPINA1, serpin family A member 1; TP63, tumour protein p63; ROC, receiver operating characteristic.

Figure 5. Analysis of clinical correlation with the prognostic models.

(A–B) Forest plots (A) and nomograms (B) of the multivariate Cox regression model for key ARGs. (C–E) Calibration curves of the 2- (C), 3- (D), and 4-year (E) survival for nomograms of multivariate Cox regression. Black lines at the top correspond to the distributions of predicted survival probabilities of samples. (F-H) 2- (F), 3- (G), and 4-year (H) DCA of multivariate Cox regression. ARGs, anoikis-related genes; DCA, decision curve analysis.

Figure 6. MRNA levels of ARGs in MCF-10A, MDA-MB 231, HCC1954, and MCF-7 cell lines.

mRNA levels of CD24 (A), KRT15 (B), MIA (C), NDRG1 (D), and TP63 (E) determined by qPCR. The mean ± SD is shown. * p < 0.05. mRNA, messenger ribonucleic acid; ARGs, anoikis-related genes; MCF, Michigan Cancer Foundation; CD24, cluster of differentiation 24; KRT15, keratin 15; MIA, melanoma-inhibiting activity; NDRG1, N-Myc downstream regulated 1; TP63, tumour protein p63; qPCR, quantitative polymerase chain reaction; SD, standard deviation.