Gene expression and miRNA data and corresponding clinical information were obtained from the TCGA- prostate adenocarcinoma (PRAD) dataset using the UCSC Xena browser (https://xenabrowser.net/). TCGA data from 333 PCa patients who underwent RP were included in the present study. TCGA samples with a high Gleason score (≥8) (n = 129) or a low Gleason score (≤6) patients (n = 27) were used as the discovery cohort to construct a gene classifier (associated with miRNA-mediated regulation) for predicting prognosis. The entire set of TCGA samples (n = 333) was used as an internal validation cohort. Another cohort of 106 samples from 100 PCa patients after RP was obtained from the GEO database (accession number: GSE54460 (Qi et al., 2014)). Duplicates from six patients were removed, and the remaining 100 samples were used as the external validation cohort. The clinical characteristics of the PCa patients after RP in the discovery cohort, internal validation cohort and external validation cohort are summarized in Table 1 .

For samples in the TCGA dataset, gene expression and miRNA profiles were converted to log₂(x + 1) of the initial expression value. Only protein-coding gene profiles were extracted for further analysis by identifying genes with protein products in the HUGO Gene Nomenclature Committee (HGNC) database (https://www.genenames.org/). Given that failures of gene and miRNA expression detection might exist, we only selected genes and miRNAs that were abundantly expressed for further analysis. The criteria were as follows: (1) expression level >0; (2) appeared in >50% of all specimens. For samples in the external validation cohort (from the GEO database), the gene expression data were preprocessed into fragment per kilobase per million reads (FPKM) values.

We constructed a gene classifier associated with miRNA-mediated regulation using the following four steps. Step 1 involved obtaining candidate genes and miRNAs. Logistic regression was applied to the discovery cohort to identify candidate genes and miRNAs (with p < 0.05) associated with the Gleason score. Step 2 involved obtaining significantly negitively correlated miRNA-gene pairs. A Pearson correlation analysis was applied to each miRNA-gene pair based on the candidate genes and miRNAs identified in step 1. A miRNA-gene pair with a Pearson correlation coefficient < -0.4 and p < 0.05 was defined as having a statistically significant negative correlation and used for further study. Step 3 involved, obtaining the target genes associated with miRNA-mediated regulation. We retrieved the potential target genes of the significantly negitively correlated miRNAs from miRWalk 3 (http://mirwalk.umm.uni-heidelberg.de/). Target genes that overlapped, i.e., were also significantly negitively correlated genes, were employed to develop the classifier. Step 4 involved, least absolute shrinkage and selector operation (LASSO) regression and a decision tree to further narrow down the variables and create two candidate classifiers. The Glmnet package was used to perform logistic regression and LASSO regression while the rpart package was used to construct the decision tree, in R software version 3.4.0.

Receiver operating characteristic (ROC) curve analysis was applied to assess the classifier’s ability to distinguish samples with a high or low Gleason score. Time-dependent ROC curve analyses were applied to assess the classifier’s ability to predict the 3-year clinical recurrence-free survival (RFS) rate and BCR-free survival (BCRFS) rate, using the survival ROC package in R software. Area under the curve (AUC) values were calculated to estimate the prediction ability of the classifier or related clinical factors. The Youden index from the time-dependent ROC curve analysis was employed to define the cutoff to split samples into high- and low-risk groups. The RFS and BCRFS were compared between the two groups according to the Kaplan-Meier method with the log-rank test.

Based on a Pearson correlation analysis of the genes in the classifier, genes with an absolute Pearson correlation coefficient >0.4 were identified as co-expressed genes (Seon-Kyu et al., 2015). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of the co-expressed genes were conducted using ClusterProfiler with p < 0.01 as the cut-off criterion (Yu et al., 2012). A network that involved the above-mentioned overlapping miRNA-gene pairs was visualized using Cytoscape 3.5.1.

After preprocessing, each abundantly expressed gene and miRNA from samples in the discovery cohort was subjected to univariate logistic regression using binomial Gleason score as the dependent variable. Consequently, 3,732 genes and 98 miRNAs ( Additional File 1 ) were found to be significant. Next, 1,235 significantly negitively correlated miRNA-gene pairs were identified using a Pearson correlation analysis ( Additional File 2 ). Theredfter, we retrieved the potential target genes of the significantly negitively correlated miRNAs from miRWalk 3, obtaining 158,948 miRNA-gene pairs ( Additional File 3 ). The miRNA-gene pairs involving genes that overlapped, i.e., were present in both sets ( Additional File 4 ) were visualized in a network, with 79 genes and 28 miRNAs involved, as shown in Figure 1 .

These 79 genes were first subjected to LASSO regression, still using the binomial Gleason score as the dependent variable. LASSO regression is a parsimonious model that involves L1 regularization. λ is the coefficient of the penalty term in L1 regularization and as λ increases, the regression coefficients approach zero. Variables with non-zero regression coefficients are the variables most strongly associated with the dependent variable. Tenfold cross-validation was performed to determine the best λ value, with the AUC as the criterion. A series of models were constructed for variable selection ( Figure 2A ), among which the model with the highest AUC value was selected as the best model ( Figure 2B ). Three genes (FAM72B, GNE, and TRIM46) had non-zero coefficients in the best model. These three genes and their coefficients formed the following prognostic index (PI), representing a candidate classifier:

The PI of samples were compared between groups with high and low Gleason scores using a violin plot, as shown in Figure 2C . To further evaluate the sensitivity and specificity of Gleason score prediction, we conducted a ROC curve analysis of PI. As shown in Figure 2D , PI exhibited a remarkable ability regarding Gleason score prediction (AUC = 0.927).

Furthermore, we used the 79 genes and the binomial Gleason score to construct a decision tree using the classification and regression tree (CART) algorithm. The decision tree with minimal cross-validation error was generated, as shown in Figure 2E , comprising two genes, FAM72B and GNE. As shown in the decision tree, FAM72B had almost complete dominance, while GNE came second. The coefficients in the PI indicated the same pattern. Therefore, the candidate classifiers obtained from the two machine-learning algorithms agreed well with each other, with the exception of TRIM46 in the LASSO-based classifier. To estimate the performance of the two classifiers, their sensitivity and specificity in classifying Gleason score were calculated, as shown in Figure 2F . The accuracy of the LASSO-based classifier was a little higher than that of the decision tree-based classifier. Thus, the PI derived from the LASSO-based classifier was identified as the final classifier.

After mapping the three genes in the classifier to the corresponding overlapping miRNA-gene pairs, we found that FAM74B was regulated by both hsa-miR-133a-3p and hsa-miR-222-3p, while GNE and TRIM46 were regulated by hsa-miR-1301-3p and hsa-miR-30c-2-3p, respectively. A heatmap of the expression of the three genes and four miRNAs was generated to gain more understanding of the correlations of individual genes and miRNAs with the Gleason score ( Figure 2G ). The heatmap shows that FAM72B, TRIM46 and hsa-miR-1301-3p were positively correlated with the Gleason score, whereas GNE, hsa-miR-133a-3p, hsa-miR-30c-2-3p, and hsa-miR-222-3p were negatively correlated with Gleason score.

As described in the Materials and Methods, the three genes in the classifier were significantly negatively correlated with their corresponding miRNAs. Scatter plots of the four miRNA-gene pairs’ expression and corresponding linear regression lines were generated, as shown in Figure 3 . The Pearson correlation coefficients were -0.475, -0.411, -0.405, and -0.551 for the miRNA-gene pairs hsa-miR-133a-3p-FAM74B, hsa-miR-222-3p-FAM74B, hsa-miR-1301-3p-GNE, and hsa-miR-30c-2-3p-TRIM46, respectively.

First, we examined the associations between each gene included in the classifier and the clinical recurrence of PCa patients after RP in the internal validation cohort. RFS differences between the low- and high-risk groups were analyzed using the Kaplan-Meier method with the log-rank test, as shown in Figures 4A–C . We found that the three genes were all significantly associated with the RFS of PCa patients after RP (p < 0.001 for FAM74B, p < 0.001 for GNE and p = 0.0134 for TRIM46). Identical analyses for the four upstream miRNAs were also conducted, as shown in Figures 4D–G . Except for hsa-miR-1301-3p (p = 0.0829), all other miRNAs were significantly associated with RFS ((p-0.0121 for hsa-miR-133a-3p, p < 0.001 for hsa-miR-222-3p and p = 0.025 for hsa-miR-30c-2-3p).

Second, we examined the classifier’s ability to predict clinical recurrence in patients after RP. With the use of the best cutoff point, the patients in this cohort were further divided into a high-risk group (n = 119) and a low-risk group (n = 173). The RFS difference was analyzed between the low- and high-risk groups ( Figure 5A ) in the internal validation cohort. The survival curves indicated a highly significant difference between the low- and high-risk groups (p < 0.001, log rank statistic = 20.7). To gain a deeper insight into the classifier, it was subjected to univariate Cox proportional hazards regression, as were four candidate clinical prognostic factors (age, Gleason score, PSA, and pathological stage) in the internal validation cohort. The results indicated that high classifier value, high Gleason score, high PSA value and Regional pathological stage were associated with significantly shorter RFS. The fmultivariate analysis confirmed that the classifier (hazard ratio [HR]: 1.708, 95%CI: 1.180–2.472, p < 0.001) and Gleason score (HR: 1.720, 95% CI: 1.089–2.716, p = 0.02) were independent risk factors for RFS. To further evaluate the sensitivity and specificity of the RFS survival prediction, we conducted time-dependent ROC curve analyses for the classifier and clinical factors. As shown in Figure 5B , the classifier exhibited a remarkable ability to predict RFS, as the AUC value was 0.733, which was higher than the AUC value for the Gleason score (0.718) ( Figure 5B ). Additionally, a combination of the classifier and three prognostic factors (Gleason score, PSA and pathological stage) achieved the best ability to predict RFS (AUC = 0.752). The results of the univariate and multivariate analyses of RFS in PCa patients after RP in the internal validation cohort are summarized in Table 2 . Moreover, the classifier performed stably regarding risk stratification in both the discovery cohort (p = 0.0133, log rank statistic = 6.13) and the subgroup of the internal validation cohort with a Gleason score of = 7 (p = 0.02, log rank statistic = 5.41) ( Figures 5C, D ).

Third, we examined the classifier’s prognostic ability to predict the BCR of patients after RP. The classifier’s output was calculated using the weight coefficients used in the above-mentioned prognostic index (PI). With the use of the best cut-off point, the patients in the external validation cohort were divided into a high-risk group (n = 48) and a low risk group (n = 48). The BCRFS difference was compared between the low- and high-risk groups ( Figure 5E ). The survival curves still indicated a significant difference (p = 0.0338, log rank statistic = 4.51). A similar result was obtained in the subgroup of the external validation cohort with Gleason score of 7 (p = 0.0227, log rank statistic = 5.12) ( Figure 5F ).

To explore the underlying moleculer mechanisms of the three genes included in the classifier, we performed a functional enrichment analysis of the three genes’ co-expressed genes in the internal vlidation cohort. Genes with an absolute Pearson correlation coefficient >0.4 were identified as co-expressed genes. There were 673, 165, and 153 genes co-expressed with FAM72B, GNE, and TRIM46, respectively ( Additional file 5 ). Finally, we obtained 894 co-expressed genes after removal of the duplicate genes. GO and KEGG enrichment analyses were conducted on these co-expressed genes, with Biological Process (BP) GO terms and pathways with p < 0.01 being defined as significantly enriched. The Top 15 of all 245 significantly enriched BP GO terms are visualized in Figure 6A , including chromosome segregation, mitotic nuclear division, nuclear chromosome segregation, sister chromatid segregation, DNA replication, and mitotic sister chromatid segregation. Significantly enriched KEGG pathways (n = 14) are shown in Figure 6B , including Cell cycle, DNA replication, Homologous recombination, Mismatch repair, Fanconi anemia pathway, Base excision repair, Oocyte meiosis, Nucleotide excision repair, Progesterone-mediated oocyte maturation and beta-Alanine metabolism, Valine, leucine and isoleucine degradation, Cellular senescence, Human T-cell leukemia virus 1 infection, and p53 signaling pathway.