Systematic analyses of a novel circRNA-related miRNAs prognostic signature for Cervical Cancer

Abstract

Accumulating evidences shed light on the important roles of Circular RNAs (circRNAs) acting as competing endogenous RNAs (ceRNAs) in cervical cancer (CC) biology. The present study aimed to identify a novel circRNA-related prognostic signature for CC. The expression data and clinical information of CC were downloaded from the Gene Expression Omnibus (GEO) datasets to identify the differential circRNAs expression. Based on the targeted miRNA prediction, circRNA-related miRNAs were detected in training group and validation group of The Cancer Genome Atlas (TCGA) dataset to construct the novel prognostic signature of CC with least absolute shrinkage and selection operator (LASSO). Moreover, the Kaplan-Meier (K-M) analysis was applied to test the model. In the present study, three differentially expressed circRNAs (hsa_circ_0001498, hsa_circ_0066147, and hsa_circ_0006948) were identified in GSE102686 and GSE107472. Then, with the criteria 25 predicted miRNAs were analyzed in TCGA datasets to calculate the prognostic signature. Furthermore, we developed a six-miRNA signature (hsa-miR-217, hsa-miR-30b-3p, hsa-miR-136-5p, hsa-miR-185-3p, hsa-miR-501-5p and hsa-miR-658) based on their expression level and coefficients. We performed a Pearson correlation analysis to screen 47 mRNAs which are negatively regulated by these six miRNAs. Functional enrichment analysis indicated these mRNAs were mainly enriched in cancer-related biology, such as regulation of transcription, signal transduction, and cell cycle. The present study provides novel insight for better understanding of circRNA-related ceRNA network in CC and facilitates the identification of potential biomarkers for prognosis.

Keywords:
Circular RNAs; competing endogenous RNAs; miRNAs; least absolute shrinkage and selection operator; prognostic signature

Introduction

Among females, cervical cancer (CC) takes the fourth place for both commonly diagnosed cancer (570,000 cases) and the leading cause of cancer death (311,000 deaths) across 20 world regions in 2018 (Bray et al., 2018). There is a 10-fold variation in the incidence of CC by regions, which largely reflects human papillomavirus (HPV), human immunosuppression, parity, smoking, and oral contraceptive use (IARC Working Group on the Evaluation of Carcinogenic Risks to Humans, 2007; Tsikouras et al., 2016; Olusola et al., 2019). In spite of the development of early diagnosis and new targeted treatments, the five years overall survival (OS) rate of advanced-stage CC cases remains less than 40% (Wu et al., 2019). Therefore, it is a great importance to identify the novel biomarkers to predict CC cases with poor prognosis.

Recent studies reveal that noncoding RNAs (ncRNAs) play important biological roles in carcinogenesis and development of cancers such as cell proliferation, differentiation and apoptosis (Romano et al., 2017; Anastasiadou et al., 2018; Chan and Tay, 2018). Circular RNAs (circRNAs), reported as by-products of spliceosome-mediated splicing errors by previous published studies (Shen et al., 2015), and have been shown to exert novel crucial functions in the pathological progressions of various diseases, such as cardiovascular disease (Altesha et al., 2019), nervous system disease (Ma et al., 2020), digestive system disease (Sun et al.,, 2019). Because of the tissue-specific manner, circRNAs has drawn intense interest in their role in human diseases (Han B et al., 2018). Increasing number of studies suggest that aberrant expressions of circRNAs have appeared to participate in the development and progression of various kinds of cancers, including CC (Zhang et al., 2018; Verduci et al., 2019; Li et al., 2020). For instance, high circRNA_101996 expression level predicted poor outcomes and promoting tumor growth and metastasis by regulating miR-8075-TPX2 axis in CC (Song et al., 2019). CircRNA_0000285 was overexpressed in CC samples and CC cells and identified as an oncogene by suppressing expression level of FUS in CC (Chen et al., 2019).

The Competing Endogenous RNA (ceRNA) hypothesis was raised as a critical regulation mechanism between ncRNA and coding RNA (Salmena et al., 2011). Recent studies have shown circRNAs can interact with the microRNA (miRNA) as miRNA sponges by miRNA-binding sites (MREs) (Hansen et al., 2013), endogenously compete with mRNAs, and regulate protein translation (Abdelmohsen et al., 2017). Moreover, the miRNAs sponge is one of the most important and widely investigated functions in circRNAs. MiRNA is another type of ncRNA that plays a key role by acting as cancer suppressors or oncogenes in various diseases (Mishra et al., 2016). For example, by competing miR-326 to up-regulate the expression of ETS transcription factor ELK1 (ELK1), circ_0000515 enhanced proliferation and invasion but supressed apoptosis and autophagy of CC cells (Tang et al., 2019). CircCLK3 was shown to promote the growth and metastasis of CC in vitro and in vivo by regulating miR-320a/FoxM1 axis (Hong et al., 2019). Although various underlying function and mechanisms of circRNAs in CC have been characterized, the majority circRNAs contributing to biological properties of CC remain inconclusive.

The aim of the present study was to identify novel circRNA-related-ceRNA signatures for CC prognosis through data mining of the public database such asGene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) and The Cancer Genome Atlas (TCGA, http://cancergenome.nih.gov). By conducting a comprehensive reanalysis of circRNA expression profile, we constructed circRNA-related ceRNA, and identified a novel circRNA-related-miRNA signature as a new candidate indicator which has the potential to predict the over survival in CC cases. We present the following article in accordance with the REMARK reporting checklist.

Material and Methods

Identification of RNA in GEO and TCGA dataset

CircRNA High-throughput sequencing and microarray data from 10 pairs CC tissues and paired-paracancerous cervical tissues were downloaded from the two public microarray CC GEO datasets (GSE102686 and GSE107472). The basic information for these two microarray circRNA profiles is summarized in Table 1. Normalized CC miRNA and mRNA raw data (level 3) and clinical data from 307 CC tissues were obtained from the TCGA Data portal. The exclusion criteria were set as follows: i) histologic diagnosis is not CC; ii) suffering of other malignancy except CC; iii) patients samples without complete data for analysis; and iv) overall survival (OS) >5 years. Overall, a total of 283 CC patients were included in this study. These data were reanalysed by using R software.

Table 1 -
Basic information of the two microarray datasets from Gene Expression Omnibus.

Profile	RNA type	Platform	Organism	Experiment type	Sample size (T/N)	Region	Year
GSE102686	CircRNA	GPL19978	Homo sapiens	Non-coding RNA profiling by array	5/5	China	2017
GSE107472	CircRNA	GPL16791	Homo sapiens	Non-coding RNA profiling by high throughput sequencing	5/5	USA	2017

• Identification of differentially expressed circRNAs and target-regulation miRNA prediction

In this study, the differential expression of circRNAs were determined by fold-change and associated P-values in GSE102686 and GSE107472, respectively. Fold change represents the difference in expression of each circRNA between CC tissues and paired-paracancerous cervical tissues. The criteria of the differentially expressed circRNAs were fold changes>1.5 or <0.67 and P<0.05. Then, the intersection of differentially expressed circRNAs was selected from GSE102686 and GSE107472. In addition, both miRanda (http://www.microrna.org/micror-na/home.do) and PITA (http://genie.weizmann.ac.il/pubs/mir07/mir07_data.html) were used to predict the circRNAs-interact miRNAs.

Construction of the CC-specific circRNA-related miRNAs prognostic signature

All CC cases in TCGA datasets were evenly randomized into training group and validation group. Based on the above circRNAs-interact miRNAs, CC-specific circRNA-related miRNAs were selected in TCGA database. The expression of miRNAs should assurance more than 60% of the amount of expression greater than zero.

In the training group, Cox regression analysis with least absolute shrinkage and selection operator (LASSO) were performed to construct risk score to indicate the relative CC progression hazard for each case (Mao et al., 2019). LASSO minimizes the sum of squares of residuals whose absolute values are less than constant. LASSO regularization contains a parameter λ to limit the number of selected features, where the larger λ retains more features. Lambda.min is used to determine the λ value, and Lambda determines the minimum average cross validation error (Birnbaum et al., 2017). Afterward, a miRNA-related signature was determined by the sum value of miRNA expression profiles. According to the cutoff of risk score, cases were divided into High-risk group and Low-risk group. The hazard ratio (HR) and 95% confidence interval (CI) were assessed.

Predictive efficiency evaluation of the risk signature

Kaplan-Meier (K-M) plotter along with log-rank methods was further applied to compare the survival distributions between high‐risk and low‐risk groups (Sui et al., 2018). The time‐dependent receiver operating characteristic (ROC) curves was plotted to assess the predictive accuracy of the risk score for time-dependent outcomes analysis using the R package “survival-ROC” (Heagerty et al., 2000). Then univariate Cox regression analysis was applied to identify risk score-related miRNAs linked to OS.

Risk score-related miRNAs functional assessment

The target genes of risk score-related miRNAs were predicted by mirTarBase (http://mirtarbase.cuhk.edu.cn/php/index.php), miRanda, Targetscan (http://www.targetscan.-org/mamm_31/), and miRWalk (http://mirwalk.umm.uniheidelberg.de/). To identify the target genes with prognostic characteristics, the univariate Cox regression analysis was also used to explore the relationship between the mRNAs and the OS of CC cases in TCGA datasets. Furthermore, we performed a Pearson correlation analysis to screen the negatively regulated mRNAs. To explore the possible pathways and biological processes of negatively regulated target genes, KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways and GO (Gene Ontology) biological processes were analyzed through DAVID (the Database for Annotation, Visualization, and Integrated Discovery, http://david.abcc.ncifcrf.gov/) platform (Sui et al., 2016).

The circRNAs and risk score-related miRNAs verification

This part was approved by the Ethics Committee of the Sir Run Run Shaw Hospital (20211018093606877). CC and normal control samples were obtained from biobank of Sir Run Run Shaw Hospital ZheJiang University School of Medicine. All patients had signed informed consent for donating their samples to biobank of Sir Run Run Shaw Hospital ZheJiang University School of Medicine. Total RNA was extraction from nine pairs CC tissues and paired-paracancerous cervical tissues with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacture's protocol. U6 was used as an internal normalized reference to confirm reliability and validity of miRNAs. Reverse transcription reactions using A214 reverse transcription system kit (Genstar, Beijing, China) was conducted in two steps according to the manufacturer's protocol. Real-time PCR was performed to detect the expression level of the candidate miRNAs with the StepOne Plus™PCR System (Applied biosystems, Foster City, CA, USA). qRT-PCR was then performed with 2×RealStar Green Fast Mixture with ROX (Genstar, A303, Beijing, China) according to the manufacturer's protocol. The primer sequences were as follows: miR-185-3p, forward, 5′‐GATCACACTCTTGTGGTAGTTGC‐3′, reverse 5′‐CTCTTCCTTGCTCGTTGTTGGTAT‐3′, miR‐501-5p, forward, 5′‐TGCGCAATCCTTTGTCCCTGGG‐3′, reverse, 5′‐CCAGTGCAGGGTCCGAGGTATT‐3′, and U6, forward 5′‐CTCGCTTCGGCAGCACA‐3′, reverse 5′‐AACGCTTCACGAATTTGCGT‐3′. The PCR reaction components were 0.4 µL of cDNA, 10 µl of 2×RealStar Green Fast Mixture, 0.8 µl (1 µl/pmol) PCR primers, and 8.8 µL RNase-free water. The reaction was performed at 95˚C for 2 min, followed by 40 cycles of 95˚C for 15 s, 60˚C for 30 s and 72˚C for 30 s. A dissociation curve was analyzed from 60 to 95˚C. The Ct-value for each sample was calculated with the ∆∆Ct method, and the results of fold change were calculated by 2^-∆∆Ct, where ∆∆Ct = (Ct_miRNAs - Ct_U6 )_tumor - (Ct_miRNAs - Ct_U6 )_{paracancerous cervical tissues}. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Statistical analysis

Statistical analyses were performed using IBM SPSS 21.0 software (Chicago, USA). The results were expressed as mean±standard deviation (SD). Statistical significance was determined using the Student's t test or One-way Analysis of Variance (ANOVA) at the probability of P<0.05.

Results

Identification of differentially expressed circRNAs

With the criteria of fold changes>1.5 or <0.67 and P<0.05, differentially expressed circRNAs profile consisting of 507 circRNAs, with 253 downregulated and 254 upregulated circRNAs, was identified from GSE102686 (Figure 1A). In GSE107472, 73 downregulated and 160 upregulated circRNAs were identified (Figure 1B). To further consolidate the data reliability, we selected three differentially expressed circRNAs (hsa_circ_0001498 downregulated, and hsa_circ_0066147 and hsa_circ_0006948 upregulated) from the intersection of the above two groups (Figure 1C).

Figure 1-
Differentially expressed circRNAs in GEO. (A) Differentially expressed circRNAs from GSE107472. (B) Differentially expressed circRNAs from GSE102686. (C) Intersection of differentially expressed circRNAs.

circRNAs target-regulation miRNA prediction

Based on the above three circRNAs (hsa_circ_0001498, hsa_circ_0066147 and hsa_circ_0006948), a total of 244 and 534 miRNAs were predicted by miRanda and PITA, respectively. After merging these two results, 146 miRNAs were screened (^{Table S1}). Furthermore, with the criteria of more than 60% of the amount of expression greater than zero, only 25 miRNAs left in TCGA datasets.

Identification of a miRNA‐related prognostic signature

283 CC cases from TCGA datasets were randomly divided into training group (n=142) and validation group (n=141). To identify the novel prognostic signature for CC cases, LASSO Cox regression model was implemented to analyze the data of miRNAs expression profiles in the training group. A batch of six miRNAs, hsa-miR-217, hsa-miR-30b-3p, hsa-miR-136-5p, hsa-miR-185-3p, hsa-miR-501-5p and hsa-miR-658, was selected as the most accurate prognostic predictor in the training group (Figure 2A). Further, we developed a six-miRNA signature based on the expression level and coefficients of these miRNAs. Risk score = 0.1944*hsa-miR-217+(-0.029*hsa-miR-30b-3p)+(-0.1975*hsa-miR-136-5p)+(-0.0179*hsa-miR-185-3p)+(-0.1869*hsa-miR-501-5p)+(-0.8713*hsa-miR-658) (Figure 2B).

Figure 2-
Coefficients of selected features are shown by lambda parameter (A) LASSO coefficient profiles of the miRNAs associated with the overall survival of cervical cancer. (B) Plots of the cross-validation error rates.

The risk score for each case in all CC cases, training group and validation group was calculated and plotted. As shown in Figure 3, the value of risk score in the outcome of alive cases was significantly lower than dead cases in all CC cases, training group and validation group (P<0.05). The robustness of the miRNA signature was detected by evaluating their ability to classify low-risk group and high-risk group in training group and validation group. Firstly, the miRNA signature-based risk score for each CC case was calculated. Then, CC cases were divided into low-risk group or high-risk group according to the median value. K-M curves were plotted along with log rank p-test, to compare the OS of the two groups. According to the results, significant differences in K-M survival analysis were similarly observed in low-risk group and high-risk group separated by the miRNA signature in training group (Figure 4A) and validation group (Figure 4B). The time-dependent ROC curves were used to determine the sensitivity and specificity of the predictive model. As shown in Figure 4C, predict accuracy of the novel miRNA signature in training group and validation group. In addition, univariate Cox proportional hazards regression showed that two (hsa-miR-185-3p and hsa-miR-501-5p) of these above six miRNAs in the TCGA datasets were identified with a significant prognostic value (Figure 5).

Figure 3-
circRNA risk score analysis of 283 cervical cancer patients in TCGA cohort.(A) Risk score of circRNA signature in total cervical cancer cases. The value of risk score in the outcome of alive cases was significantly lower than dead cases (P<0.05). (B) Risk score of circRNA signature in training group. The value of risk score in the outcome of alive cases was significantly lower than dead cases (P<0.05). (C) Risk score of circRNA signature in validation group. The value of risk score in the outcome of alive cases was significantly lower than dead cases (P<0.05).283 cervical cases were randomly divided into training group (n=142) and validation group (n=141).

Figure 4-
The six-miRNA prognostic signature of cervical cancer for outcome in cervical cancer patients in TCGA cohort. (A) The Kaplan-Meier curves for cervical cancer risk groups obtained from the TCGA cohort (training group, n=142) divided by the median cutoff point. Patients with high scores had poor outcome in terms of OS (P<0.05). (B) The Kaplan-Meier curves for cervical cancer risk groups obtained from the TCGA cohort (validation group, n=141) divided by the median cutoff point. Patients with high scores had poor outcome in terms of OS (P<0.05).(C) The risk score shown by the time-dependent ROC curve for predicting 5-year survival in training group and validation group. 283 cervical cases were randomly divided into training group (n=142) and validation group (n=141).

Figure 5 -
The Kaplan-Meier curves for miRNA obtained from 283 cervical cancer patients in TCGA cohort. (A) The Kaplan-Meier curves for miR-185-3p obtained from the TCGA cervical cancer cohort (n=283) divided by the median cutoff point. Patients with high expression of miR-185-3p had better outcome in terms of OS (P<0.05). (B) The Kaplan-Meier curves for miR-501-5p obtained from the TCGA cervical cancer cohort (n=283) divided by the median cutoff point. Patients with high expression of miR-501-5p had better outcome in terms of OS (P<0.05).

Prediction of miRNA targets

To inquire into the KEGG pathways and GO biological processes relevant to the six-prognostic signature-related miRNAs, we first predicted the targeted correlation between miRNA and mRNAs with mirTarBase, miRanda, Targetscan, and miRWalk. 1226, 12690, 7083, and 8356 targeted mRNAs were identified respectively (^{Figure S1}). The results of univariate Cox regression between these above mRNAs and the OS of CC cases in TCGA datasets were shown in ^{Figure S2}. Then, based on negative regulation theory (Fabian et al., 2010), we performed a Pearson correlation analysis to screen 47 negatively regulated mRNAs (^{Table S2}). MSMO1, OSBPL2, TANGO2, and WHAMM in these 47 mRNAs were significantly associated with OS (P<0.05) (Figure 6).

Figure 6-
The Kaplan-Meier test of negatively regulated mRNAs (MSMO1, OSBPL2, TANGO2 and WHAMM) for the OS.

Functional enrichment analysis for miRNAs signature

Furthermore, these above 47 negatively regulated mRNAs were utilized for KEGG pathways and GO biological processes enrichment analysis by DAVID. Pathway enrichment analysis indicated three significantly enriched pathways were Oxytocin signaling pathway, Protein processing in endoplasmic reticulum, and Proteoglycans in cancer (Figure 7A). For GO analysis, when considering CC, the mRNAs were enriched in regulation of transcription, signal transduction, regulation of immune response, and cell cycle (Figure 7B).

Figure 7 -
GO and KEGG term analysis of potential genes related to negatively regulated mRNAs. (A)Three top enrichment of KEGG pathways. (B) Twenty top enrichment of GO terms. The Rich factor shows the degree of enrichment, which was calculated by the formula: (the number of selected genes in a term/total number of selected genes)/(the total number of genes in a term of the database/ the total number of genes in the database). The node size represents the number of selected genes, and color represents the P-value of the enrichment analysis.

The expression of miRNA verification.

Finally, hsa-miR-185-3p was downregulated expressed in CC tissues, while hsa-miR-501-5p was upregulated expressed in CC tissues detected by qRT-PCR (9 pairs of CC and adjacent tissues were tested) (Figure 8).

Figure 8-
Results of qRT-PCR validation of miR-185-3p (A) and miR-501-5p (B) in cervical cancer tissues and paired-paracancerous tissues (upper panel: histogram, lower panel: violin plot).

Discussion

In the present study, the circRNAs expression profiles of CC cases downloaded from GEO datasets were reanalyzed and three CC-specific circRNAs were identified. According to the targeted-miRNAs prediction of these three circRNAs, CC cases in TCGA datasets were randomly separated to training group and validation group to discover a prognostic six‐miRNA signature related to OS by LASSO Cox regression model. LASSO is considered to be a novel high‐dimensional variable selection method for regression (Tibshirani, 1997). LASSO is extremely important to help select parsimonious variables in the application of increasing the efficiency of CC prognosis prediction. Furthermore, LASSO can also improve the prediction accuracy and enhance numerical stability and model interpretability. Risk score was calculated with the regression coefficients that generated by the Cox regression analysis. Based on the prognostic signature, the risk score was calculated for each CC case and the training group was divided into low‐risk group and high‐risk group. Then prognostic value of the six‐miRNA signature was verified in the validation group. K-M analysis showed that distinct separation between the survival curves of the high‐risk group and low‐risk group that classified by the same criteria in both training group and validation group was observed. These results indicated this miRNA signature in CC cases is high reproducible.

To further assess the predictive efficiency based on risk score, time‐dependent ROC curves was introduced into our study. The results showed that the novel prognostic signature has a better predictive accuracy. Moreover, the risk score based on these six-miRNAs could be used as a new indicator for the prognosis of CC cases. Generally, miRNAs act as the ceRNAs which regulate gene expression through mRNA processing and degradation (Wang et al., 2017). Their related negatively regulated mRNAs were predicted in various software and TCGA datasets to analyze the possible biological processes. The functional enrichment analysis suggested that these above mRNAs were enriched in cancer‐related signaling pathways such as proteoglycans in cancer and protein processing in endoplasmic reticulum, and biological process such as signal transduction, regulation of immune response, and cell cycle. The results of functional enrichment clarified the possible effects of the CC-specific circRNAs related miRNAs in the prognostic signature.

It has been demonstrated that non-coding RNAs play an important role in the progression of cervical cancer. For instance, for 158 cervical specimens, including 38 normal, 52 cervical intraepithelial neoplasia and 68 CC tissues, the result showed a remarkable increase of miR-25, miR-92a, and miR-378 with lesion progression but no obvious change of miR-22, miR-29a, and miR-100 among the HPV-infected tissues (Wang et al., 2014). It was also reported that miR-18a targets the tumor suppressor STK4/MST1, which is necessary for HPV-positive cervical cancer transformation (Morgan et al., 2020). Long noncoding RNAs (lncRNAs) are also involved in the regulation of cellular processes that are dysregulated during carcinogenesis (Anastasiadou et al., 2018). For example, the long non-coding RNA DINO is acutely expressed in human papillomavirus-positive cervical cancer cells and reactivates the dormant TP53 tumor suppressor through ATM/CHK2 signaling (Sharma and Munger, 2020a). In addition, the expression of long non-coding RNA DINO is mediated by KDM6A, which leads to the stability of TP53 tumor suppressor in cells expressing human papillomavirus 16 E7 (Sharma and Munger, 2020b).

Previous studies have demonstrated that miR-185-3p is downregulated in CC tissues and regulated the proliferation and self-renewal ability of CC stem cells by repressing FOXP3 expression in CC cells (Zhang et al., 2021). Besides, unexpressed miR‐185 promoted apoptosis, suppressed cell migration and invasion, and related with better OS (Gao et al., 2020). Upregulated miR-501-5p was detected in gastric cancer tissues and cell lines with poor OS in gastric cancer patients by targeted and suppressed multiple repressors of the WNT/β-catenin signaling (Fan et al., 2016). A few studies have reported the role of miR-217 in CC. Zhu et al. (2019) found miR-217 was downregulated in CC tissues and CC cells. Upregulated miR-217 suppressed the viability, migration and invasion in CC cells (Dong et al., 2018; Yin and Ren, 2019; Zhu et al., 2019). MiR-30b-3p has been reported to be a potential targeted regulatory of GPNMB to regulate the prognosis of gastric cancer patients (Yao et al., 2020). MiR-136 has been shown to be downregulated in CC tissues (Han M-S et al., 2018), and regulated by lncRNA FOXP4-AS1 to involve the CC progression (Zhao et al., 2020). Xu et al. (2020) reported circNFIC suppresses breast cancer progression by sponging miR-658. It has been demonstrated that circRNAs act as competitive endogenous RNAs (ceRNAs) and play an important role in the progression of cervical cancer. CircRNA8924 is highly expressed in CC tissue and can be considered a competitive endogenous RNA of the miR-518d-5p/519-5p family to promote the malignant biological behavior of CC cells (Liu et al., 2018). circ_0067934 and circ_0005576 could promote CC progression via miR-545/EIF3C axis and circ_0005576/miR-153-3p/KIF20A axis, respectively (Hu et al., 2019; Ma et al., 2019). However, there is no report on the prognosis of these six miRNAs with CC. Thus, the six-miRNAs prognostic signature might be a novel prognostic biomarker in CC treatment and survival evaluation.

Although the results of the present study might have remarkable prognostic implications, some limitations should also be considered. First, a longer follow-up time should be verified in this study. Second, in addition to the data in the GEO and TCGA datasets, other experimental methods should apply to validate the results in the present study. Third, the function of the three circRNAs and six miRNAs should be elucidated in further studies.

In conclusion, the present study identified a CC-specific circRNA-related six-miRNA signature could be used as a potential outcome predictor for CC cases via re-analyzing the genome-wide RNA expression profile datasets from GEO and TCGA database with a ceRNA network. Future functional studies are required to elucidate the roles of these circRNAs and miRNAs in CC.

Acknowledgements

All people collaborated in the collection and interpretation of the data and contributed to the manuscript. The results present in the manuscript are mainly based on bioinformatics analysis and scientific computation. Limited experiments were performed based on proprietary conditions.

References

• Abdelmohsen K, Panda AC, Munk R, Grammatikakis I, Dudekula DB, De S, Kim J, Noh JH, Kim KM, Martindale JL et al (2017) Identification of HuR target circular RNAs uncovers suppression of PABPN1 translation by CircPABPN1. RNA Biol 14:361-369.
• Altesha M-A, Ni T, Khan A, Liu K and Zheng X (2019) Circular RNA in cardiovascular disease. J Cell Physiol 234:5588-5600.
• Anastasiadou E, Jacob LS and Slack FJ (2018) Non-coding RNA networks in cancer. Nat Rev Cancer 18:5-18.
• Birnbaum DJ, Finetti P, Lopresti A, Gilabert M, Poizat F, Raoul J-L, Delpero J-R, Moutardier V, Birnbaum D, Mamessier E et al (2017) A 25-gene classifier predicts overall survival in resectable pancreatic cancer. BMC Med 15:170.
• Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA and Jemal A (2018) Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 68:394-424.
• Chan JJ and Tay Y (2018) Noncoding RNA:RNA regulatory networks in cancer. Int J Mol Sci 19:1310.
• Chen R-X, Liu H-L, Yang L-L, Kang F-H, Xin L-P, Huang L-R, Guo Q-F and Wang Y-L (2019) Circular RNA circRNA_0000285 promotes cervical cancer development by regulating FUS. Eur Rev Med Pharmacol Sci 23:8771-8778.
• Dong J, Wang M, Ni D, Zhang L, Wang W, Cui X, Fu S and Yao S (2018) MicroRNA-217 functions as a tumor suppressor in cervical cancer cells through targeting Rho-associated protein kinase 1. Oncol Lett 16:5535-5542.
• Fabian MR, Sonenberg N and Filipowicz W (2010) Regulation of mRNA translation and stability by microRNAs. Annu Rev Biochem 79:351-379.
• Fan D, Ren B, Yang X, Liu J and Zhang Z (2016) Upregulation of miR-501-5p activates the wnt/beta-catenin signaling pathway and enhances stem cell-like phenotype in gastric cancer. J Exp Clin Cancer Res 35:177.
• Gao P, Wang Z, Hu Z, Jiao X and Yao Y (2020) Circular RNA circ_0074027 indicates a poor prognosis for NSCLC patients and modulates cell proliferation, apoptosis, and invasion via miR-185-3p mediated BRD4/MADD activation. J Cell Biochem 121:2632-2642.
• Han B, Chao J and Yao H (2018) Circular RNA and its mechanisms in disease: From the bench to the clinic. Pharmacol Ther 187:31-44.
• Han M-S, Lee JM, Kim S-N, Kim JH and Kim H-S (2018) Human papillomavirus 16 oncoproteins downregulate the expression of miR-148a-3p, miR-190a-5p, and miR-199b-5p in cervical cancer. Biomed Res Int 2018:1942867.
• Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK and Kjems J (2013) Natural RNA circles function as efficient microRNA sponges. Nature 495:384-388.
• Heagerty PJ, Lumley T and Pepe MS (2000) Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics 56:337-344.
• Hong H, Zhu H, Zhao S, Wang K, Zhang N, Tian Y, Li Y, Wang Y, Lv X, Wei T et al (2019) The novel circCLK3/miR-320a/FoxM1 axis promotes cervical cancer progression. Cell Death Dis 10:950.
• Hu C, Wang Y, Li A, Zhang J, Xue F and Zhu L (2019) Overexpressed circ_0067934 acts as an oncogene to facilitate cervical cancer progression via the miR-545/EIF3C axis. J Cell Physiol 234:9225-9232.
• IARC Working Group on the Evaluation of Carcinogenic Risks to Humans (2007) Human papillomaviruses. IARC Monogr Eval Carcinog Risks Hum 90:1-636.
• Li R, Jiang J, Shi H, Qian H, Zhang X and Xu W (2020) CircRNA: A rising star in gastric cancer. Cell Mol Life Sci 77:1661-1680.
• Liu J, Wang D, Long Z, Liu J and Li W (2018) CircRNA8924 promotes cervical cancer cell proliferation, migration and invasion by competitively binding to MiR-518d-5p /519-5p family and modulating the expression of CBX8. Cell Physiol Biochem 48:173-184.
• Ma H, Tian T, Liu X, Xia M, Chen C, Mai L, Xie S and Yu L (2019) Upregulated circ_0005576 facilitates cervical cancer progression via the miR-153/KIF20A axis. Biomed Pharmacother 118:109311.
• Ma N, Zhang W and Wan J (2020) Research progress on circRNA in nervous system diseases. Curr Alzheimer Res 17:687-697.
• Mao Y, Fu Z, Zhang Y, Dong L, Zhang Y, Zhang Q, Li X and Wang C (2019) A six-microRNA risk score model predicts prognosis in esophageal squamous cell carcinoma. J Cell Physiol 234:6810-6819.
• Mishra S, Yadav T and Rani V (2016) Exploring miRNA based approaches in cancer diagnostics and therapeutics. Crit Rev Oncol Hematol 98:12-23.
• Morgan EL, Patterson MR, Ryder EL, Lee SY, Wasson CW, Harper KL, Li Y, Griffin S, Blair GE, Whitehouse A et al (2020) MicroRNA-18a targeting of the STK4/MST1 tumour suppressor is necessary for transformation in HPV positive cervical cancer. PLoS Pathog 16:e1008624.
• Olusola P, Banerjee HN, Philley JV and Dasgupta S (2019) Human papilloma virus-associated cervical cancer and health disparities. Cells 8:662.
• Romano G, Veneziano D, Acunzo M and Croce CM (2017) Small non-coding RNA and cancer. Carcinogenesis 38:485-491.
• Salmena L, Poliseno L, Tay Y, Kats L and Pandolfi PP (2011) A ceRNA hypothesis: The rosetta stone of a hidden RNA language? Cell 146:353-358.
• Sharma S and Munger K (2020a) Expression of the long noncoding RNA DINO in human papillomavirus-positive cervical cancer cells reactivates the dormant TP53 tumor suppressor through ATM/CHK2 signaling. mBio 11:e01190-20.
• Sharma S and Munger K (2020b) KDM6A-mediated expression of the long noncoding RNA DINO causes TP53 tumor suppressor stabilization in human papillomavirus 16 E7-expressing cells. J Virol 94:e02178-19.
• Shen T, Han M, Wei G and Ni T (2015) An intriguing RNA species--perspectives of circularized RNA. Protein Cell 6:871-880.
• Song T, Xu A, Zhang Z, Gao F, Zhao L, Chen X, Gao J and Kong X (2019) CircRNA hsa_circRNA_101996 increases cervical cancer proliferation and invasion through activating TPX2 expression by restraining miR-8075. J Cell Physiol 234:14296-14305.
• Sui J, Li Y-H, Zhang Y-Q, Li C-Y, Shen X, Yao W-Z, Peng H, Hong W-W, Yin L-H, Pu Y-P et al (2016) Integrated analysis of long non-coding RNAassociated ceRNA network reveals potential lncRNA biomarkers in human lung adenocarcinoma. Int J Oncol 49:2023-2036.
• Sui J, Miao Y, Han J, Nan H, Shen B, Zhang X, Zhang Y, Wu Y, Wu W, Liu T et al (2018) Systematic analyses of a novel lncRNA-associated signature as the prognostic biomarker for hepatocellular carcinoma. Cancer Med 7:3240-3256.
• Sun Z, Chen C, Su Y, Wang W, Yang S, Zhou Q, Wang G, Li Z, Song J, Zhang Z et al (2019) Regulatory mechanisms and clinical perspectives of circRNA in digestive system neoplasms. J Cancer 10:2885-2891.
• Tang Q, Chen Z, Zhao L and xu H (2019) Circular RNA hsa_circ_0000515 acts as a miR-326 sponge to promote cervical cancer progression through up-regulation of ELK1. Aging (Albany NY) 11:9982-9999.
• Tibshirani R (1997) The Lasso method for variable selection in the Cox model. Stat Med 16:385-395.
• Tsikouras P, Zervoudis S, Manav B, Tomara E, Iatrakis G, Romanidis C, Bothou A and Galazios G (2016) Cervical cancer: Screening, diagnosis and staging. J BUON 21:320-325.
• Verduci L, Strano S, Yarden Y and Blandino G (2019) The circRNA-microRNA code: Emerging implications for cancer diagnosis and treatment. Mol Oncol 13:669-680.
• Wang D, Xin L, Lin J-H, Liao Z, Ji J-T, Du T-T, Jiang F, Li Z-S and Hu L-H (2017) Identifying miRNA-mRNA regulation network of chronic pancreatitis based on the significant functional expression. Medicine (Baltimore) 96:e6668.
• Wang X, Wang H-K, Li Y, Hafner M, Banerjee NS, Tang S, Briskin D, Meyers C, Chow LT, Xie X et al (2014) microRNAs are biomarkers of oncogenic human papillomavirus infections. Proc Natl Acad Sci U S A 111:4262-4267.
• Wu W, Sui J, Liu T, Yang S, Xu S, Zhang M, Huang S, Yin L, Pu Y and Liang G (2019) Integrated analysis of two-lncRNA signature as a potential prognostic biomarker in cervical cancer: A study based on public database. PeerJ 7:e6761.
• Xu G, Ye D, Zhao Q, He R, Ma W, Li Y, Tang S, Zhou Z, Li X and Zhang Z (2020) circNFIC suppresses breast cancer progression by sponging miR-658. J Cancer 11:4222-4229.
• Yao K, Wei L, Zhang J, Wang C, Wang C, Qin C and Li S (2020) Prognostic values of GPNMB identified by mining TCGA database and STAD microenvironment. Aging (Albany NY) 12:16238-16254.
• Yin Z and Ren W (2019) MicroRNA-217 acts as a tumor suppressor and correlates with the chemoresistance of cervical carcinoma to cisplatin. Onco Targets Ther 12:759-771.
• Zhang H-D, Jiang L-H, Sun D-W, Hou J-C and Ji Z-L (2018) CircRNA: A novel type of biomarker for cancer. Breast Cancer 25:1-7.
• Zhang L, Guo C, Ji T and Chen X (2021) SOX2 Regulates lncRNA CCAT1/MicroRNA-185-3p/FOXP3 axis to affect the proliferation and self-renewal of cervical cancer stem cells. Nanoscale Res Lett 16:2.
• Zhao J, Yang T and Li L (2020) LncRNA FOXP4-AS1 is involved in cervical cancer progression via regulating miR-136-5p/CBX4 axis. Onco Targets Ther 13:2347-2355.
• Zhu L, Yang S and Wang J (2019) miR-217 inhibits the migration and invasion of HeLa cells through modulating MAPK1. Int J Mol Med 44:1824-1832.

Supplementary material

The following online material are available for this article:

Table S1 - The intersected miRNAs predicted by targeted circRNAs. Table S2 - miRNAs targeting CC-specific mRNAs. Figure S1 - Predicted intersection of targeted mRNA with mirTarBase, miRanda, Targetscan, and miRWalk. Figure S2 - The Kaplan-Meier test of six-miRNA targeted mRNAs.

Publication Dates

• Publication in this collection
  24 June 2022
• Date of issue
  2022

History

• Received
  22 Dec 2021
• Accepted
  17 May 2022