Construction autophagy-related prognostic risk signature combined with clinicopathological validation analysis for survival prediction of kidney renal papillary cell carcinoma patients

Little data is available on prognostic biomarkers and effective treatment options for Kidney Renal Papillary Cell Carcinoma (KIRP) patients, to find potential prognostic biomarkers and new targets was an urgent mission for KIRP therapy. The differentially expressed autophagy-related genes (DEARGs) were screened out according to the RNA sequencing data in The Cancer Genome Atlas database, then identified survival-related DEARGs to establish a prognostic model for survival predicting of KIRP patients. Then we verified the robustness and validity of the prognostic risk model through clinicopathological data. At last, we evaluate the prognostic value of genes that formed the prognostic risk model individually. We analyzed the expression of 232 autophagy-related genes (ARGs) in 289 KIRP and 32 non-tumor tissue cases, and 40 mRNAs were screened out as DEARGs. The functional and pathway enrichment analysis was done and protein–protein interaction network was constructed for all DEARGs. To sift candidate DEARGs associated with KIRP patients’ survival and create an autophagy-related risk prognostic model, univariate and multivariate Cox regression analysis were did separately. Eventually 3 desirable independent prognostic DEARGs (P4HB, NRG1, and BIRC5) were picked out and used for construct the autophagy-related risk model. The accuracy of the prognostic risk model for survival prediction was assessed by Kaplan–Meier plotter, receiver-operator characteristic curve, and clinicopathological correlational analyses. The prognostic value of above 3 genes was verified individually by survival analysis and expression analysis on mRNA and protein level. The autophagy-related prognostic model is accurate and applicable, it can predict OS independently for KIRP patients. Three independent prognostic DEARGs can benefit for facilitate personalized target treatment too.


Background
Renal cell carcinoma (RCC) is the sixth/eighth most common tumor in men and women. About 73,750 new cases of RCC and 14,830 RCC-related deaths happened yearly in the United States [1,2]. Kidney Renal Papillary Cell Carcinoma (KIRP) is the second most commonly diagnosed subtype of RCC which accounting for 15-20% of RCC cases, and the most common subtype is clear cell renal carcinoma (ccRCC) [3,4]. It has been widely assumed that KIRP has a significantly better prognosis than ccRCC in organ-confined stages [5,6]. Most KIRP were manifested as localized disease and treated by partial nephrectomy, but substantial numbers of patients will eventually relapse [7]. About 1/3 of KIRP were manifested as metastatic KIRP (m-KIRP) [8]. Most m-KIRP patients needed systemic treatment eventually. In general, m-KIRP had a worse prognosis than metastatic ccRCC [9,10]. As most KIRP patients had a good prognosis before metastasis, and m-KIRP patients are not so many in contrast to ccRCC, there is little data on the efficacy of available treatment options and few prognostic molecular markers have been discovered. Therefore, there is a need to explore potential prognostic markers and new molecular targets for KIRP therapy.
Traditionally, the TNM stage has been used to evaluate the risk of tumor recurrence for all RCC subtypes. However, it has limited accuracy [11]. Now, some prognostic factors such as grade and pathology stage are also used to evaluate the prognosis. However, these prognostic models were often established for ccRCC only or all RCC subtypes [12,13]. Therefore, there is a need to refine the prognostic risk model of KIRP and build a more accurate approach for managing this second commonest subtype of RCC.
Autophagy is a non-specific, lysosome-mediated degradation. The process is beneficial for cells internally break down, clearance of damaged or superfluous proteins, and recycle cellular components. Autophagy-Related Genes (ARGs) participates in autophagy, for example, ATG7 involved in energy metabolism is an ARG. Lots of researchers proved that autophagy is associated with the progress of RCC [14][15][16]. For example, inhibiting autophagy in RCC increases the efficacy of many therapies [17,18]. However, whether the expression level of ARGs has prognostic value is unknown. Hence, this research utilized ARGs to establish the prognostic risk model of KIRP. In this study, the relevance between differentially expressed autophagy-related genes (DEARGs) and clinicopathological parameters in 321 KIRP patients from TCGA database were examined, and an autophagy-related risk prognostic model was constructed as an independent predictor for overall survival of KIRP patients. We verified our risk score model from several aspects and confirmed it's available, and we hope to provide more helpful guidance for evaluated prognosis and targeted treatment of KIRP with this novel prognostic risk model.

