Identification of a Robust Five-Gene Risk Model in Prostate Cancer: A Robust Likelihood-Based Survival Analysis

Abstract

Aim. In this paper, we aimed to develop and validate a risk prediction method using independent prognosis genes selected robustly in prostate cancer. Method. We considered 723 samples obtained from TCGA (the Cancer Genome Atlas), GSE46602, and GSE21032. Prostate cancer prognosis-related genes with $P<0.05$ were selected using Univariable Cox regression analysis. We then built the lowest AIC (Akaike information criterion score) optimal gene model using the “Rbsurv” package in TCGA train set. The coefficients were obtained by Multivariable Cox regression analysis. We named the new prognosis method CMU5. The CMU5 risk score was verified in TCGA test set, GSE46602, and GSE21032. Results. FAM72D, ARHGAP33, TACR2, PLEK2, and FA2H were identified as independent prognosis factors in prostate cancer patients. We built the computing model as follows: CMU5 risk score = 1.158$\ast$FAM72D + 1.737$\ast$ARHGAP33 − 0.737$\ast$TACR2 − 0.651$\ast$PLEK2 − 0.793$\ast$FA2H. The AUC of DFS was 0.809 in the train set (274 samples), 0.710 in the test set (273 samples), and 0.768 in the complete set (547 samples). The benign prediction capacity of CMU5 was verified by GSE46602 (36 samples; $\text{AUC}=0.6039$) and GSE21032 GPL5188 (140 samples; $\text{AUC}=0.7083$). Using the cut-off point of 2.056, a significant difference was shown between high- and low-risk groups. Conclusion. A prognosis-related risk score formula named CMU5 was built and verified, providing reliable prediction of prostate cancer outcome. This signature might provide a basis for individualized treatment of prostate cancer.

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.

1. Introduction

Prostate cancer (PCa) is the second most common male malignant tumor [1]. The mortality of PCa patients was reported as 40% over ten years, and the overall biological recurrence rate was 30.4% [2, 3]. PCa has strong heterogeneity. Its incidence is affected by factors such as age, ethnicity, and genetics. Tumor biological characteristics and prognosis vary greatly among individuals. Some slow-growing, weakly aggressive, low- and medium-risk tumors do not affect life expectancy [4]. Active local treatment of such patients may increase the occurrence of complications and affect quality of life; instead, active monitoring and other treatment methods can be adopted [5]. By contrast, other prostate cancer patients display high degrees of invasiveness and rapid progress. Therefore, it is important to stratify PCa patients with reasonable risk according to clinical and pathological parameters and to make clinical decisions based on life expectancy, health status, and subjective desires, then to formulate individualized treatment and follow-up plans. After radical prostatectomy, prostate cancer patients were treated with antiandrogenic drugs, and the PSA level was monitored trimonthly. Although the pathological stage and Gleason score were lower during surgery in patients, PSA increased quickly after surgery. Therefore, we urgently need an independent prognostic prediction method to assist us in grouping high- and low-risk patients in different stages and guide medication such as antiandrogens.

Advances in high-throughput sequencing and open source databases of tumors such as TCGA (the Cancer Genome Atlas) have enabled us to investigate the relationship between genes and prognosis. For example, HOXB5, GPC2, PGA5, and AMBN were used to establish an overall survival scoring model with $\text{AUC}=0.904$ [6]. SMIM22, NINL, NRG2, TOP2A, REPS2, and TPCN2 were shown to be biological recurrence prediction factors [7]. A methylation score formula consisting of HSPB1, CCND2, TIG1, DPYS, PITX2, and MAL was formulated to predict overall survival with $\text{AUC}=0.710$ [8]. Another study demonstrated that 60 miRNAs, 1578 mRNAs, and 61 lncRNAs were differentially expressed by a coexpression network [9]. These findings predicted prognosis using various methods and models. In this paper, we aimed to conduct the most reliable disease-free survival prediction model of prostate cancer. In addition, the degree of freedom of the prediction model should be limited, which reduced prediction costs. The robust method selects genes using the partial likelihood of the Cox model and identifies the optimal model using the lowest AICs (Akaike information criterion scores). Therefore, the robust method has more clinical significance, compared with the multivariable Cox proportional hazard method.

In this article, we have constructed a prediction method for the prognostic risk of prostate cancer patients, which is more accurate than TPSA. The prediction method was applied in different pathological stages and Gleason score subgroups and can effectively distinguish patients with different prognostic risk, providing a new method for actively monitoring prostate cancer patients.

2. Materials and Methods

2.1. Sample Source

