Recurrence-Associated Multi-RNA Signature to Predict Disease-Free Survival for Ovarian Cancer Patients

Zhang, Yu; Ye, Qingjian; He, Junxian; Chen, Peigen; Wan, Jing; Li, Jing; Yang, Yuebo; Li, Xiaomao

doi:https://doi.org/10.1155/2020/1618527

BioMed Research International

On this page

Abstract Introduction Materials and Methods Results Discussion Conclusion Data Availability Conflicts of Interest Authors’ Contributions Supplementary Materials References Copyright Related Articles

Research Article | Open Access

Volume 2020 | Article ID 1618527 | https://doi.org/10.1155/2020/1618527

Recurrence-Associated Multi-RNA Signature to Predict Disease-Free Survival for Ovarian Cancer Patients

Yu Zhang,¹Qingjian Ye,¹Junxian He,¹Peigen Chen,¹Jing Wan,¹Jing Li,¹Yuebo Yang,¹and Xiaomao Li¹

Academic Editor: Maria Barbolina

Received30 Oct 2019

Revised12 Dec 2019

Accepted23 Dec 2019

Published15 Feb 2020

Abstract

Ovarian cancer (OvCa) is an intractable gynecological malignancy due to the high recurrence rate. Several molecular biomarkers have been previously screened for early identifying patients with a high recurrence risk and poor prognosis. However, all the known studies focused on a single type of RNAs, not integrating various types. This study was to construct a new multi-RNA-based model to predict the recurrence and prognosis for OvCa patients by using the messenger RNA (mRNA, including long noncoding RNA (lncRNA)) and microRNA (miRNA) sequencing data of The Cancer Genome Atlas database. After univariate Cox regression and least absolute shrinkage and selection operator analyses, a multi-RNA-based signature (2 miRNAs: hsa-miR-508, hsa-miR-506; 1 lncRNA: TM4SF1-AS1; 11 mRNAs: MAGI3, SLAMF7, GLI2, PDK1, ARID3A, PLEKHG4B, TNFAIP8L3, C1QTNF3, NDUFAF1, CH25H, TMEM129) was generated and used to establish a risk score model. The high- and low-risk patients classified by the median risk score exhibited significantly different recurrence risks (89% versus 61%, ) and survival time (the area under the receiver operating characteristic curve (AUC) = 0.901 for 5-year disease-free survival (DFS)). This risk model was independent of other clinical features and superior to pathologic staging for DFS prediction (AUC, 0.906 versus 0.524; C-index, 0.633 versus 0.510). Furthermore, some new interaction axes were revealed to explain the possible functions of these RNAs (competing endogenous RNA: TM4SF1-AS1-miR-186-STEAP2, LINC00536-miR-508-STEAP2, LINC00475-miR-506-TMEM129; coexpression: LINC00598-PLEKHG4B). In conclusion, this multi-RNA-based risk model may be clinically useful to stratify OvCa patients with different recurrence risks and survival outcomes and included RNAs may be potential therapeutic targets.

1. Introduction

Ovarian cancer (OvCa) is one of the frequently diagnosed, lethal malignancies in the female genital system worldwide. It was estimated that there were 22,240 new cases and 14,070 deaths in the United States in 2018 [1] while 52,100 new cases and 22,500 mortalities were reported in China in 2015 [2]. Surgical resection combined with chemotherapy is the first-line treatment for OvCa, which has been demonstrated to improve the prognosis of patients. However, the overall five-year survival rate remains lower (approximately 30%) [3, 4] because more than 70% of patients [5] may experience recurrence, ultimately resulting in the treatment failure. Therefore, early identification of patients with a high recurrence risk is necessary in order to schedule personalized therapies and improve prognostic outcomes.

Recently, with the rapid developments in molecular biology, several scholars devote themselves to identify molecular biomarkers for prediction of recurrence and survival for OvCa patients. For example, Yang et al. used the microarray data of GSE9891 and GSE30161 retrieved from the Gene Expression Omnibus (GEO) database and LASSO (least absolute shrinkage and selection operator) penalized regression analysis to identify a six-long noncoding RNA (lncRNA: including RUNX1-IT1, MALAT1, H19, HOTAIRM1, LOC100190986, and AL132709.8) recurrence signature. This signature was found to predict the disease-free survival (DFS) for OvCa patients, with the area under the receiver operating characteristic (ROC) curve (AUC) of 0.813 (GSE9891, training), 0.697 (GSE9891, internal validation), and 0.711 (GSE30161, external validation) [6]. Dong et al. applied the support vector machine (SVM) method to screen a 19-microRNA (miRNA) signature that could well distinguish the recurrent from the nonrecurrent samples. Six of them (miR-193b, miR-211, miR-218, miR-505, miR-508, and miR-514) were found to be independently related to prognosis and were used to establish a risk score which was proved to significantly discriminate the recurrence-free survival (RFS) between the high-risk and low-risk groups (AUC = 0.961 for The Cancer Genome Atlas (TCGA) sequencing dataset; AUC = 0.922 for GSE25204 microarray dataset; AUC = 0.921 for GSE27290 microarray dataset) [7]. Zhou et al. also utilized the SVM algorithm to screen 39 feature genes. They distinguished the recurrence samples from nonrecurrence samples, with the prediction accuracies of 92.7%, 93.3%, 96.6%, and 90.4% for GSE17260, GSE44104, GSE51088, and TCGA datasets, respectively. The Kaplan–Meier (K-M) survival curve revealed that the survival time of patients with predicted nonrecurrent OvCa was significantly longer than that of patients with predicted recurrent OvCa [8]. Zhao et al. screened KCNN4 and S100A14 combination as the perfect recurrence prediction model (AUC = 0.5442–0.9524) [9]. However, all these studies focused on the prognostic values of single RNA type, not a combination of them [10]. The study on colorectal cancer indicated that a multi-RNA-based classifier model (AUC = 0.83) seemed to be more effective in stratifying patients who are at a high risk of mortality than lncRNAs (AUC = 0.56, ), miRNAs (AUC = 0.67, ), or mRNAs alone (AUC = 0.76, ) [11]. Therefore, it is of clinical significance to develop a multi-RNA-based prognostic signature for OvCa.

