Development of a prognostic model based on the ceRNA network in Triple-Negative Breast cancer

Collections of data from Gene Expression Omnibus and The Cancer Genome Atlas

The eligible GEO should conform to the following criteria: (1) studies with TNBC and normal tissues; (2) the information of platforms and studies were qualified for analysis; (3) the dataset type was microarray.

We retrieved five associated microarray datasets from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/). There were one circRNA expression dataset (GSE101123), two lncRNA expression datasets (GSE64790 and GSE115275), and two mRNA expression datasets (GSE38959 and GSE53752). One hundred forty-four samples were identified, including 94 TNBC and 50 para-cancer tissue samples.

To further analyze the impact of differentially expressed genes (DEGs) on TNBC patients, RNA sequencing expression data as well as clinical information for mRNA, miRNA, and lncRNA were downloaded from The Cancer Genome Atlas (TCGA) database for 115 TNBC cases, 989 non-TNBC cases, and 113 paracancerous tissues (https://www.cancer.gov/ccg/research/genome-sequencing/tcga).

The selection criteria for clinical information of TNBC patients in the TCGA cohort included: (1) histologically confirmed TNBC with negative ER, PR, and HER2 status; (2) complete clinical follow-up data ≥ 90 days; (3) clinical information such as age, tumor stage, treatment history, and survival outcomes (overall survival time and status).

Analysis of the differently expressed lncRNA, circ-RNAs, and mRNAs by microarrays

The differential expression analysis was conducted using multiple microarray datasets from the GEO database, where circRNA analysis utilized dataset GSE101123 with annotations from the GPL19978 platform (Agilent-069978 Arraystar Human CircRNA microarray V1). LncRNA expression was analyzed using two independent datasets: GSE64790 (annotated using GPL19612; Agilent-062978 Human lncRNA v4 Microarray) and GSE115275 (annotated using GPL21827; Agilent-079487 Human lncRNA v4 Microarray), while mRNA expression profiling integrated data from GSE38959 and GSE53752, annotated using GPL4133 (Agilent-014850 Whole Human Genome Microarray 4x44K G4112F) and GPL7264 (Agilent-012097 Human 1A Microarray G4110B) platforms respectively. To mitigate batch effects and technical variations arising from different platforms and experimental time points, we implemented batch correction using the Sva package in R (version 4.3.0; R Core Team, 2023), where the ComBat function was applied with default parameters to harmonize the data across different batches while preserving biological variation. Differential expression analysis was conducted with the limma package in R-Bioconductor, including log2 transformation of expression values, linear model fitting using the lmFit function, and empirical Bayes statistics using the eBayes function with default parameters. The —log FC— > 1 and FDR-P values < 0.05 were considered as the available threshold for identifying differential expression genes.

Construction of the competing endogenous RNA (circRNA/miRNA/mRNA/LncRNA)

We systematically predicted RNA interactions to construct a ceRNA network based on the differential expression results from GEO datasets. The DEmRNA-miRNA interaction pairs were predicted through the integrated analysis of TargetScan (version 7.2, https://www.targetscan.org/vert_80/) and miRDB (version 5.0, https://mirdb.org/), with the intersection of both databases’ predictions being used to minimize false positives. The DElncRNA-miRNA interactions were identified using miRcode (version 11; http://www.mircode.org/). For DEcircRNA-miRNA interaction pairs, StarBase (Version 2.0; https://rnasysu.com/encori/) was employed to identify interaction pairs. The ceRNA network was visualized using Cytoscape (version 3.6.1) with RNA type-specific color coding (circRNAs: red circles, miRNAs: light blue diamond, mRNAs: green triangle, lncRNAs: Blue square).

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analysis of differently expressed mRNA in the ceRNA network

We performed a comprehensive functional enrichment analysis of the differentially expressed mRNAs identified in the ceRNA network using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations. The analysis was conducted using the clusterProfiler R package with Bioconductor 3.18. For GO analysis, we utilized the org.Hs.eg.db annotation package with GO database version, examining all three GO categories: Biological Process (BP), Molecular Function (MF), and Cellular Component (CC). We employed the enrichKEGG function with parameters organism = “hsa” and default parameters for KEGG analysis. The results were visualized using the enrichplot package. P values < 0.05 were considered the threshold for identifying significant GO and KEGG.

Screening the co-DEGs both in ceRNA and TCGA

Initial data preprocessing included filtering to retain genes with mean expression values > 0.2. The R-Bioconductor limma package was applied to analyze the expression of genes from TCGA, and the | log FC| > 1 and FDR-P values < 0.05 were considered as the available threshold for identifying differential expression genes between TNBC/non-TNBC cancer tissues and para-cancer tissues. Venn intersection analysis determined the co-DEGs in TCGA and ceRNA.

Development of a prognostic model

Clinical information and expression data for the co-differentially expressed genes were collected from the TCGA. The clinical data included age, tumor stage, treatment history, and survival outcomes (overall survival time and status). Data preprocessing involved removing samples with missing survival information or a follow-up time of less than 90 days. Multivariate and univariate Cox proportional hazards regression analysis was performed using the R survival and gtsummary package to develop a prognostic index. Patients were stratified into high-risk as well as low-risk groups using the median risk score as the cutoff. Survival analysis was conducted using the survival package and survminer. The survivalROC and PRROC package was employed to examine the predictive accuracy of the prognostic model. Additionally, we developed three distinct prognostic models: Model I (9-gene signature): SH3BGRL2, CA12, LRP8, NAV3, GFRA1, DCDC2, CDC7, ABAT, and NPTX1. Model II (6-gene signature): SH3BGRL2, CA12, NAV3, GFRA1, ABAT, and NPTX1.Model III (5-gene signature): SH3BGRL2, CA12, NAV3, CDC7, and ABAT.

Gene set enrichment analysis of DEGs associated with a prognostic model

To investigate the underlying biological mechanisms of the prognostic index, we performed Gene Set Enrichment Analysis (GSEA) using GSEA software version 4.0.3. The analysis used the Molecular Signatures Database (MSigDB) version 6.0, specifically the c2.cp.kegg.v6.0.symbols.gmt dataset, which contains 186 curated KEGG pathway gene sets. Gene set size filters ≥ 15. Statistical significance was assessed using normalized enrichment score (NES), nominal P-value < 0.05, and false discovery rate (FDR) < 0.25. The results were visualized using enrichment plots showing the running enrichment score and positions of gene set members on the ranked list.

Statistical analysis

Categorical variables were presented as frequencies and percentages (%). The normality of continuous variables was assessed using the Shapiro–Wilk test. Normally distributed continuous variables were expressed as mean ± standard error (SE), while non-normally distributed continuous variables were presented as median with interquartile range (IQR, 25th–75th percentiles). Between-group comparisons for categorical variables were performed using Pearson’s chi-square test or Fisher’s exact test. An independent samples t-test was used for continuous variables with normally distributed data, while the Mann–Whitney U test was applied for data that were not non-normally distributed. All statistical tests were two-sided, and P < 0.05 was considered statistically significant. All statistical analyses and visualizations were conducted using the R (version 4.3.0).

Analysis of differentially expressed mRNA/lncRNA/circRNAs in TNBC compared to normal breast tissue from the GEO database

After downloading the raw data and platform information from the GEO database, two different mRNAs/lncRNAs datasets were normalized in batch. Then, they merged to increase the number of patients and data reliability. Based on the cutoff value with —log FC— > 1, P values < 0.0. One hundred forty differentially expressed circRNAs, 493 differentially expressed mRNAs, and 2,998 differently expressed LncRNAs were identified, and the heatmap demonstrated distinct clustering patterns between TNBC and paired non-tumorous tissues based on the expression profiles of circRNAs, mRNAs, and LncRNAs (Fig. 1).

Construction of CeRNA (CircRNAs/LncRNAs/miRNAs/mRNAs) network

To elucidate the complex regulatory interactions among differentially expressed RNAs in TNBC, we constructed a ceRNA network. The initial prediction analysis using the miRcode database identified significant interactions between 34 differentially expressed lncRNAs and 203 miRNAs (Table S1). Further analysis through the integration of TargetScan and miRDB databases revealed regulatory relationships between differentially expressed 138 mRNAs and 193 miRNAs (Table S2). Using the StarBase database, we further identified additional interactions involving 20 differentially expressed circRNAs and 167 miRNAs (Table S3). After integrating these results and filtering for common interactions, we finalized a ceRNA network comprising 15 circRNAs, 34 lncRNAs, 78 miRNAs, and 107 mRNAs. This interconnected regulatory network was visualized using Cytoscape software to illustrate the complex RNA-RNA interactions in TNBC, while the connection degree of each gene was calculated to demonstrate its contribution in the ceRNA network (Fig. 2 and Table S5).

GO and KEGG pathway analysis of DEGs in CeRNA

To investigate the biological functions and signaling pathways associated with the DEGs in the ceRNA network, we performed Gene GO term and KEGG pathway analyses. GO analysis revealed significant enrichment in several key biological processes (P < 0.05), including negative regulation of protein phosphorylation, negative regulation of phosphorylation, cellular response to alcohol, oxygen level sensing, metal ion response, hypoxia adaptation, and negative regulation of the MAPK cascade (Fig. 3A). KEGG pathway analysis identified four significantly enriched pathways (P < 0.05): the MAPK signaling pathway, microRNAs in cancer, transcriptional misregulation in cancer, and the FOXO signaling pathway (Fig. 3B). Notably, the MAPK signaling pathway was the most significantly enriched (P = 0.005), suggesting its potential central role in TNBC progression through the ceRNA network regulation.

Identification of co-expressed DEGs in CeRNA and TCGA

To identify TNBC-specific gene signatures, we conducted differential expression analysis comparing TNBC samples with non-TNBC breast cancer tissues and normal tissues using RNA-sequencing data from the TCGA database (Figs. 4A and 4B). Through filtering criteria (—log FC— > 1, FDR-adjusted P < 0.05), a total of 1,981 DEGs exclusively expressed in TNBC were identified (Fig. 4C). Subsequently, a Venn diagram analysis was performed to identify genes in our constructed ceRNA network that were differentially expressed in TNBC samples. This intersection analysis revealed nine key genes: SH3BGRL2, CA12, LRP8, NAV3, GFRA1, DCDC2, CDC7, ABAT, and NPTX1 (Fig. 4D), suggesting their potential role as TNBC-specific molecular markers.

Development of a prognostic model utilizing the nine DEGs

One hundred and four TNBC patients with an overall survival time exceeding 90 days were included for subsequent analysis. As demonstrated in Table 1, significant differences were observed in the clinical stage (P = 0.005) and chemotherapy history (P = 0.02) exhibited between the survival and mortality groups. Furthermore, univariate and multivariate Cox regression analysis suggested that risk score derived from nine DEGs (P = 0.008) emerged as a statistically significant independent predictive factor (Table 2). The risk score formula was calculated using the following formula: Risk Score = [−0.79 × SH3BGRL2] + [0.40 × CA12] + [0.16 × LRP8] + [1.98 × NAV3] + [0.37 × GFRA1] + [0.23 × DCDC2] + [−0.55 ×CDC7] + [−1.18 × ABAT] + [−0.90 × NPTX1]. Patients were stratified into two groups using the median risk score as the threshold value (Fig. 5A). As shown in Fig. 5B, Kaplan–Meier survival analysis indicated that the low-risk score group exhibited significantly improved prognosis (P < 0.01). The area under the curve (AUC) of receiver operating characteristic (ROC) was 0.90, indicating that this prognostic signature demonstrates superior predictive capability for clinical outcomes in TNBC patients compared to alternative models with AUCs of 0.86 and 0.83 (Fig. 6A). Similarly, in the precision-recall (PR) Curve analysis, model I achieved an AUC of 0.46, outperforming Model II and matching the performance of Model III (Fig. S1). Moreover, elevated risk scores corresponded to increased mortality risk (Fig. 6B).

GSEA of the DEGs based on the prognostic model

To identify the signal pathway associated with the prognostic model, we empolyed the GSEA method to analyze the expression profile of the nine genes in the model. The results revealed that the high expression of CA12, GFRA1, and NPTX1 was significant on the common MAPK signal pathway (Fig. 7), which consisted of the KEGG analysis of CeRNA (Fig. 3B). In addition, both the GEO and TCGA datasets demonstrated significantly reduced expression of CA12 (P < 0.01), GFRA1 (P < 0.01), and NPTX1 (P < 0.01) in TNBC tissue (Fig. 8).

Identification of the CircRNAs associated with MAPK signal pathway

Based on the GSEA analysis result, we constructed a ceRNA network to identify circRNAs associated with the MAPK signal pathway (Fig. 9A and Table S4). As illustrated in Fig. 9B, hsa_circ_0005455, hsa_circ_000632, hsa_circ_0001666, and hsa_circ_0000069 emerged as key regulators of CA12, GFRA1, NPTX1, which are critical components of the MAPK signal pathway. Specifically, hsa_circ_0005455 regulates NPTX1 through hsa-miR-135a/b-5p, hsa-miR-130a/b-3p, hsa-miR-454-3p, and hsa-miR-301a/b-3p. Similarly, hsa_circ_0001666 modulated GFRA1 through hsa-miR-300 and hsa-miR-381-3p, while hsa_circ_0000632 regulated CA12 through the hsa-miR-520 family (a/b/c/d/e), hsa-miR-302 family (a/b/d/e), and hsa-miR-372/373-3p. Furthermore, cirRNAs such as hsa_circ_0005455 and hsa_circ_00000069 targets hsa-miR-135b-5p, thereby co-regulating NPTX1 expression. Our analysis confirmed the significant differential expression of seven key lncRNAs: ADAMTS9-AS1, MAGI2-AS3, ADAMTS9-AS2, DIO3OS, WDFY3-AS2, MEG3, and TPRG1-AS1 (Fig. 10). Among the 22 miRNAs identified in our network, three miRNAs (miR-139-5p, miR-130a-5p, and miR-135b-5p) demonstrated significant associations with clinical outcomes in TNBC patients (Fig. 11). Notably, validation using the TCGA database provided strong evidence for the clinical relevance of these regulatory networks. These TCGA-based findings substantially support the clinical significance of our identified regulatory network components.