Data acquisition
There are 232 genes presently known associated with autophagy were downloaded from HADb (Human Autophagy Database, http://autophagy.lu/). The RNAseq data and the corresponding clinical data of 289 Kidney Renal Papillary Cell Carcinoma (KIRP) patients and 32 non-tumor samples were obtained from TCGA database (The Cancer Genome Atlas database, https:// www.cancer.gov/about-nci/organization/ccg/research/ structural-genomics/tcga).

Pathway enrichment analysis and Functional annotation for all DEARGs
To reveal the involved pathways and biological function of DEARGs, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis and Gene Ontology (GO) analysis with the clusterProfiler package of R (version 4.0.1), and p-value < 0.05 was used as a strict cutoff.

Protein-protein interaction (PPI) networks construction
The functional protein-protein interaction (PPI) analysis is performed by STRING database (Search Tool for the Retrieval of Interacting Genes, https://string-db.org/) for all DEARGs. The cut-off criterion of interaction score is 0.4. We make use of the Cytoscape software to search hub genes and achieve two-dimensional (2D) visualization of PPI networks.

Identify the prognostic DEARGs
To identify the prognostic DEARGs whose expression profiles had notable correlation with the overall survival (OS) of patients with KIRP, the univariate Cox regression model was constructed. P < 0.05 is the threshold. These selected DEARGs were regarded as candidate genes had a correlation with patients' survival.

Prognostic risk model construction and risk score calculation
The prognostic DEARGs identified by make use of the univariate Cox regression analysis were subjected to a multivariate Cox proportional hazards model to remove the genes that might not be an independent indicator in prognosis predicting. After that, several optimal independent survival-related DEARGs were obtained and the risk score composed of expression value of these genes was established. We calculated the risk score for each patient utilizing the regression coefficients of the individual DEARGs obtained from the multivariate Cox Fig. 1 Identification of differentially expressed autophagy-related genes (DEARGs) from KIRP tissues and non-tumor kidney specimens. a Volcano plot of 242 autophagy-related genes (ARGs). Up-regulated genes are marked red, they are DEARGs which |Log 2 FoldChange| > 1.0 in mRNA level; Down-regulated genes are marked green and they are also DEARGs whose |Log 2 FoldChange| > 1.0. The genes whose |Log 2  The risk score was calculated based on a linear combination of the relative gene expression level multiplied regression coefficients. The regression coefficients are obtained from the multiple Cox analysis and represents the relative weight of the genes. The risk score is a measure of prognostic risk for KIRP patients. Patients were divided into 2 groups by the median risk score as the critical value. High-risk score group had worse prognosis than low-risk score group.

Assessment of prognostic risk model
To verify the robustness and validity of the prognostic risk model, we plotted the survival curves and assessed the differences in the survival rates between high-risk and low-risk groups using the log-rank test. Then, we evaluated the survival prediction accuracy of the prognostic risk model using receiver-operator characteristic (ROC) curve. The area under the curve (AUC) of ROC curve is a discrimination criterion, it ranges from 0.5 to 1.0, the higher the value, the more accurate the model.
To explore whether the autophagy-related prognostic risk model could be an independent predictor of OS not rely on other clinicopathological parameters, we performed cox proportional hazard regression analysis. The association between risk score and clinical traits were explored. To validate the prognostic value of the risk score model, we took age, sex, pathological stage, tumor grade and T classification (lymphatic metastasis excluded) as candidate risk factors for univariate and multivariate Cox regression analyses.

Evaluation of the prognostic value of 3 prognostic-related DEARGs
As described above, to validate the availability of the prognostic risk score model, we compared the survival differences between high-risk group and low-risk group, which grouping is based on risk scores. Afterwards, we studied the association between the expression level of 3 prognostic-related DEARGs and KIRP patients' survival individually. KIRP patients' survival data in TCGA were used for Kaplan-Meier survival analyses.
The expression of 3 independent prognostic-related DEARGs (P4HB, NRG1 and BIRC5) were compared and validated between normal kidney tissues and KIRP tissues in mRNA level and protein level. The mRNA expression of 3 independent prognostic-related DEARGs in kidney and KIRP tissue was analyzed utilizing cancer

Identification of DEARGs
We downloaded the mRNA sequencing data and corresponding clinical data of 289 KIRP tissue samples and 32 normal kidney samples from TCGA database. The gene expression profile and clinical follow-up information of 265 KIRP patients were involved in our subsequent analysis. We extracted expression profile of 232 ARGs. Finally, 40 DEARGs were screened out, involved 31 upregulated ARGs and 9 downregulated ARGs with |log 2 (FoldChange)| > 1 and FDR < 0.05 as filter criteria. The flow chart of the overall procedures in this manuscript was showed in Supplementary  Fig. 1.
In Fig. 1a, the volcano plot exhibited the distribution of all DEARGs. X-axis of the volcano plot is log 2 Fold-Change and Y-axis represents false discovery rate. The fold change patterns of 40 DEARGs in 32 non-tumor tissues and 289 KIRP tissues were showed in a heat map in   Fig. 1c visualized expression of 40 DEARGs between KIRP and normal tissues. In Table 1, we provide a detailed source of information on all DEARGs, including log 2 FoldChange and statistical significance.

PPI network establishment and function annotation of all DEARGs
The interaction of all DEARGs were visualized in Fig. 2a, and there are 12 hub genes arranged in a circle are DEARGs with interaction degree > 15. Biological processes (BP) annotation reminded us that DEARGs had a strong association with enzyme-related process, such as autophagy, regulation of peptidase activity and macroautophagy. In the aspect of the molecular function (MF), DEARGs seems played vital roles in some protein binding related functions, for example, peptidase regulator activity, ubiquitin protein ligase binding and protease binding. Regarding the cellular components (CC), the DEARGs encoded proteins constituted autophagosome membrane, vacuolar membrane, endoplasmic reticulum-Golgi intermediate compartment and so on (Fig. 2b). In Fig. 2c, the mainly pathways that had a positive relation with screened DEARGs was showed, including hepatitis B, pathogenic Escherichia coli infection, human cytomegalovirus infection and Kaposi sarcoma-associated herpesvirus infection. Z-scores of enriched pathways were all > 0, it means that all the pathways were more likely to be enhanced.

Establishment of autophagy-related risk signature
To identify the relationship between the expression of 40 DEARGs and overall survival in KIRP patients, we constructed the univariate Cox proportional hazards model. The results showed that there are 14 DEARGs significantly related to the prognosis of KIRP patients (p < 0.05) (Fig. 3a). In order to raise the robustness, the screened 14 prognostic-related DEARGs were further included in the subsequent multivariate Cox regression analysis. At last, 3 DEARGs (P4HB, NRG1 and BIRC5) were filtered out and used for autophagy-related risk model construction (Fig. 3b). The risk score for each patient was calculated according to the following formula: risk score = (0.8658× expression value of P4HB) +  0.3379× expression value of NRG1) + (1.1201× expression value of BIRC5). Patients were divided into highrisk (n = 132) and low-risk group (n = 133) by the median risk score as the critical value. The risk score was calculated for each patient and list of they belong to low or high-risk group (Supplementary Table 1).

Validation of the risk signature
To measure the accuracy of the autophagy-related risk model to predict the prognosis of KIRP patients, we draw Kaplan-Meier plotter to compare the survival time difference between high-risk and low-risk group. Low-risk group patients had more survival probability (p = 4.406E-05) (Fig. 4a). Next, the ROC curves were employed to determine the predictive performance of the prognostic risk model. In Fig. 4b showed, the AUC value of risk score was 0.923, it was larger than AUC values of other indicators except for the pathologic stage, which confirmed that the autophagyrelated prognostic model is an excellent and independent prognostic predictor comparing with other clinicopathology indicators. The risk scores of all KIRP patients were visualized from small to large (Fig. 4c). With the increase of the risk score, the death number of KIRP patients is more (Fig. 4d). The expression patterns of 3 prognostic-related DEARGs in different risk groups was exhibited in the heatmap (Fig. 4e).

Clinical verification of prognostic model
To determine the relationship between the autophagyrelated prognostic risk model and clinicopathological features in KIPR patients, we put several familiar clinicopathological factors and risk score to do univariate and multivariate cox regression analyses (Fig. 5). Since AUC values are often used to assess the performance of an individual clinicopathological indicator, and the larger of the AUC value, the more accurate of the indicator to predict prognosis. In our study, the AUC values of the clinicopathological features including age and sex to predict OS is less than 0.5 (Fig. 4b), demonstrated that age or sex alone was unable to predict prognosis as an individual indicator. The relationship between risk scores and age/sex was listed in Fig. 5a and Fig. 5b. No difference in risk score was observed between elder patients and younger patients (Fig. 5a). Instead, the AUC values of the clinicopathological features (pathological stage) and risk score is more than 0.9 (Fig. 4b), illustrated that no matter pathological stage or risk score can make a comparatively accurate prediction KIRP patients' prognosis. Risk scores were lower in pathological stage I than in pathological stage II-IV (p = 6.676e-05) (Fig. 5c), and lower in T classification T1-2 than in T3-4 (p = 5.622e-04) (Fig. 5d). The relationship between the expression of 3 genes composed the risk model and 4 clinicopathological features in KIRP are exhibited in Fig. 5a, b, c, d. The raw TCGA data that containing basic information of all KIRP patients was listed in Supplementary Table 2. In Table 2, pathological stage, T classification, and risk score had obviously positive correlation with prognosis of KIRP patients in univariate Cox analysis, in addition, risk score and pathological stage were independent prognostic predictor of KIRP patients in multivariate Cox analysis. All above results demonstrate that the autophagy-related risk signature can be an excellent prognostic predictor ifor KIRP patients.

Validation of the function of 3 prognostic-related DEARGs in KIRP
According to our results, 3 prognostic-related DEARGs including P4HB, NRG1 and BIRC5 were identified to develop the survival-related risk prognostic model. We analyzed the correlation between 3 prognostic-related DEARGs that composed the prognostic model and survival probability of KIRP patients. The results of the Kaplan-Meier analysis showed that the upregulation of P4HB was obviously associated with the low survival probability of KIRP patients. Also, NRG1 or BIRC5 overexpression leads to worse OS (Fig. 6a). To further compare the expression difference of 3 prognostic-related DEARGs between KIRP and normal tissues, we performed a clinical study using cancer microarray database  of Oncomine and human proteome database Human Protein Atlas. The mRNA level of 3 prognostic-related DEARGs was verified and showed in Fig. 6b. The expression trend of 3 prognostic-related DEARGs is in accordance with our previous results obtained from TCGA database which showed in Fig. 1. In Fig. 7a, results of immunohistochemistry (IHC) confirmed the expression of P4HB, NRG1 and BIRC5 protein are stronger in KIRP tissues than normal kidney tissues. Single-gene GSEA of the 3 prognostic-related DEARGs explored the potential roles of 3 prognostic DEARGs in KIRP (Fig. 7b).

Discussion
Autophagy is a dynamic and conserved process that can maintain cellular homeostasis [19]. Many researchers had proved autophagy played an important function in cancer [20][21][22]. Some targeted agents aimed at autophagy had been applied in the clinics in RCC patients [23][24][25]. Therefore, we decided to develop an autophagy-related prognostic risk model for prognosis predicting of KIRP patients.
The autophagy-related prognostic model is comprised of 3 genes, including P4HB, BIRC5, and NRG1. The 3 genes are all closely relevant to clinicopathological features of KIRP patients, the expression of them is either correlated with pathological stage or T classification of KIRP patients (p < 0.05). The risk score obtained according to the prognostic model is also significantly associated with clinicopathological features of KIRP patients. Several assessment methods confirmed the autophagyrelated signature can be an independently prognosis predictor for KIRP patients. Individual assessment of the function of 3 prognostic-related DEARGs in KIRP further proof P4HB, BIRC5, and NRG1 are up-regulated in KIRP, and overexpression of them is associated with worse survival of KIRP patients. All results proved that risk signature constructed through P4HB, BIRC5, and NRG1 to evaluate the prognosis of KIRP patients is clinically practicable. P4HB encoded protein disulfide isomerase, it is an autophagy-related gene, Xie et al proved that P4HB was a novel biomarker for ccRCC diagnosis and prognosis predicting [26]. BIRC5 also called survivin, the encoding products of BIRC5 play a significant role in negative regulation of apoptosis or programmed cell death. Philipp et al demonstrated that BIRC5 is of importance for renal pathophysiology and pathology [27]. NRG1 is a growth factor of the epidermal growth factor family, the relationship between NRG1 and RCC had not been explained clearly, but Sushma et al thought that NRG1 fusion is a low frequency event in most tumor types, including RCC [28]. We found that most of prognostic molecular indicators and therapeutic targets were identified and verified in ccRCC only or all RCC subtypes, there is an urgent need to explore possible prognostic molecular indicators and novel targets for KIRP therapy.
Functional enrichment analysis exhibited that 40 DEARGs were mainly involved in infection related pathway. It is generally acknowledged that infection and inflammation are all correlated with carcinogenesis. Kaymakcalan et al reported that RCC patients treated with mTOR inhibitors had a risk of infection [29], Alexander et al reported that urinary tract infection history is positively associated with RCC development [30]. Hence, the affection of infection to KIRP patients should be carefully assessed and managed.
This research constructed a novel prognostic risk model for prognosis predicting of KIRP patients, and proved it is steady and credible by verification with molecular signature combined with clinical features. Several assessment methods confirmed the prognostic risk model is obviously an independently prognosis predictor for KIRP patients. We believe, apart from traditional clinicopathological features (including pathologic stage, T classification and so on), risk score derived from the autophagy-related genes signature could also be incorporated into the clinical evaluation indicators to better predict clinical outcomes. Individual assessment of the function of 3 prognostic-related DEARGs in KIRP further proved P4HB, BIRC5, and NRG1 all play significant roles in KIRP. The 3 prognostic-related DEARGs can benefit personalized target therapy also. All results proved that the risk score calculated according to expression level of P4HB, BIRC5, and NRG1 to evaluate the prognosis of KIRP patients is reliable.

Conclusions
This research analyzed mRNA sequencing data of 289 KIRP tissue specimens and 32 non-tumor specimens and assessed of 232 ARGs' expression difference in the two groups. We screened out 9 down-regulated DEAR Gs and 31 up-regulated DEARGs in KIRP with the threshold of |log 2 FC| > 1.0 and P < 0.05. From 40 DEAR Gs, 3 prognostic DEARGs (P4HB, NRG1, BIRC5) were determined to establish a prognostic risk model, and the risk score was calculated according to expression of the