The gene expression matrix and clinical follow-up information of PCa were obtained from TCGA (the Cancer Genome Atlas) database (https://www.cancer.gov/) [10]. In total, 547 samples were applied to this study, including 52 nontumor tissue samples and 495 tumor tissue samples. The days of new tumor events were considered DFS (disease-free survival) data. Patients were randomly assigned to a train set ($n=274$) containing 22 nontumor tissue samples and 252 tumor tissue samples and a test set ($n=273$). The final gene expression data were transformed by log2(exp+1); 36 samples with biochemical recurrence time were obtained from GSE46602 on the GPL570 Affymetrix Human Genome U133 Plus 2.0 Array [11]. A total of 140 samples with disease-free survival time and clinical stage were obtained from GSE21032 on the GPL5188 [12] platform in cBioPortal (http://www.cbioportal.org/) [13].

2.2. Prognosis-Related Gene Selection

The variation genes in each sample was identified as follows: the median and variance of the expression levels of a gene were greater than 20% of all genes. Subsequently, the relevance level between gene expression and disease-free survival status was evaluated in the train set. Univariable Cox regression analysis between gene expression and the disease-free survival state was performed using R with the “Survival” package [14]. The prognosis-related genes were determined with $P<0.05$ using “Survdiff” commands in R, and the prognosis correlations were analyzed using Kaplan–Meier survival curves and the ROC curve in GraphPad Prism 8.0 [15]. The Human Protein Atlas (HPA) (http://www.proteinatlas.org/) is an open source database. Expression of independent prognosis factors was evaluated on transcriptional and translational levels [16].

2.3. Pathway and Function Enrichment Analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID, v6.8) is a function enrichment tool that supplies biological explanations of gene lists and proteomic studies obtained from high-throughput sequencing [17]. Enrichment analysis for Gene Ontology [18] and KEGG pathway [19] was performed using DAVID, v6.8. Histogram was performed using the “ggplot2” package in R to show results [20]. $P<0.05$ indicated significance. GSEA (http://software.broadinstitute.org/gsea/index.jsp) was applied to show different pathways enriched in high- and low-risk groups.

2.4. Robust Selection of Prognostic-Related Genes

To establish the most reliable prognostic assessment model with the lowest degree of freedom, the robust principle and AICs were used to identify the best prognostic-related genes. The Rbsurv package in R was used to conduct robust likelihood-based survival analysis among survival-associated genes with the parameter as follows: $\text{iteration}\text{times}=100$ and $\max \text{concern}\text{genes}=20$ [21, 22]. Cluster analysis and expression levels of the best survival-associated genes were performed in the train set by running the “pheatmap” package [23]. Kaplan–Meier survival curves were used to evaluate survival differences between the two clusters [24]. Subsequently, prognosis factors were evaluated in TCGA using box plots between disease-free and adverse events. ROC curves and Kaplan–Meier survival analyses were applied to illustrate the independent prognosis values of the prognostic factors.

2.5. Risk Scoring System Establishment and Validation

CMU5 was established using the best prognosis-related genes. The estimated regression coefficients of each gene were calculated using multivariate Cox proportional hazard regression with the method “enter.” $\text{Log}\left[h\left(ti\right)/h0\left(ti\right)\right]=\beta 1\text{X}1+\beta 2\text{X}2+\beta 3\text{X}3+\cdots \beta \text{KXK}$, where $h\left(ti\right)$ is called the hazard function, and $h0\left(ti\right)$ is the baseline hazard. Terms $\text{X}1,\text{X}2,\text{X}3,\cdots \text{X}k$ are covariates and $\beta 1,\beta 2,\cdots \beta k$ are the corresponding regression coefficients. Based on the score formula analysis, the ROC curve was performed to evaluate the formula and calculate the best cut-off score with the maximal sensitivity and specificity in the train set [25]. Using same best cut-off point, ROC and Kaplan–Meier curves were generated in the test and complete sets to validate the risk score model and the cut-off point. The various risk groups were compared using the log-rank test [26]. Subsequently, the train set in TCGA, GSE46602 in the GEO database, GSE21032, and GPL5188 were used to verify the predictive ability of CMU5. Finally, the risk score formula CMU5 was fitted in various Gleason scores and stages in the complete set (547), and we validated the evaluative ability of the risk score formula using the log-rank test.

3. Results

3.1. Identification Prognosis-Related Genes

The overall process is presented in Figure 1. In total, 2707 protein coding genes were confirmed as prognosis-related genes using Cox proportional hazard regression with $P<0.05$ in TCGA train set (Table 1). The result of the KEGG enrichment analysis is shown in Figure 2(a), ordered by gene ratio. The 2707 genes were significantly enriched in the MAPK signaling pathway, neuroactive ligand-receptor interaction, focal adhesion, calcium signaling pathway, and cell cycle.

[figure omitted; refer to PDF]

Table 1

Univariable Cox regression of top 20 genes related to DFS survival.

Gene symbol	HR	HR.95L	HR.95H	$P$ value
PLEK2	0.490	0.391	0.614	5.80$E$-10
SPRED3	2.586	1.852	3.611	2.45$E$-08
TACR2	0.567	0.459	0.699	1.14$E$-07
RSPH10B	2.301	1.687	3.138	1.42$E$-07
TRIM73	2.033	1.549	2.670	3.22$E$-07
AMZ1	1.891	1.481	2.414	3.23$E$-07
FA2H	0.588	0.479	0.721	3.24$E$-07
ARHGAP33	2.373	1.694	3.324	4.94$E$-07
CPNE9	1.472	1.263	1.716	7.23$E$-07
C20orf203	1.564	1.309	1.868	8.63$E$-07
DPP4	0.699	0.605	0.807	1.01$E$-06
ASIC4	1.641	1.340	2.010	1.69$E$-06
CCDC180	1.950	1.483	2.564	1.71$E$-06
SEC61A2	3.051	1.925	4.837	2.09$E$-06
AL157935.2	2.018	1.509	2.698	2.15$E$-06
MXD3	1.933	1.471	2.539	2.23$E$-06
KMT5C	2.749	1.808	4.180	2.24$E$-06
SOX8	1.804	1.412	2.305	2.35$E$-06
FAM72D	1.751	1.378	2.224	4.52$E$-06

HR: hazard rate; DFS: disease-free survival.

[figures omitted; refer to PDF]

3.2. Identification of Robust Prognosis-Related Genes

To generate an optimal model with survival associated genes that were selected robustly, we selected 20 genes with the largest values of negative log-likelihoods. We obtained 20 prognosis-related gene signatures based on these 20 genes. The first model was generated using gene A with the largest value of negative log-likelihoods; the second model was generated using A+B, with B being the gene with the largest value of negative log-likelihoods except for that of A. The third model was generated by A+B+C, and others. The AICs [27] were calculated for each signature. The signature with the lowest AICs was selected, and it was considered to be the most reliable and feasible model with the minimum degree of freedom. The result is shown in Table 2, where the genes in the optimal signature are marked as ($\ast$). Tachykinin receptor 2 (TACR2), a family with sequence similarity 72-member D (FAM72D), pleckstrin-2 (PLEK2), fatty acid 2-hydroxylase gene (FA2H), and rho GTPase activating protein 33 (ARHGAP33) were strictly selected. Based on the expression level of the five genes, cluster analysis was performed in the train set to show the expression levels of the five prognosis-related genes (Figure 2(b)). Two clusters were identified based on the expression level, and the survival-related analysis showed that patients in cluster 2 had significantly higher risk than those in cluster 1 (Figure 2(c)). Subsequently, the box plots of the five genes between different statuses were generated in TCGA (547 samples). We found that FAM72D and ARHGAP33 were overexpressed in the adverse event group, and TACR2, FA2H, and PLEK2 were expressed in low amounts in the adverse event group (Figures 3(a)–3(e)). The AUC of TACR2, FAM72D, FA2H, PLEK2, and ARHGAP33 were 0.6155, 0.6615, 0.6613, 0.6613, and 0.6767, respectively (Figures 3(f)–3(j)). Kaplan–Meier survival analysis of the five genes showed higher survival risk in the high-expression group, suggesting that these factors are independent prognosis prediction factors (Figures 3(k)–3(o)). The five genes showed clinical correlation with the clinical stage and Gleason score (Figures 3(p)–3(y)). Based on the analysis mentioned in Figure 3, FAM72D and ARHGAP33 were considered to be positive risk factors and TACR2, FA2H, and PLEK2 were considered to be protective factors. The HPA database indicated that the color intensity of FAM72D and ARHGAP33 was higher in cancer tissue than in normal tissues, while TACR2, FA2H, and PLEK2 were the opposite (http://www.proteinatlas.org/) (Figure 4).

Table 2

The best prognosis-related model results selected by “Rbsurv” package in R.

Order	Gene	nloglik	AICs	Selected
0	0	275.87	551.73
1	TACR2	261.70	525.41	$$\ast$$
2	FAM72D	255.02	514.03	$$\ast$$
3	PLEK2	251.31	508.62	$$\ast$$
4	FA2H	248.07	504.15	$$\ast$$
5	ARHGAP33	243.38	496.76	$$\ast$$
6	TRIM74	242.55	497.11
7	TRIM73	242.03	498.05
8	SCNN1D	242.01	500.02
9	KRTAP5-1	241.66	501.32
10	CCDC180	240.98	501.97
11	MXD3	240.97	503.94
12	GPC2	239.75	503.5
13	SSPO	239.74	505.48
14	CPLX1	239.59	507.19
15	AL157935.2	237.49	504.99
16	SOX8	236.62	505.24
17	FGF17	233.77	501.53
18	SPRED3	233.21	502.42
19	SEC61A2	232.82	503.63

AIC: Akaike information criterion score; nloglik: negative log-likelihood.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

3.3. Risk Score Formula Establishment

To determine the relationship between signature and prognosis status, we built the risk score formula CMU5 using Cox proportional hazard regression with the method “enter”: $\log \left[h\left(ti\right)/h0\left(ti\right)\right]\left(\text{CMU}5\right)=1.158\ast FAM72D+1.737\ast ARHGAP33-0.737\ast TACR2-0.651\ast PLEK2-0.793\ast FA2H$. Each patient risk score was calculated in both the train and test sets in TCGA. The risk score curves and heatmaps of the train and test sets are shown in Figure 5. An ROC curve was generated to evaluate the risk score in the train set. We found that the AUC of CMU5 was 0.809 (Figure 6(a)), and the optimal threshold score with the maximal sensitivity and specificity was 2.0559. Subsequently, the ROC curve of TPSA, Gleason score, age, and stage were calculated in the train set (Figure 6(a)). The AUC of CMU5 was higher than that of the TPSA, Gleason score, age, and stage, suggesting that the CMU5 method better predicted the prognosis of PCa patients. Based on the threshold of 2.0559, the patients in the train set (274) were divided into a high-risk group (168) and a low-risk group (106). The Kaplan–Meier curve and log-rank test showed significant differences between the two groups ($P=3.7979e^{-12}$, $\text{HR}=6.604$) (Figure 6(b)). The univariable Cox proportional hazard analyses of age, Gleason score, and stage were performed in the train set (Table 3). We found that stage was a prognostic factor related to DFS. No significant survival risk was found for either the Gleason score or age. Based on these findings, we evaluated the predictive ability of CMU5 among various Gleason scores and stages in the train set. CMU5 risk increased as the stage and Gleason score increased (Figures 6(c) and 6(d)). CMU5 was also found to have good clinical relevance and prognostic ability in various subgroups ($P<0.05$) (Figures 6(e)–6(j)).

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

Table 3

Univariable Cox regression of age Gleason score and stage.

Term	Count	HR (95% CI)	$P$ value
Age
<60	121	1
≥60	153	0.821 (0.491–1.373)	0.450
Gleason score
6	20	1
7	122	0.740 (0.187–2.938)	0.632
8	29	2.049 (0.647–6.486)	0.268
9–10	81	2.820 (1.243–6.397)	0.072
$N$	22
Stage
Normal	22	1
II	78	3.714 (0.565–24.40)	0.172
III	164	3.249 (1.764–6.538)	0.003
IV	5	12.830 (1.253–131.327)	0.032
NA	5
Risk
Low	168	1
High	106	6.604 (3.842–11.351)	<0.001

CI: confidence interval; HR: hazard ratio; DFS: disease-free-survival.

3.4. Risk Score Formula Evaluation

To evaluate the risk score formula and threshold score, the complete set was applied to evaluate the results. The AUC of the complete set (547) was 0.768 (Figure 7(a)). Using the same cut-off point as 2.0559, the complete set was divided into the low-risk group and the high-risk groups. We found a significant survival risk difference ($P<0.0001$, $\text{HR}=4.269$, Figure 7(b)). In addition, CMU5 showed good prognostic ability and close clinical relevance in various subgroups in TCGA complete set (Figures 7(c)–7(j)).

[figures omitted; refer to PDF]

3.5. Test Set and External Data Validation

The AUC of the test set (273 samples) was 0.710 (Figure 8(a)). Using the cut-off point of 2.0559, the test set was divided into low-risk and high-risk groups. There was a significant survival risk difference (Figure 8(b)) ($P=0.0007$, $\text{HR}=2.649$). Subsequently, the CMU5 was shown to have a positive correlation with the stage and Gleason score (Figures 9(c) and 9(d)). Based on this analysis of TCGA data, CMU5 showed good prognostic value. To further confirm the significance of CMU5, we adopted two external datasets. The ROC curve ($\text{AUC}=0.7083$) and Kaplan–Meier curve ($P=0.0004$, $\text{HR}=3.773$) of GSE21032 (140 samples) are shown in Figures 8(e) and 8(f). The ROC curve ($\text{AUC}=0.6039$) and Kaplan–Meier curve ($P=0.0073$, $\text{HR}=2.976$) of GSE46602 (36 samples) are shown in Figures 8(i) and 8(j). The clinically relevant phenotypic analysis also illustrated the positive correlation of CMU5 with the stage and Gleason score in GSE21032 and GSE46602 (Figures 8(g), 8(h), 8(k), and 8(l)). Based on the verification of these two external data sets, we had once again identified the prognostic value of the risk evaluation model CMU5 in PCa patients. Taken together, the data suggest that the CMU5 method showed good prognostic prediction ability in PCa patients.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

3.6. GSEA

To investigate the changes of the pathway in the low-risk and high-risk groups, GSEA analysis was used. The results are shown in Figure 8. The homologous-recombination pathway, the DNA-replication pathway, the mismatch-repair pathway, the cell-cycle pathway, and the base excision repair pathway were significantly related to the high-risk group, suggesting an active cell proliferation process occurring in the high-risk group (Figures 9(a)–9(e)). In the low-risk group, the arginine and proline metabolism pathway, the butanoate-metabolism pathway, the glycosaminoglycan-degradation pathway, the propanoate-metabolism pathway, and the valine and isoleucine-degradation pathway were significantly enriched, suggesting that low metabolic levels might contribute to better prognosis compared with the high-risk group (Figures 9(f)–9(j)).

4. Discussion

PCa is the most common malignant tumor of the male genitourinary system. According to the 2018 GLOBOCAN statistics of the World Health Organization, the incidence of PCa ranks the second among all male malignancies worldwide, second only to lung cancer [28]. PSA testing is recommended for patients with a life expectancy of more than 10 years, and further risk assessment should be conducted for asymptomatic patients with normal DRE and a $\text{PSA}\text{level}<10\text{ng}/\text{ml}$ [29]. In this paper, we constructed a method to supplement the prognostic risk of patients and a prediction method with higher accuracy than TPSA. Our scoring method was effective in differentiating patients with different prognostic risks through three validation correlations, providing a new method for active surveillance of PCa.

With the maturity of bioinformatic analysis in recent years, there have appeared many methods to predict the risk of PCa based on gene expression. We summarize the existing prediction models and improve their shortcomings. Xu et al. built a prediction model of overall survival including four mRNA ($\text{AUC}=0.904$) [6] and conducted a projection for overall survival analysis. Statistical indicators were significant; however, there was lack of external validation analysis of data sets. In addition, we found in the follow-up data of TCGA death events which occurred in only nine samples; therefore, in our article, we selected disease-free survival analysis. The CAPRA score constructed by Ahmad et al. [8] predicted the risk of early PCa with AUC 0.710 for 10 years; however, we obtained a more accurate prediction model ($\text{AUC}=0.809$) and verified it in two other cohorts. Therefore, to the best of our knowledge, we have obtained a prediction model with the lowest degree of freedom, the highest accuracy, and consistently good prediction in various cohorts and subgroups.

TACR2, FAM72D, PLEK2, FA2H, and ARHGAP33 were first proposed as independent predictors of PCa in this paper. We built the CMU5 score based on these five protein coding genes. We applied the days to new tumor events as the parameters of disease-free survival, which were related to tumor recurrence and another adverse events. The robust method was applied because it builds multiple gene models sequentially with survival-associated genes selected robustly. The risk score formula and the best cut-off point were both verified using the Kaplan–Meier curve and log-rank tests in the test and complete sets.

TACR2, PLEK2, and FA2H were considered protective factors in PCa. Tachykinin receptor 2 (TACR2), also called NK2R, is one of the family of genes that encodes receptors for tachykinins and interacts with G proteins and seven hydrophobic transmembrane regions. Tachykinins are modulators of the immune system, related to the generation, activation, development, and migration of immune cells [30]. Tachykinins also mediate T cell differentiation; Zhang et al. found that CD8⁺ T cells were significantly decreased after treatment with tachykinin antagonist CD8⁺ T cells which play crucial roles in cellular immunity, providing protection from tumor cell infiltration. Pleckstrin-2 (PLEK2) is associated with membrane-bound phosphatidylinositols generated by phosphatidylinositol 3-kinase. Bach et al. suggested that pleckstrin-2 binds to membrane-associated phosphatidylinositols regulated by PI3K, thereby promoting the actin cytoskeleton in lymphocyte spreading and immune synapse formation [31]. Fatty acid 2-hydroxylase (FA2H) was shown to play a crucial role in regulating hedgehog signaling and the suppression of gastric tumor growth. Downregulation of the hedgehog signaling pathway also suppressed PCa cell proliferation and invasion [32]. These findings suggest that TACR2, PLEK2, and FA2H provide protection from tumor invasion; they were applied as protective factors in our risk score method.

FAM72D and ARHGAP33 are risk factors for PCa. The family with sequence similarity 72 member D (FAM72D) is also known as GCUD2; it is a poor prognostic gene of myeloma and control cell proliferation and survival in the FOXM1 transcription factor network [33]. FAM72 paralogs are upregulated in tumor cells and are related to mitotic cell cycle genes that promote the formation of centrosomes and mitotic spindles and act as prognostic biomarkers for glioblastoma [34]. Rho GTPase activating protein 33 (ARHGAP33) is a high-affinity receptor for the brain-derived neurotrophic factor [35]. Chen et al. [36] suggested that ARHGAP9, 15, 18, 19, 25, and 30 were associated with breast cancer. To our best knowledge, ours is the first study to identify ARHGAP33 and PLEK2 as PCa prognosis factors. Nevertheless, the mechanisms underlying the effects of these genes on prognosis in PCa require further research.

In this study, we established a risk score called CMU5 that divides PCa patients into different groups, and we provided the disease-free survival prediction time in high-risk and low-risk groups. With the CMU5 score support, we can distinguish high-risk patients with low Gleason scores to provide patients with individualized treatment. The CMU5 score was verified to be reliable in two other external datasets. These data suggest that both CMU5 and the threshold value make sense in terms of disease-free survival time and status. Nevertheless, because of the limitations of our research methods, there was no in-depth study mechanism of action of the factors in the scoring model, and the scoring algorithm requires further verification based on basic science research.

5. Conclusions

We developed a five-gene signature for survival prediction in PCa patients from TCGA. A five-gene signature (TACR2, FAM72D, PLEK2, FA2H, and ARHGAP33) named CMU5, with genes selected robustly, was identified using the “Rbsurv” package. Based on the cut-off of 2.056, high-risk and low-risk groups were identified. Based on the verification of the benign nature and evaluation effect, the CMU5 score might have potential prognostic and therapeutic implications for PCa patients.

Acknowledgments

This work was supported by the China Medical University Youth Backbone Support Program (QGZD2018029). We want to thank TCGA, GEO, and cBioPortal for the free usage.

References

[1] B. Farhood, K. Mortezaee, H. Haghi-Aminjan, N. Khanlarkhani, E. Salehi, M. S. Nashtaei, M. Najafi, A. Sahebkar, "A systematic review of radiation-induced testicular toxicities following radiotherapy for prostate cancer," Journal of Cellular Physiology, vol. 234 no. 9, pp. 14828-14837, DOI: 10.1002/jcp.28283, 2019.

[2] D. Ma, Z. Zhou, B. Yang, Q. He, Q. Zhang, X. H. Zhang, "Association of molecular biomarkers expression with biochemical recurrence in prostate cancer through tissue microarray immunostaining," Oncology Letters, vol. 10 no. 4, pp. 2185-2191, DOI: 10.3892/ol.2015.3556, 2015.

[3] M. Van Hemelrijck, Y. Folkvaljon, J. Adolfsson, O. Akre, L. Holmberg, H. Garmo, P. Stattin, "Causes of death in men with localized prostate cancer: a nationwide, population-based study," BJU International, vol. 117 no. 3, pp. 507-514, DOI: 10.1111/bju.13059, 2016.

[4] F. H. Schröder, J. Hugosson, M. J. Roobol, T. L. Tammela, M. Zappa, V. Nelen, M. Kwiatkowski, M. Lujan, L. Määttänen, H. Lilja, L. J. Denis, F. Recker, A. Paez, C. H. Bangma, S. Carlsson, D. Puliti, A. Villers, X. Rebillard, M. Hakama, U. H. Stenman, P. Kujala, K. Taari, G. Aus, A. Huber, T. van der Kwast, R. van Schaik, H. J. de Koning, S. M. Moss, A. Auvinen, ERSPC Investigators, "Screening and prostate cancer mortality: results of the European Randomised Study of Screening for Prostate Cancer (ERSPC) at 13 years of follow-up," The Lancet, vol. 384 no. 9959, pp. 2027-2035, DOI: 10.1016/S0140-6736(14)60525-0, 2014.

[5] A. Bill-Axelson, L. Holmberg, M. Ruutu, H. Garmo, J. R. Stark, C. Busch, S. Nordling, M. Häggman, S. O. Andersson, S. Bratell, A. Spångberg, J. Palmgren, G. Steineck, H. O. Adami, J. E. Johansson, SPCG-4 Investigators, "Radical prostatectomy versus watchful waiting in early prostate cancer," The New England Journal of Medicine, vol. 364 no. 18, pp. 1708-1717, DOI: 10.1056/NEJMoa1011967, 2011.

[6] N. Xu, Y. P. Wu, H. B. Yin, X. Y. Xue, X. Gou, "Molecular network-based identification of competing endogenous RNAs and mRNA signatures that predict survival in prostate cancer," Journal of Translational Medicine, vol. 16 no. 1,DOI: 10.1186/s12967-018-1637-x, 2018.

[7] F. Li, J. P. Ji, Y. Xu, R. L. Liu, "Identification a novel set of 6 differential expressed genes in prostate cancer that can potentially predict biochemical recurrence after curative surgery," Clinical & Translational Oncology, vol. 21 no. 8, pp. 1067-1075, DOI: 10.1007/s12094-018-02029-z, 2019.

[8] A. S. Ahmad, N. Vasiljević, P. Carter, D. M. Berney, H. Møller, C. S. Foster, J. Cuzick, A. T. Lorincz, "A novel DNA methylation score accurately predicts death from prostate cancer in men with low to intermediate clinical risk factors," Oncotarget, vol. 7 no. 44, pp. 71833-71840, DOI: 10.18632/oncotarget.12377, 2016.

[9] Y. Ye, S. L. Li, S. Y. Wang, "Construction and analysis of mRNA, miRNA, lncRNA, and TF regulatory networks reveal the key genes associated with prostate cancer," PLoS One, vol. 13 no. 8, article e0198055,DOI: 10.1371/journal.pone.0198055, 2018.

[10] K. Tomczak, P. Czerwińska, M. Wiznerowicz, "Review The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge," Współczesna Onkologia, vol. 1A, pp. 68-77, DOI: 10.5114/wo.2014.47136, 2015.

[11] M. M. Mortensen, S. Høyer, A. S. Lynnerup, T. F. Ørntoft, K. D. Sørensen, M. Borre, L. Dyrskjøt, "Expression profiling of prostate cancer tissue delineates genes associated with recurrence after prostatectomy," Scientific Reports, vol. 5 no. 1, article 16018,DOI: 10.1038/srep16018, 2015.

[12] B. S. Taylor, N. Schultz, H. Hieronymus, A. Gopalan, Y. Xiao, B. S. Carver, V. K. Arora, P. Kaushik, E. Cerami, B. Reva, Y. Antipin, N. Mitsiades, T. Landers, I. Dolgalev, J. E. Major, M. Wilson, N. D. Socci, A. E. Lash, A. Heguy, J. A. Eastham, H. I. Scher, V. E. Reuter, P. T. Scardino, C. Sander, C. L. Sawyers, W. L. Gerald, "Integrative genomic profiling of human prostate cancer," Cancer Cell, vol. 18 no. 1, pp. 11-22, DOI: 10.1016/j.ccr.2010.05.026, 2010.

[13] J. Gao, B. A. Aksoy, U. Dogrusoz, G. Dresdner, B. Gross, S. O. Sumer, Y. Sun, A. Jacobsen, R. Sinha, E. Larsson, E. Cerami, C. Sander, N. Schultz, "Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal," Science Signaling, vol. 6 no. 269, article pl1,DOI: 10.1126/scisignal.2004088, 2013.

[14] J. O'Quigley, T. Moreau, "Cox's regression model: computing a goodness of fit statistic," Computer Methods and Programs in Biomedicine, vol. 22 no. 3, pp. 253-256, DOI: 10.1016/0169-2607(86)90001-5, 1986.

[15] J. Kishore, M. K. Goel, P. Khanna, "Understanding survival analysis: Kaplan-Meier estimate," International Journal of Ayurveda Research, vol. 1 no. 4, pp. 274-278, DOI: 10.4103/0974-7788.76794, 2010.

[16] M. Uhlen, L. Fagerberg, B. M. Hallstrom, C. Lindskog, P. Oksvold, A. Mardinoglu, A. Sivertsson, C. Kampf, E. Sjostedt, A. Asplund, I. Olsson, K. Edlund, E. Lundberg, S. Navani, C. A. K. Szigyarto, J. Odeberg, D. Djureinovic, J. O. Takanen, S. Hober, T. Alm, P. H. Edqvist, H. Berling, H. Tegel, J. Mulder, J. Rockberg, P. Nilsson, J. M. Schwenk, M. Hamsten, K. von Feilitzen, M. Forsberg, L. Persson, F. Johansson, M. Zwahlen, G. von Heijne, J. Nielsen, F. Ponten, "Tissue-based map of the human proteome," Science, vol. 347 no. 6220,DOI: 10.1126/science.1260419, 2015.

[17] D. W. Huang, B. T. Sherman, Q. Tan, J. R. Collins, W. G. Alvord, J. Roayaei, R. Stephens, M. W. Baseler, H. C. Lane, R. A. Lempicki, "The DAVID gene functional classification tool: a novel biological module-centric algorithm to functionally analyze large gene lists," Genome Biology, vol. 8 no. 9,DOI: 10.1186/gb-2007-8-9-r183, 2007.

[18] H. Ogata, S. Goto, W. Fujibuchi, M. Kanehisa, "Computation with the KEGG pathway database," Bio Systems, vol. 47 no. 1-2, pp. 119-128, DOI: 10.1016/S0303-2647(98)00017-3, 1998.

[19] M. Ashburner, C. A. Ball, J. A. Blake, D. Botstein, H. Butler, J. M. Cherry, A. P. Davis, K. Dolinski, S. S. Dwight, J. T. Eppig, M. A. Harris, D. P. Hill, L. Issel-Tarver, A. Kasarskis, S. Lewis, J. C. Matese, J. E. Richardson, M. Ringwald, G. M. Rubin, G. Sherlock, "Gene Ontology: tool for the unification of biology," Nature Genetics, vol. 25 no. 1, pp. 25-29, DOI: 10.1038/75556, 2000.

[20] ggplot2, Elegant Graphics for Data Analysis, 2016.

[21] H. Cho, A. Yu, S. Kim, J. Kang, S. M. Hong, "Robust likelihood-based survival modeling with microarray data," Journal of Statistical Software, vol. 29 no. 1, 2009.

[22] G. Renaud, U. Stenzel, T. Maricic, V. Wiebe, J. Kelso, "deML: Robust demultiplexing of Illumina sequences using a likelihood-based approach," Bioinformatics, vol. 31 no. 5, pp. 770-772, DOI: 10.1093/bioinformatics/btu719, 2015.

[23] R. Kolde, "pheatmap: Pretty Heatmaps. R package version 1.0.12," . 2019, https://CRAN.R-project.org/package=pheatmap

[24] R. Porcher, "CORR insights(®): Kaplan-Meier survival analysis overestimates the risk of revision arthroplasty: a meta-analysis," Clinical Orthopaedics and Related Research, vol. 473 no. 11, pp. 3443-3445, DOI: 10.1007/s11999-015-4291-0, 2015.

[25] T. Fawcett, "An introduction to ROC analysis," Pattern Recognition Letters, vol. 27 no. 8, pp. 861-874, DOI: 10.1016/j.patrec.2005.10.010, 2006.

[26] J. Sun, Q. Zhao, X. Zhao, "Generalized log-rank tests for interval-censored failure time data," Scandinavian Journal of Statistics, vol. 32 no. 1, pp. 49-57, DOI: 10.1111/j.1467-9469.2005.00414.x, 2005.

[27] E. J. Wagenmakers, S. Farrell, "AIC model selection using Akaike weights," Psychonomic Bulletin & Review, vol. 11 no. 1, pp. 192-196, DOI: 10.3758/BF03206482, 2004.

[28] F. Bray, J. Ferlay, I. Soerjomataram, R. L. Siegel, L. A. Torre, A. Jemal, "Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries," CA: a Cancer Journal for Clinicians, vol. 68 no. 6, pp. 394-424, DOI: 10.3322/caac.21492, 2018.

[29] T. B. L. Lam, S. MacLennan, P. P. M. Willemse, M. D. Mason, K. Plass, R. Shepherd, R. Baanders, C. H. Bangma, A. Bjartell, A. Bossi, E. Briers, A. Briganti, K. T. Buddingh, J. W. F. Catto, M. Colecchia, B. W. Cox, M. G. Cumberbatch, J. Davies, N. F. Davis, M. de Santis, P. Dell’Oglio, A. Deschamps, J. F. Donaldson, S. Egawa, C. D. Fankhauser, S. Fanti, N. Fossati, G. Gandaglia, S. Gillessen, N. Grivas, T. Gross, J. P. Grummet, A. M. Henry, A. Ingels, J. Irani, M. Lardas, M. Liew, D. W. Lin, L. Moris, M. I. Omar, K. H. Pang, C. C. Paterson, R. Renard-Penna, M. J. Ribal, M. J. Roobol, M. Rouprêt, O. Rouvière, G. Sancho Pardo, J. Richenberg, I. G. Schoots, J. P. M. Sedelaar, P. Stricker, D. Tilki, S. Vahr Lauridsen, R. C. N. van den Bergh, T. van den Broeck, T. H. van der Kwast, H. G. van der Poel, G. J. L. H. van Leenders, M. Varma, P. D. Violette, C. J. D. Wallis, T. Wiegel, K. Wilkinson, F. Zattoni, J. M. O. N’Dow, H. van Poppel, P. Cornford, N. Mottet, "EAU-EANM-ESTRO-ESUR-SIOG prostate cancer guideline panel consensus statements for deferred treatment with curative intent for localised prostate cancer from an international collaborative study (DETECTIVE study)," European Urology, vol. 76 no. 6, pp. 790-813, DOI: 10.1016/j.eururo.2019.09.020, 2019.

[30] S. Khorasani, N. Boroumand, F. Lavi Arab, S. I. Hashemy, "The immunomodulatory effects of tachykinins and their receptors," Journal of Cellular Biochemistry, vol. 121 no. 5-6, pp. 3031-3041, DOI: 10.1002/jcb.29668, 2020.

[31] T. L. Bach, W. T. Kerr, Y. Wang, E. M. Bauman, P. Kine, E. L. Whiteman, R. S. Morgan, E. K. Williamson, E. M. Ostap, J. K. Burkhardt, G. A. Koretzky, M. J. Birnbaum, C. S. Abrams, "PI3K regulates pleckstrin-2 in T-cell cytoskeletal reorganization," Blood, vol. 109 no. 3, pp. 1147-1155, DOI: 10.1182/blood-2006-02-001339, 2006.

[32] P. Xu, F. Cai, X. Liu, L. Guo, "LKB1 suppresses proliferation and invasion of prostate cancer through hedgehog signaling pathway," International Journal of Clinical and Experimental Pathology, vol. 7 no. 12, pp. 8480-8488, 2014.

[33] F. Chatonnet, A. Pignarre, A. A. Sérandour, G. Caron, S. Avner, N. Robert, A. Kassambara, A. Laurent, M. Bizot, X. Agirre, F. Prosper, J. I. Martin-Subero, J. Moreaux, T. Fest, G. Salbert, "The hydroxymethylome of multiple myeloma identifies FAM72D as a 1q21 marker linked to proliferation," Haematologica, vol. 105 no. 3, pp. 774-783, DOI: 10.3324/haematol.2019.222133, 2020.

[34] C. S. Rahane, A. Kutzner, K. Heese, "A cancer tissue-specific FAM72 expression profile defines a novel glioblastoma multiform (GBM) gene-mutation signature," Journal of Neuro-Oncology, vol. 141 no. 1, pp. 57-70, DOI: 10.1007/s11060-018-03029-3, 2019.

[35] T. Nakazawa, R. Hashimoto, K. Sakoori, Y. Sugaya, A. Tanimura, Y. Hashimotodani, K. Ohi, H. Yamamori, Y. Yasuda, S. Umeda-Yano, Y. Kiyama, K. Konno, T. Inoue, K. Yokoyama, T. Inoue, S. Numata, T. Ohnuma, N. Iwata, N. Ozaki, H. Hashimoto, M. Watanabe, T. Manabe, T. Yamamoto, M. Takeda, M. Kano, "Emerging roles of ARHGAP33 in intracellular trafficking of TrkB and pathophysiology of neuropsychiatric disorders," Nature Communications, vol. 7 no. 1, article 10594,DOI: 10.1038/ncomms10594, 2016.

[36] W.-X. Chen, M. Lou, L. Cheng, Q. Qian, L.-Y. Xu, L. Sun, Y.-L. Zhu, H. Dai, "Bioinformatics analysis of potential therapeutic targets among ARHGAP genes in breast cancer," Oncology Letters, vol. 18 no. 6, pp. 6017-6025, DOI: 10.3892/ol.2019.10949, 2019.