In this study, we aimed to collect the lncRNA, miRNA, and mRNA expression data from the TCGA database and construct a multi-RNA-based risk score model to predict the recurrence and DFS for OvCa patients. In addition, the underlying functions of lncRNAs, miRNAs, and mRNAs were predicted by constructing the coexpression, competing endogenous RNA (ceRNA) [12], and protein-protein interaction (PPI) [8] networks, which was also not simultaneously analyzed in previous studies on the recurrence signature of OvCa. Accordingly, we hypothesize that our study may provide a more effective, function-clear signature for the prediction of recurrent prognosis in OvCa patients.

2. Materials and Methods

2.1. Patient Datasets

The mRNA (including lncRNA, n = 379) and miRNA (n = 498) expression data (level 3) with their clinical information were collected by searching the TCGA database portal (https://portal.gdc.cancer.gov/) on August 6, 2019, with the keyword of ovarian cancer. These expression profiles were obtained by sequencing on the Illumina HiSeq 2000 platform. A total of 322 samples were retained for the following analysis because they satisfied the following inclusion criteria: (1) mRNAs and miRNAs were both sequenced; (2) information of recurrence (n = 80, nonrecurrent; n = 242, recurrent) and survival status was clearly described. Due to the fact that no other datasets simultaneously met the above two inclusion criteria, these 322 patients, obtained from the TCGA database, were randomly divided to training and validation sets to achieve internal validation. Their detailed clinical characteristics are presented in Table 1.

2.2. Differential Expression Analyses

The mRNAs, lncRNAs, and miRNAs in RNA-seq profiles were annotated according to the annotation information of the HUGO Gene Nomenclature Committee (HGNC; http://www.genenames.org/) that includes 4,422 lncRNAs, 19,223 protein-coding genes, and 1,914 miRNAs [13]. The differentially expressed genes (DEGs), lncRNAs (DELs), and miRNAs (DEMs) between recurrent and nonrecurrent samples were identified using the Linear Models for Microarray Data (LIMMA) method (version 3.34.7; https://bioconductor.org/packages/release/bioc/html/limma.html) [14] in R package (version 3.4.1; http://www.R-project.org/). RNAs with false discovery rate (FDR) < 0.05 and |logFC(fold change)| > 0.5 were considered to be significant. Heat map of differentially expressed RNAs (DERs) was plotted with “pheatmap” package (version: 1.0.8; https://cran.r-project.org/web/packages/pheatmap) based on centered Pearson’s correlation.

2.3. Construction of Multi-RNA-Based Prognostic Signature

Univariate Cox regression analysis was first performed to preliminarily screen the DEGs, DELs, and DEMs that were associated with DFS using the “survival” package in R (version, 2.41-1; http://bioconductor.org/packages/survivalr/). Subsequently, a multivariate Cox regression analysis was conducted to confirm their independence. Log-rank p value < 0.05 was set as the statistical threshold for these two analyses. In order to further optimize the prognostic signature, an L1-penalized estimation-based Cox proportional hazards regression model (that is, LASSO) in the penalized package (version, 0.9-5; http://bioconductor.org/packages/penalized/) [15, 16] was applied for the independent prognostic DERs, in which the value of penalty parameter lambda was chosen via cross-validation likelihood routine (1000 iterations). Then, the prognostic risk score was established according to the expression levels of the RNAs () and prognostic coefficients ():

2.4. Assessment of the Prognostic Performance of the Risk Score Model

Using the median risk score as the cutoff point, the OvCa patients were classified into high-risk and low-risk groups. K-M survival curve was drawn by using the “survival” package in R to observe the DFS differences between high- and low-risk groups. Furthermore, the ROC curves with the calculation of AUC were also plotted using the pROC in R (version, 1.14.0; https://cran.r-project.org/web/packages/pROC/index.html) to estimate the prediction accuracy of this risk score for the 1-, 3-, and 5-year survival rate of high- and low-risk patients. These analyses were performed for each dataset, respectively.

To verify that the multi-RNA signature was an independent prognostic factor, univariate and multivariate Cox regression analyses were performed for DFS using the risk score and several clinical features (including age, neoplasm histologic grade, pathologic stage, lymphovascular invasion, and vascular invasion) in each dataset. The superiority of the risk score model to clinical factors in survival prediction was also validated by K-M curve analysis for subgroups, calculation of AUC, and Harrell’s concordance index (C-index) for each classifier. C-index was calculated using survcomp package in R (http://www.bioconductor.org/packages/release/bioc/html/survcomp.html) [17].

2.5. Function Prediction of the Prognostic RNAs

Coexpression, ceRNA, and PPI networks were constructed to predict the possible functions of signature lncRNAs, miRNAs, and mRNAs, respectively. Pearson correlation coefficients (PCC) between DELs and DEGs were calculated to represent their possible coexpression relationships using cor.test function (https://stat.ethz.ch/R-manual/R-devel/library/stats/html/cor.test.html) in R. Potential interactions between DEMs and DELs as well as between DEMs and DEGs were predicted using DIANA-LncBase (version, 2.0; http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=lncbasev2/index-predicted) [18] and starBase database (version, 2.0; http://starbase.sysu.edu.cn/) [19], respectively. Only the DEL-DEM and DEM-DEG interactors that had an opposite expression trend were used to construct the ceRNA network [20]. The interactions between DEGs were predicted by mapping the genes to the Search Tool for the Retrieval of Interacting Genes (STRING; version, 10.0; http://string-db.org/) database [21]. PPI network was constructed based on interaction pairs with a combined score >0.4. All networks were visualized using Cytoscape (version 3.6.1; http://www.cytoscape.org/).

Function enrichment analysis was performed for the genes in all three networks using the Enrichr online tool (http://amp.pharm.mssm.edu/Enrichr/) [22], in which Kyoto Encyclopedia of Genes and Genomes (KEGG), BioCarta, Reactome Pathways and Gene Ontology (GO) biological process (BP), and molecular function (MF) terms could be obtained. A p value < 0.05 was considered statistically significant.

3. Results

3.1. Differential Expression Analysis

A total of 12,728 mRNAs, 1,214 lncRNAs, and 547 miRNAs were annotated based on HGNC. Under the threshold of |log2FC| > 0.5 and FDR < 0.05, 540 of them were screened as DERs between recurrent and nonrecurrent samples (Figure 1(a)), including 451 DEGs, 68 DELs, and 21 DEMs (Figure 1(b)). The heat map showed that these DERs can obviously distinguish the recurrent from the nonrecurrent samples (Figure 1(c)).

(a)

(b)

(c)

Figure 1

Differentially expressed RNAs. (a) Scatter diagram to show the distribution of differentially expressed RNAs (blue dot). The horizontal dotted line is FDR < 0.05; the vertical dotted line is |logFC| >0.5; (b) upregulated and downregulated number and ratio of each RNA type; (c) heat map of all differentially expressed RNAs. The horizontal axis is the sample, and the vertical axis is the expression level. Red, high expression; green, low expression. FC, fold change; FDR, false discovery rate.

3.2. Construction of a Prognostic Signature

Univariate Cox regression analysis identified 160 DERs to be significantly associated with DFS, including 140 DEGs, 10 DELs, and 10 DEMs. Multivariate Cox regression analysis filtered out 61 DERs as independent prognostic biomarkers, including 53 DEGs, 4 DELs, and 4 DEMs. Fourteen of them (2 DEMs: hsa-miR-508, hsa-miR-506; 1 DEL: TM4SF1-AS1; 11 DEGs: MAGI3, SLAMF7, GLI2, PDK1, ARID3A, PLEKHG4B, TNFAIP8L3, C1QTNF3, NDUFAF1, CH25H, TMEM129) were further selected as the optimal prognostic combination after LASSO analysis (Table 2; Figure 2(a)). Also, the results of the expression, univariate, and LASSO analyses were consistent, showing that hsa-miR-506, MAGI3, SLAMF7, GLI2, PDK1, ARID3A, and PLEKHG4B may be tumor suppressor genes and high expression of them predicted excellent prognosis while hsa-miR-508, TM4SF1-AS1, TNFAIP8L3, C1QTNF3, NDUFAF1, CH25H, and TMEM129 may be oncogenic genes and poor prognosis would occur in patients with high levels of them (Table 2). The prognostic risk score model was then established based on the expression of these biomarkers and their LASSO coefficients (Table 2; Figure 2(a)): prognostic risk score = (−0.20 hsa-miR-506) + (−0.27 MAGI3) + (−0.29 SLAMF7) + (−0.12 GLI2) + (−0.08 PDK1) + (−0.02 ARID3A) + (−0.03 PLEKHG4B) + (0.12 hsa-miR-508) + (0.04 TM4SF1-AS1) + (0.29 TNFAIP8L3) + (0.13 C1QTNF3) + (0.15 NDUFAF1) + (0.16 CH25H) + (0.14 TMEM129).

3.3. Assessment of the Prognostic Performance of This Signature

The risk score was calculated for each patient in each dataset (Supplemental Information 1), and the patients in each dataset were divided into a high-risk group and low-risk group according to their median risk score. As shown in Figure 2(b), the patients with high-risk scores were shown to be at a high risk of recurrence (72/81, 89% versus 49/80, 61%; Chi-square = 16.5, ) and had shorter DFS than those with the low-risk scores (Figure 2(c)). Also, patients with high-risk scores tended to express oncogenic RNAs, whereas tumor suppressor RNAs inclined to be highly expressed in patients with low-risk scores (Figure 2(d)).

To further confirm the prognostic performance of this signature, the K-M survival curve was drawn for the training (Figure 3(a)), validation (Figure 3(c)), and entire sets (Figure 3(e)). As expected, the DFS of the patients in the high-risk group was significantly shorter than that of the low-risk group in all sets. ROC curve analysis also validated that the predictive accuracy of this signature seemed to be relatively high, with the AUC of 0.901, 0.933, and 0.987 for the 5-year survival rate in the training (Figure 3(b)), validation (Figure 3(d)), and entire sets (Figure 3(f)), respectively.

(a)

(b)

(c)

(d)

(e)

(f)

Univariate and multivariate Cox regression analyses demonstrated that the prognostic value of this multi-RNA-based signature was independent of other clinical features in the training and entire datasets (Table 3). Furthermore, the pathologic stage was also found to be an independent prognostic factor for OvCa patients, with a higher stage having poor DFS (Table 3; Figure 4). To investigate whether our multi-RNA-based signature added additional prognostic values to the commonly used pathologic stage, stratification analyses were also conducted. The results showed this multi-RNA-based signature can further distinguish the survival of patients at the same stage 3 (Figure 4). The survival of patients with stages 2 and 4 not significantly separated by this multi-RNA-based signature may be related to small sample size in these subgroups.

(a)

(b)

(c)

(d)

The superiority of the risk score model to clinical factors in survival prediction was also validated by comparison of AUC and C-index for each classifier. As anticipated, the results revealed that the AUC of the multi-RNA-based model was relatively higher than that of the pathologic stage and single RNA, but nearly similar to the combined model using the multi-RNA and pathologic stage (Figure 5). These findings implied that our multi-RNA based signature may be a promising biomarker for clinical use to predict the DFS of OvCa patients.

3.4. Functional Analysis for Prognostic Signature

To explore the underlying molecular mechanisms of these 14 RNAs in the signature, coexpression network between DELs and DEGs, ceRNA network among DELs, DEMs, and DEGs, and PPI networks between DEGs were constructed.

A total of 1086 coexpression pairs were identified between 67 DELs and 191 DEGs and were used to establish the coexpression network (Figure 6), in which 14 prognostic genes were also found to have positive coexpression with several lncRNAs, such as LINC00598-GLI2/ARID3A/PLEKHG4B/MAGI3/PDK1/SLAMF7, LINC00659-TMEM129, LINC00475-TNFAIP8L3/C1QTNF3/NDUFAF1, and LINC00189-CH25H (all PCC > 0.4). Function analysis showed that ARID3A was associated with Generic Transcription Pathway_Homo sapiens_R-HSA-212436; GLI2 was enriched in the Hedgehog signaling pathway, Basal cell carcinoma, and positive regulation of DNA replication (GO:0045740); and C1QTNF3 was enriched in regulation of cytokine secretion (GO:0050707) (Table 4).

A total of 38 DELs were predicted to interact with 18 DEMs by the DIANA-LncBase database, which formed 90 interaction pairs, while 13 DEMs were predicted to interact with 66 DEGs by the StarBase database to constitute 104 interaction pairs. These DEL-DEM and DEM-DEG interaction pairs were integrated to construct the ceRNA network (Figure 7), in which prognostic signature lncRNA TM4SF1-AS1 may function as a ceRNA for miR-186 to regulate STEAP2; prognostic signature miR-508 and miR-506 related ceRNA axes were LINC00536-miR-508-STEAP2 and LINC00475-miR-506-TMEM129, respectively. Function enrichment analysis showed that STEAP2 was involved in mineral absorption, transmembrane transport of small molecules_Homo sapiens_R-HSA-382551, copper ion import (GO:0015677), copper ion transport (GO:0006825), and iron ion import across plasma membrane (GO:0098711); TMEM129 was involved in ubiquitin-protein ligase activity involved in ERAD pathway (GO:1904264); PFKFB4 was enriched in positive regulation of glycolytic process (GO:0045821) and positive regulation of coenzyme metabolic process (GO:0051197) (Table 5).

A total of 265 PPI pairs were predicted between 223 DEGs using the STRING database, which was used to construct the PPI network (Figure 8). In this network, prognostic GLI2, ARID3A, MAGI3, PDK1, and TMEM129 were included. Function enrichment analysis showed that GLI2 was involved in Hedgehog signaling pathway, Basal cell carcinoma, Hippo signaling pathway, pathways in cancer, and regulation of transcription, DNA-templated (GO:0006355); ARID3A was enriched in transcription regulatory region sequence-specific DNA binding (GO:0000976) and RNA polymerase II regulatory region sequence-specific DNA binding (GO:0000977); MAGI3 could interact with TNK2 to participate in positive regulation of protein phosphorylation (GO:0001934) and growth factor receptor binding (GO:0070851) (Table 6). Additionally, PDK1 could interact with PFKFB4 and thus may take part in the GO terms as described in the function enrichment results of the ceRNA network.

4. Discussion

In the present study, we, for the first time, constructed a multi-RNA-based signature (consisting of 11 mRNAs, 2 miRNAs, and 1 lncRNA) to predict the recurrence and DFS for OvCa patients. The high- and low-risk patients classified by the median risk score were shown to exhibit significantly different recurrence risks and DFS time. The AUC was 0.975, 0.912, and 0.901 for prediction of 1-, 3- and 5-year DFS in the training set, respectively, which seemed to be higher compared with other single RNA signatures reported previously, including Yang et al. (6-lncRNA, AUC = 0.813 at 3 year) [6], Bagnoli et al. (35-miRNA, AUC = 0.72) [23], and Zhao et al. screened (2-mRNA, AUC = 0.6075) [9]. This superiority of multi-RNA risk score to single RNA was also validated in our study according to the AUC (0.906 versus 0.621 for lncRNA, 0.666 for miRNAs, and 0.843 for mRNAs) and C-index results (0.633 versus 0.502 for lncRNA, 0.545 for miRNAs, and 0.618 for mRNAs) and was in line with the study of Xiong et al. [11].

Clinically, the pathologic staging is commonly considered as the gold standard for the assessment of the recurrence risk [24] and prognosis [25] of OvCa patients. This conclusion was also confirmed in our study, showing that pathologic staging was an independent factor for the prediction of DFS (that is, the prognosis was the worst in stage 4) after multivariate Cox regression analysis. However, some studies revealed that patients at the same pathological stage also had different recurrence risks and prognostic outcomes [26, 27]. Thus, there is a strong need to explore more effective prognostic biomarkers to independently or jointly identify the recurrence risk and prognosis of patients. Recently, several studies on other cancers had demonstrated that the risk score established by molecular signature had a similar or even higher accuracy for prognostic prediction than clinical features and the AUC (and/or C-index) was the largest for the combination of them [11, 28, 29]. Hereby, we also calculated the AUC (and/or C-index) of our multi-RNA-based risk score and independent pathologic stage factor. In accordance with previous studies [11, 28, 29], we also found that the AUC (and/or C-index) of our multi-RNA-based risk score was relatively higher than that of the pathologic stage, and the predicted and actually observed survival was significantly correlated using the risk score (), but not pathologic stage (). The stratification analysis also showed that the risk score can further distinguish the survival of patients at the same stage 3. These findings implied the insufficiency to use pathologic staging for prognosis prediction in the clinic and the advantage of our risk score. Of course, the optimum was still the combined model, with the AUC and C-index of 0.913 and 0.634, respectively.

There was no study to investigate the roles of lncRNA TM4SF1-AS1 and it may be a novel signature gene for cancer. However, we can indirectly speculate its functions by its interacted miR-186-STEAP2 axis. Similar to mRNAs, lncRNAs also harbor miRNA-response elements (MREs) and thus they can completely bind to miRNAs, leading to fewer miRNAs interacted with mRNAs, while it is well known that miRNAs can negatively regulate gene expression by interfering their translation or stability (ceRNA hypothesis). Accordingly, the upregulation of TM4SF1-AS1 identified in our recurrent samples may result in the overexpression of STEAP2 and possible lower expression of miR-186 in OvCa. This theory was validated in our study and previous studies: extensive evidence had suggested that miR-186 was a tumor suppressor gene, with a downregulated expression level in cancer tissues, cells [30, 31], and blood of recurrent oral squamous cell carcinoma patients compared with controls [32], including OvCa [33]. Overexpression of miR-186 into cancer cells significantly inhibited proliferation, invasion, metastasis [29, 34, 35], and induced mesenchymal-to-epithelial transition, G1 cell cycle arrest, and cell apoptosis [33]. STEAP2 (Seven Transmembrane Epithelial Antigen of the Prostate 2) is a gene primarily expressed in the prostate and the ovary is the only other area with a significant expression [36]. The roles of STEAP2 in OvCa may be similar to prostate cancer in which STEAP2 was observed to be upregulated and associated with advanced stage and Gleason score [37, 38]. The overexpression of STEAP2 induced the normal prostate cells to gain an ability to migrate and invade [36] while the knockdown of STEAP2 significantly reduced proliferation and migration of prostate cancer cells [39]. However, no study was performed to verify the interaction between miR-186-STEAP2 and TM4SF1-AS1, which may be the direction of our future search.

In addition to miR-186, we found that prognostic miR-508 also targeted STEAP2 and, hereby, miR-508 may be possibly downregulated in OvCa, which was also proved in our study and another literature [40]. Transfection of cancer cells with miR-508 mimics significantly suppressed cell proliferation, migration, and invasion while the use of miR-508 inhibitors resulted in an opposite trend [40, 41]. More interestingly, LINC00536 may be an upstream gene to influence the interaction between miR-508 and STEAP2 by acting as a ceRNA in the present study. Although their interaction was not reported, the opposed expression trend and roles of LINC00536 indirectly supported its potential relationship with miR-508. LINC00536 was identified to be highly expressed in bladder cancer tissues compared with controls and negatively associated with patients’ survival rate. Function assays indicated that the knockdown of LINC00536 attenuated the malignant cell phenotypes in vitro and inhibited bladder cancer growth in vivo [42].

The present study showed that prognostic miR-506 was involved in OvCa by impacting the LINC00475-miR-506-TMEM129 ceRNA axis. This was also a novel mechanism explained for OvCa because no study reported their interactions previously, except that the expression level of each gene was preliminarily investigated as follows: Nam et al. revealed that the expression of miR-506-3p was significantly decreased in recurrent samples of OvCa patients compared with normal ovarian tissue and primary tumors [43]. The high level of miR-506 was positively associated with the early FIGO stage and longer survival in OvCa [44]. The introduction of miR-506 in OvCa cells can inhibit its proliferation and dissemination and promote senescence [44, 45]. Mo et al. used the lncRNA microarray to identify the high expression of LINC00475 in gastric cancer tissues compared with paired nontumor tissues [46]. The study of Hou et al. revealed that LINC00475 was associated with the overall survival of patients with clear cell renal cell carcinoma [47]. Silencing of LINC00475 suppressed cell proliferation, migration, and invasion in renal cell carcinoma [48] and glioma cells [49] in vitro. TMEM129 encodes an E3 ubiquitin ligase that was reported to mediate endoplasmic reticulum-associated HLA class I degradation [50] while tumors with downregulation of MHC class I conferred a significantly poorer prognosis in OvCa patients [51]. These findings suggested that TMEM129 may be upregulated in OvCa. In line with these studies, we also found that miR-506 was downregulated and LINC00475 and TMEM129 were upregulated in recurrent samples compared with nonrecurrent controls.

As for the prognostic genes not explored in the ceRNA network, most of them had been shown to be associated with the development and progression of OvCa or other cancers previously in the literature. For example, C1QTNF3 (C1q and TNF related 3) was found to be upregulated in bowel metastasis samples than primary OvCa [52]. ARID3A (AT-rich interactive domain 3A) positivity was implied to be correlated with longer DFS and cancer-specific survival in patients with residual rectal cancer [53]. MAGI3 (PDZ domain-containing protein membrane-associated guanylate kinase inverted 3) was identified to be downregulated in glioma samples. Overexpression of MAGI3 inhibited proliferation, migration, and cell cycle progression of glioma cells and decreased subcutaneous tumor growth in mice by inactivation of Wnt/β-catenin signaling pathway. High expression of MAGI3 was also significantly associated with excellent overall survival [54]. TNFAIP8L3 (TNF alpha-induced protein 8 like 3, also known as TIPE3), a transfer protein of phosphoinositide second messengers, was detected to be significantly upregulated in breast cancer tissues (especially invasive or metastasized type) as compared with adjacent nontumor tissues. Inhibition of TNFAIP8L3 may block cancer cell proliferation, migration, and invasion in vitro and in vivo by inactivation of the AKT and NF-κB signaling pathways [55, 56]. In consistence with these studies, we also found that C1QTNF3 and TNFAIP8L3 were upregulated, while ARID3A and MAGI3 were downregulated in recurrent samples, and they were respectively an oncogenic or tumor suppressor gene for prediction of DFS. Although there were also studies to confirm the roles of Gli2 (Glioma-associated oncogene family member 2) [57], PDK1 (pyruvate dehydrogenase kinase) [58], NDUFAF1 (NADH dehydrogenase 1 alpha subcomplex assembly factor 1) [59], SLAMF7 (SLAM family member 7) [60], and CH25H (cholesterol 25-hydroxylase) [61], their conclusions seemed to be opposite with our study, which may be attributed to the differences in cancer type or sample size and thus further validation experiments are needed. PLEKHG4B was not investigated in any study of cancer, but we speculated that its expression may be positively regulated by LINC00598 because of their PCC of 0.43. Knockdown of LINC00598 was previously proved to induce G0/G1 cell cycle arrest and inhibit proliferation, indicating its tumor-promoting roles [62]. Accordingly, we speculated that PLEKHG4B was also highly expressed in OvCa, which was verified in our study as expected.

5. Conclusion

In the present study, we used the recurrence-associated genes to construct a multi-RNA-based risk score model which was shown to effectively stratify OvCa patients into groups at low risk and high risk of shorter DFS. This model was not only independent of clinical features but also superior to commonly used pathologic staging in the clinic for prognostic prediction. Furthermore, we predicted several new interaction axes to explain the possible mechanisms of these RNAs in OvCa, such as TM4SF1-AS1-miR-186-STEAP2, LINC00536-miR-508-STEAP2, LINC00475-miR-506-TMEM129, and LINC00598-PLEKHG4B. However, some limitations still exist. First, quantitative PCR experiments should be performed to validate the expression and prognosis accuracy of our signature genes in clinical samples because there may be potential differences with TCGA sequencing data results as reported for Gli2, PDK1, NDUFAF1, SLAMF7, and CH25H. Second, in vitro and in vivo experiments are also essential to validate the ceRNA or coexpression mechanisms we identified.

Data Availability

All data can be available in TCGA database (https://portal.gdc.cancer.gov/).

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Yu Zhang and Qingjian Ye contributed equally to this work. YZ, QJY, YBY, and XML conceived the study and participated in its design. YZ and QJY collected the data and performed bioinformatic analyses. JXH, PGC, JW, and JL searched the references and were involved in the interpretation of the results. YZ and QJY drafted the manuscript. YBY and XML revised the manuscript. All authors have read and approved the final version of the article.

Supplementary Materials

Supplemental Information 1: risk score calculated by each patient in three datasets. (Supplementary Materials)

References

R. L. Siegel, K. D. Miller, and A. Jemal, “Cancer statistics, 2018,” CA: A Cancer Journal for Clinicians, vol. 68, no. 1, pp. 7–30, 2018.
View at: Publisher Site | Google Scholar
W. Chen, R. Zheng, P. D. Baade et al., “Cancer statistics in China, 2015,” CA: A Cancer Journal for Clinicians, vol. 66, no. 2, pp. 115–132, 2016.
View at: Publisher Site | Google Scholar
M. K. Edwards, M. C. Noer, C. D. Sperling et al., “Survival of ovarian cancer patients in Denmark: results from the Danish gynaecological cancer group (DGCG) database, 1995–2012,” Acta Oncologica, vol. 55, no. 2, pp. 36–43, 2007.
View at: Publisher Site | Google Scholar
M. Timmermans, G. S. Sonke, K. K. Van de Vijver, M. A. van der Aa, and R. F. P. M. Kruitwagen, “No improvement in long-term survival for epithelial ovarian cancer patients: a population-based study between 1989 and 2014 in The Netherlands,” European Journal of Cancer, vol. 88, pp. 31–37, 2018.
View at: Publisher Site | Google Scholar
I. Vistad, L. Bjørge, O. Solheim et al., “A national, prospective observational study of first recurrence after primary treatment for gynecological cancer in Norway,” Acta Obstetricia et Gynecologica Scandinavica, vol. 96, no. 10, pp. 1162–1169, 2017.
View at: Publisher Site | Google Scholar
K. Yang, Y. Hou, A. Li et al., “Identification of a six-lncRNA signature associated with recurrence of ovarian cancer,” Scientific Reports, vol. 7, no. 1, p. 752, 2017.
View at: Publisher Site | Google Scholar
J. Dong and M. Xu, “A 19-miRNA Support Vector Machine classifier and a 6-miRNA risk score system designed for ovarian cancer patients,” Oncology Reports, vol. 41, pp. 3233–3243, 2019.
View at: Publisher Site | Google Scholar
J. Zhou, L. Li, L. Wang, X. Li, H. Xing, and L. Cheng, “Establishment of a SVM classifier to predict recurrence of ovarian cancer,” Molecular Medicine Reports, vol. 18, pp. 3589–3598, 2018.
View at: Publisher Site | Google Scholar
H. Zhao, E. Guo, T. Hu et al., “KCNN4 and S100A14 act as predictors of recurrence in optimally debulked patients with serous ovarian cancer,” Oncotarget, vol. 7, no. 28, pp. 43924–43938, 2016.
View at: Publisher Site | Google Scholar
S. Shukla, J. R. Evans, R. Malik et al., “Development of a RNA-seq based prognostic signature in lung adenocarcinoma,” Journal of the National Cancer Institute, vol. 109, no. 1, Article ID djw200, 2017.
View at: Publisher Site | Google Scholar
Y. Xiong, R. Wang, L. Peng et al., “An integrated lncRNA, microRNA and mRNA signature to improve prognosis prediction of colorectal cancer,” Oncotarget, vol. 8, no. 49, pp. 85463–85478, 2017.
View at: Publisher Site | Google Scholar
X. Wang, L. Han, L. Zhou, L. Wang, and L. M. Zhang, “Prediction of candidate RNA signatures for recurrent ovarian cancer prognosis by the construction of an integrated competing endogenous RNA network,” Oncology Reports, vol. 40, no. 5, pp. 2659–2673, 2018.
View at: Publisher Site | Google Scholar
S. Povey, R. Lovering, E. Bruford, M. Wright, M. Lush, and H. Wain, “The HUGO gene Nomenclature committee (HGNC),” Human Genetics, vol. 109, no. 6, pp. 678–680, 2001.
View at: Publisher Site | Google Scholar
M. E. Ritchie, B. Phipson, D. Wu et al., “Limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic Acids Research, vol. 43, no. 7, p. e47, 2015.
View at: Publisher Site | Google Scholar
J. J. Goeman, “L1 penalized estimation in the Cox proportional hazards model,” Biometrical Journal, vol. 52, no. 1, pp. 70–84, 2010.
View at: Publisher Site | Google Scholar
R. Tibshirani, “The lasso method for variable selection in the Cox model,” Statistics in Medicine, vol. 16, no. 4, pp. 385–395, 1997.
View at: Publisher Site | Google Scholar
M. S. Schröder, A. C. Culhane, J. Quackenbush, and B. Haibe-Kains, “survcomp: an R/Bioconductor package for performance assessment and comparison of survival models,” Bioinformatics, vol. 27, no. 22, pp. 3206–3208, 2011.
View at: Publisher Site | Google Scholar
M. D. Paraskevopoulou, G. Georgakilas, N. Kostoulas et al., “DIANA-LncBase: experimentally verified and computationally predicted microRNA targets on long non-coding RNAs,” Nucleic Acids Research, vol. 41, no. D1, pp. D239–D245, 2013.
View at: Publisher Site | Google Scholar
J.-H. Yang, J.-H. Li, P. Shao, H. Zhou, Y.-Q. Chen, and L.-H. Qu, “starBase: a database for exploring microRNA-mRNA interaction maps from Argonaute CLIP-Seq and Degradome-Seq data,” Nucleic Acids Research, vol. 39, no. 1, pp. D202–D209, 2011.
View at: Publisher Site | Google Scholar
J. Xu, J. Zhang, F. Shan, J. Wen, and Y. Wang, “SSTR5-AS1 functions as a ceRNA to regulate CA2 by sponging miR-15b-5p for the development and prognosis of HBV-related hepatocellular carcinoma,” Molecular Medicine Reports, vol. 20, no. 6, pp. 5021–5031, 2019.
View at: Publisher Site | Google Scholar
D. Szklarczyk, A. Franceschini, S. Wyder et al., “STRING v10: protein-protein interaction networks, integrated over the tree of life,” Nucleic Acids Research, vol. 43, no. D1, pp. D447–D452, 2015.
View at: Publisher Site | Google Scholar
M. V. Kuleshov, M. R. Jones, A. D. Rouillard et al., “Enrichr: a comprehensive gene set enrichment analysis web server 2016 update,” Nucleic Acids Research, vol. 44, no. W1, pp. W90–W97, 2016.
View at: Publisher Site | Google Scholar
M. Bagnoli, S. Canevari, D. Califano et al., “Development and validation of a microRNA-based signature (MiROvaR) to predict early relapse or progression of epithelial ovarian cancer: a cohort study,” The Lancet Oncology, vol. 17, no. 8, pp. 1137–1146, 2016.
View at: Publisher Site | Google Scholar
C. Chang, J. Chen, W.-A. Chen, S.-P. Ho, W. S. Liou, and A. J. Chiang, “Assessing the risk of clinical and pathologic factors for relapse of borderline ovarian tumours,” Journal of Obstetrics and Gynaecology, vol. 37, no. 2, pp. 233–237, 2016.
View at: Publisher Site | Google Scholar
E. Kempf, G. Desamericq, B. Vieites et al., “Clinical and pathologic features of patients with non-epithelial ovarian cancer: retrospective analysis of a single institution 15-year experience,” Clinical and Translational Oncology, vol. 19, no. 2, pp. 173–179, 2017.
View at: Publisher Site | Google Scholar
S.-J. Baek, J.-Y. Park, D.-Y. Kim et al., “Stage IIIC epithelial ovarian cancer classified solely by lymph node metastasis has a more favorable prognosis than other types of stage IIIC epithelial ovarian cancer,” Journal of Gynecologic Oncology, vol. 19, no. 4, p. 223, 2008.
View at: Publisher Site | Google Scholar
S. Gouy, M. Saidani, A. Maulard et al., “Characteristics and prognosis of stage I ovarian mucinous tumors according to expansile or infiltrative type,” International Journal of Gynecologic Cancer, vol. 28, no. 3, pp. 493–499, 2018.
View at: Publisher Site | Google Scholar
J. Long, L. Zhang, X. Wan et al., “A four-gene-based prognostic model predicts overall survival in patients with hepatocellular carcinoma,” Journal of Cellular and Molecular Medicine, vol. 22, no. 12, pp. 5928–5938, 2018.
View at: Publisher Site | Google Scholar
X. Tian, X. Zhu, T. Yan et al., “Recurrence-associated gene signature optimizes recurrence-free survival prediction of colorectal cancer,” Molecular Oncology, vol. 11, no. 11, pp. 1544–1560, 2017.
View at: Publisher Site | Google Scholar
T. Huang, G. Wang, L. Yang et al., “MiR-186 inhibits proliferation, migration, and invasion of non-small cell lung cancer cells by downregulating Yin Yang 1,” Cancer Biomarkers, vol. 21, no. 1, pp. 221–228, 2017.
View at: Publisher Site | Google Scholar
K. Yao, L. He, Y. Gan, Q. Zeng, Y. Dai, and J. Tan, “MiR-186 suppresses the growth and metastasis of bladder cancer by targeting NSBP1,” Diagnostic Pathology, vol. 10, no. 1, p. 146, 2015.
View at: Publisher Site | Google Scholar
J. Ries, C. Baran, F. Wehrhan et al., “Prognostic significance of altered miRNA expression in whole blood of OSCC patients,” Oncology Reports, vol. 37, no. 6, pp. 3467–3474, 2017.
View at: Publisher Site | Google Scholar
X. Zhu, H. Shen, X. Yin et al., “miR-186 regulation of Twist1 and ovarian cancer sensitivity to cisplatin,” Oncogene, vol. 35, no. 3, pp. 323–332, 2016.
View at: Publisher Site | Google Scholar
X. Hua, Y. Xiao, W. Pan et al., “miR-186 inhibits cell proliferation of prostate cancer by targeting GOLPH3,” American Journal of Cancer Research, vol. 6, no. 6, pp. 1650–1660, 2016.
View at: Google Scholar
J. Li, L. Xia, Z. Zhou et al., “MiR-186-5p upregulation inhibits proliferation, metastasis and epithelial-to-mesenchymal transition of colorectal cancer cell by targeting ZEB1,” Archives of Biochemistry and Biophysics, vol. 640, pp. 53–60, 2018.
View at: Publisher Site | Google Scholar
K. P. Porkka, M. A. Helenius, and T. Visakorpi, “Cloning and characterization of a novel six-transmembrane protein STEAP2, expressed in normal and malignant prostate,” Laboratory Investigation, vol. 82, no. 11, pp. 1573–1582, 2002.
View at: Publisher Site | Google Scholar
S. E. A. Burnell, S. Spencer-Harty, S. Howarth et al., “Utilisation of the STEAP protein family in a diagnostic setting may provide a more comprehensive prognosis of prostate cancer,” PLoS One, vol. 14, Article ID e0220456, 2019.
View at: Publisher Site | Google Scholar
H. Whiteland, S. Spencer-Harty, C. Morgan et al., “A role for STEAP2 in prostate cancer progression,” Clinical & Experimental Metastasis, vol. 31, no. 8, pp. 909–920, 2014.
View at: Publisher Site | Google Scholar
S. E. A. Burnell, S. Spencerharty, S. Howarth et al., “STEAP2 knockdown reduces the invasive potential of prostate cancer cells,” Scientific Reports, vol. 8, no. 1, p. 6252, 2018.
View at: Publisher Site | Google Scholar
L. Hong, Y. Wang, W. Chen, and S. Yang, “MicroRNA-508 suppresses epithelial-mesenchymal transition, migration, and invasion of ovarian cancer cells through the MAPK1/ERK signaling pathway,” Journal of Cellular Biochemistry, vol. 119, no. 9, pp. 7431–7440, 2018.
View at: Publisher Site | Google Scholar
W. Wang, W. Hu, Y. Wang, J. Yang, and Z. Yue, “MicroRNA-508 is downregulated in clear cell renal cell carcinoma and targets ZEB1 to suppress cell proliferation and invasion,” Experimental and Therapeutic Medicine, vol. 17, no. 17, pp. 3814–3822, 2019.
View at: Publisher Site | Google Scholar
R. Li, L. Zhang, Z. Qin et al., “High LINC00536 expression promotes tumor progression and poor prognosis in bladder cancer,” Experimental Cell Research, vol. 378, no. 1, pp. 32–40, 2019.
View at: Publisher Site | Google Scholar
E. J. Nam, S. Kim, T. S. Lee et al., “Primary and recurrent ovarian high-grade serous carcinomas display similar microRNA expression patterns relative to those of normal ovarian tissue,” Oncotarget, vol. 7, no. 43, pp. 70524–70534, 2016.
View at: Publisher Site | Google Scholar
Y. Sun, L. Hu, H. Zheng et al., “MiR-506 inhibits multiple targets in the epithelial-to-mesenchymal transition network and is associated with good prognosis in epithelial ovarian cancer,” The Journal of Pathology, vol. 235, no. 1, pp. 25–36, 2015.
View at: Publisher Site | Google Scholar
G. Liu, Y. Sun, P. Ji et al., “MiR-506 suppresses proliferation and induces senescence by directly targeting the CDK4/6-FOXM1 axis in ovarian cancer,” The Journal of Pathology, vol. 233, no. 3, pp. 308–318, 2014.
View at: Publisher Site | Google Scholar
X. Mo, Y. Wu, L. Chen et al., “Global expression profiling of metabolic pathway-related lncRNAs in human gastric cancer and the identification of RP11-555H23.1 as a new diagnostic biomarker,” Journal of Clinical Laboratory Analysis, vol. 33, no. 2, Article ID e22692, 2019.
View at: Publisher Site | Google Scholar
W. Hou, Q. Tang, and F. Bi, “Comprehensive analysis of the aberrantly expressed profiles of lncRNAs, miRNAs and the regulation network of the associated ceRNAs in clear cell renal cell carcinoma,” Sheng Wu Yi Xue Gong Cheng Xue Za Zhi, vol. 36, no. 2, pp. 267–273, 2019.
View at: Publisher Site | Google Scholar
H. Chen, Y. Cai, Z. Zheng, Z. Chen, and W. Xie, “Clinical significance of LINC0047 in renal cell carcinoma and its effect on 786-O cells,” Lingnan Modern Clinics in Surgery, vol. 2, pp. 162–166, 2018.
View at: Google Scholar
D. L. Zhang, H. G. Zhou, J. Liu, J. Liu, and J. Mao, “Effect of LncRNA LINC00475 on biological behaviors of glioma cells,” Journal of Third Military Medical University, vol. 41, pp. 213–218, 2019.
View at: Google Scholar
D. J. H. van den Boomen, R. T. Timms, G. L. Grice et al., “TMEM129 is a Derlin-1 associated ERAD E3 ligase essential for virus-induced degradation of MHC-I,” Proceedings of the National Academy of Sciences, vol. 111, no. 31, pp. 11425–11430, 2014.
View at: Publisher Site | Google Scholar
N. Leffers, A. J. A. Lambeck, P. de Graeff et al., “Survival of ovarian cancer patients overexpressing the tumour antigen p53 is diminished in case of MHC class I down-regulation,” Gynecologic Oncology, vol. 110, no. 3, pp. 365–373, 2008.
View at: Publisher Site | Google Scholar
A. Mariani, C. Wang, A. L. Oberg et al., “Genes associated with bowel metastases in ovarian cancer,” Gynecologic Oncology, vol. 154, no. 3, pp. 495–504, 2019.
View at: Publisher Site | Google Scholar
G. Yoon, J. Y. Park, H. J. Kim et al., “ARID3A positivity correlated with favorable prognosis in patients with residual rectal cancer after Neoadjuvant chemoradiotherapy,” Anticancer Research, vol. 39, no. 6, pp. 2845–2853, 2019.
View at: Publisher Site | Google Scholar
Q. Ma, Y. Yang, D. Feng et al., “MAGI3 negatively regulates Wnt/β-catenin signaling and suppresses malignant phenotypes of glioma cells,” Oncotarget, vol. 6, no. 34, pp. 35851–35865, 2015.
View at: Publisher Site | Google Scholar
Y. Fan, Y. Shi, Z. Lin et al., “miR-9-5p suppresses malignant biological behaviors of human gastric cancer cells by negative regulation of TNFAIP8L3,” Digestive Diseases and Sciences, vol. 64, no. 10, pp. 2823–2829, 2019.
View at: Publisher Site | Google Scholar
K. Lian, C. Ma, C. Hao et al., “TIPE3 protein promotes breast cancer metastasis through activating AKT and NF-κB signaling pathways,” Oncotarget, vol. 8, no. 30, pp. 48889–48904, 2017.
View at: Publisher Site | Google Scholar
H. Zhang, Y. Wang, T. Chen et al., “Aberrant activation of Hedgehog signalling promotes cell migration and invasion via matrix metalloproteinase-7 in ovarian cancer cells,” Journal of Cancer, vol. 10, no. 4, pp. 990–1003, 2019.
View at: Publisher Site | Google Scholar
S. Sradhanjali and M. M. Reddy, “Inhibition of pyruvate dehydrogenase kinase as a therapeutic strategy against cancer,” Current Topics in Medicinal Chemistry, vol. 18, no. 6, pp. 444–453, 2018.
View at: Publisher Site | Google Scholar
J. D. Yang, S.-Y. Seol, S.-H. Leem et al., “Genes associated with recurrence of hepatocellular carcinoma: integrated analysis by gene expression and methylation profiling,” Journal of Korean Medical Science, vol. 26, no. 11, pp. 1428–1438, 2011.
View at: Publisher Site | Google Scholar
J. Kikuchi, M. Hori, H. Iha et al., “Soluble SLAMF7 promotes the growth of myeloma cells via homophilic interaction with surface SLAMF7,” Leukemia, vol. 34, no. 1, pp. 180–195, 2019.
View at: Publisher Site | Google Scholar
T. Tsujioka, A. Yokoi, Y. Itano et al., “Five-aza-2′-deoxycytidine-induced hypomethylation of cholesterol 25-hydroxylase gene is responsible for cell death of myelodysplasia/leukemia cells,” Scientific Reports, vol. 5, no. 1, p. 16709, 2015.
View at: Publisher Site | Google Scholar
O. S. Jeong, Y. C. Chae, H. Jung et al., “Long noncoding RNA linc00598 regulates CCND2 transcription and modulates the G1 checkpoint,” Scientific Reports, vol. 6, no. 1, p. 32172, 2016.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2020 Yu Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies