A total of 73 specimens were obtained from patients who underwent surgery at Beijing Tiantan Hospital (Beijing, China) between October 2007 and July 2014. The patients included 34 males and 39 females, with a median age of 55 years (range 25–73 years) and a median follow-up of 78 months (range 36–121 months). Tumors with Hardy-Wilson classification grade IV and/or Knosp classification grades III and IV (32,33) were defined as invasive. Table I includes the clinical and pathological features of the adenomas. All tissue samples were immediately snap-frozen and stored in liquid nitrogen or at −80°C. The tumour type was determined based on pathological diagnosis and the preoperative clinical and biochemical examination results. This study was ethically approved by the Tiantan Hospital Ethics Committee of Beijing Tiantan Hospital, and written informed consent was obtained from each patient prior to the study.

Total RNA from a tumour was isolated and purified using a phenol-free mirVana™ miRNA Isolation kit (cat. no. AM1561; Ambion; Thermo Fisher Scientific, Inc., Waltham, MA, USA), according to the manufacturer's protocol, and the purity and quantity were measured using a Nanodrop spectrophotometer. Subsequently, these 73 total RNA samples (1.8 < Optical Density 260 < 2.1) were used to produce fluorescence-labelled cRNA targets for the SBC human ceRNA array V1.0 (4×180 K). An Agilent Microarray scanner (Agilent Technologies, Inc., Santa Clara, CA, US) was used to scan the slides and hybridize the labelled cRNA targets. Feature Extraction software 10.7 (Agilent Technologies, Inc.) was used to extract the data. The R 3.2.3 software (algorithm ‘mean’ and ‘var’; www.r-project.org/) was used to normalize the raw data and for data processing.

The present study used the R 3.2.3 (www.r-project.org/) program algorithm ‘sample’ to randomly divide the dataset into training (n=37) and test groups (n=36). The association between the continuous expression level of each circRNA and patient PFS in the training set was determined using univariable Cox regression analysis: A gene with a univariate Cox coefficient >0 was considered to be a risk factor, <0 was considered to be a protective factor and -log₁₀ P>1.30 (-log₁₀0.05) was regarded as significant.

Next, a model for selecting circRNA signatures with prognostic indicators was developed as follows:

Risk score = ∑^Ni = 1 (Explg*Coef) (30,34).

Where N is the number of predictive circRNAs, Explg indicates the expression of circRNAs and Coef is the evaluated regression coefficient of the circRNAs determined from the univariable Cox regression analysis.

circRNAs were screened using a random survival forest-variable hunting (RSFVH) algorithm, and then a total of 9 circRNAs were selected (30) as fewer circRNAs provide a greater value. All 9 of the circRNAs were analysed. The total of all 9 node permutations and combinations was 2⁹−1=511 and the sensitivity and specificity of the survival prediction of the risk score of the 2⁹−1=511 signatures in the training dataset were compared with time-dependent receiver operating characteristic (ROC) curves, and the area under the curve (AUROC) values were obtained. In the training test, the marginal value was obtained from the median risk score. According to the median risk score, the patients with NFPA were separated into a high-risk group and a low-risk group in each dataset.

The equality of survival distributions in different groups for the training and test NFPA cohorts was assessed using Kaplan-Meier survival analyses, and two-sided log-rank tests were used to evaluate the statistical significance. Furthermore, the association with clinical features was analysed via χ² tests. The risk grade and other clinical features in the available data when independent were tested using multivariate Cox regression analysis and data stratification analysis. P<0.05 was considered to indicate a statistically significant difference. The R 3.2.3 program was used for the analysis of all data and to generate plots; and the pROC, survival, ggplot2 and random forest SRC packages were downloaded from Bioconductor (www.r-project.org) (35,36). The analysis sequence used for development of a risk scoring model and to verify the validity of signatures is presented in Fig. 1.

Pearson's correlation coefficients were used to determine the co-expression association between prognostic circRNAs and protein-coding genes. Gene ontology (GO) (37,38) and Kyoto Encyclopaedia Gene and Genome (KEGG) (39–41) enrichment analyses using the Cytoscape plugin ClueGo (version 3.2.3) were performed to determine the co-expressed protein-coding genes and prognostic circRNAs for predicting the biological function of prognostic circRNAs (42). Cytoscape is a functional annotation tool used for assessing the over-representation of a category of interest gene concentration. Enrichment analyses were performed using functional annotation maps and functional annotation clustering options and were restricted to the KEGG pathway and GO terminology in the ‘biological process’ category. Functional annotations with P<0.05 were considered to indicate a statistically significant difference.

The present study examined the pituitary tumour tissues of 73 patients using an Agilent microarray. From this analysis, a total of 88,750 circRNAs were identified. Next, the patients were randomly divided into two groups: A training set (n=37) and a test set (n=36). The training set was used to select the circRNA prognostic signature. The test set was used to evaluate the discriminative power of the circRNA signature, which was obtained independently. In this model selection section, all analyses were based on the training set.

First, univariate Cox proportional hazards regression analysis of circRNA expression using recurrence as the dependent variable was performed. A total of 7,481 circRNAs (data not shown) were determined to be correlated (P<0.05) with patient recurrence/progression and a volcano plot was used to present the univariate Cox results (Fig. 2A). Next, a RSFVH algorithm was used to select 9 circRNAs that were most closely associated with progression/relapse in the set of 7,481 circRNAs, according to the replacement vital score of the RSFVH algorithm (Fig. 2B).

The application of time-dependent ROC curves and the comparison of the sensitivity and specificity of the risk score survival prediction may result in better predictive characteristics. The risk score is based on the risk score model of 2⁹−1=511 combinations of the circRNA training set (data not shown) for each circRNA of each patient in the training dataset. All risk scores for the circRNA signatures were calculated using this method. Then, the circRNA combination of hsa_circ_0000066 and hsa_circ_0069707, which had the largest AUROC, was selected (Fig. 2C and D; Table II). A risk score for the combination of hsa_circ_0000066 and hsa_circ_0069707 was obtained as follows: Risk score=(2.38 × expression value of RP11-702B10.1) + (8.27 × expression value of hsa_circ_0069707). The AUROC for the circRNA signature in the model was 0.87 in the training set. Thus, the optimal survival prediction value was obtained based on the expression-based risk score.

Subsequently, the predictive power of the circRNA signature was assessed. The test set (n=36) was utilized for an independent evaluation, and the training set (n=37) was evaluated in parallel for comparison.

In the training set, based on the median value of the risk score of the circRNA signature being used as the cutoff value, the patients were separated into a high-risk group (n=18) and a low-risk group (n=19). The PFS time was significantly shorter in the high-risk patients compared with the low-risk patients (median survival: 56 months vs. >80.2 months; P<0.001; Fig. 3A). The rate of no recurrence or progression in high-risk patients was <1% and that in low-risk patients was >90%.

In the independent test set, the patients were separated into two groups (n=18). Verification of the recurrence prediction ability of the circRNA signature was performed using Kaplan-Meier curves for the high-risk and low-risk groups in the test dataset; the data were plotted and are presented in Fig. 3B, revealing a significant difference in PFS between the two groups (median PFS time: 70.2 months vs. >60 months; log-rank sum test, P=0.047). The rate of no recurrence/progression was ~16.6% in the high-risk group and 50% in the low-risk group.

The clinical significance of the circRNA signature in NFPAs was obtained by associating the signature with a range of clinical pathological parameters from the dataset (n=73). As presented in Table III, except for age, vision and visual field disorders, and headache, there was no significant association between the circRNA markers and clinical pathological variables, including sex and histological grade (Table III).

To assess whether the ability of the circRNA markers in predicting the progression of recurrence was independent of other clinical features, multivariable Cox regression analysis was performed using the NFPA dataset and a signature-based risk score based on two circRNAs and other clinical features. The analysis indicated that the circRNA marker risk score was independent of the clinical characteristics in predicting survival, but that the signature was significantly associated with survival (high-risk group vs. low-risk group; hazard ratio=10.07, 95% confidence interval, 2.92–34.77; P<0.001, n=73; Table IV).

In order to compare the sensitivity and specificity of age and the circRNA-specific survival predictions, ROC analysis was performed as a larger AUROC usually indicates a better predictive model. In the dataset (n=37), the prognostic power of the circRNA marker was significantly stronger compared with age for the training and test groups (AUROC=0.87 and 0.69 for signature, respectively, vs. AUROC=0.65 and 0.60 for age, respectively; Fig. 3C and D). These data further confirm that the circRNA marker is a novel predictive marker with a high accuracy and thus has major clinical value.

Pearson's correlation coefficient was used for 73 patients to calculate the co-expression association between the two circRNAs (hsa_circ_0000066 and hsa_circ_0069707) and protein-coding genes to determine the potential biological effects of these molecular markers. Of the two circRNAs, at least one may be associated with the expression levels of 623 protein-coding genes, and the interaction network was determined with Cytoscape (Pearson's correlation coefficient >0.50; P<0.05; Figs. 4 and 5). Co-expressed protein-coding genes with the whole human genome as the background were obtained via GO and KEGG pathway function enrichment analysis. A total of 623 protein-coding genes were enriched in 57 GO/KEGG terms according to GO functional annotation (data not shown). These important GO/KEGG terms were used in the ClueGo plugin in Cytoscape to organize an interaction network with similar functionality. Several groups of functionally associated GO terms were analysed in the present study and should be further investigated in the future. Based on the present analysis, the two circRNAs likely to influence protein-coding gene interactions and major biological processes, including vesicle transport and the cellular response to unfolded proteins, may be involved in tumorigenesis (Fig. 6).