Prognostic risk signature based on the expression of three m6A RNA methylation regulatory genes in kidney renal papillary cell carcinoma

Fourteen out of twenty m6A RNA methylation regulatory genes are differentially expressed in KIRP tissues

We analyzed the expression levels of 20 m6A RNA methylation regulating proteins in KIRP (n = 289) and normal samples (n = 32) from the TCGA database. The heatmap showed that the expression of 14 m6A methylation regulators (IGF2BP3, IGF2BP1, HNRNPC, YTHDF2, KIAA1429, YTHDF3, METTL14, ZC3H13, ALKBH5, IGF2BP2, RBM15B, YTHDF1, RBMX and HNRNPA2B1) were differentially expressed in KIRP tissues compared to normal kidney tissues (Figure 1A). We observed that IGF2BP3, RBMX, YTHDF1, IGF2BP2, HNRNPA2B1, RBM15B and HNRNPC were significantly upregulated, and YTHDF2, IGF2BP1, METTL14, ALKBH5, KIAA1429, YTHDF3 and ZC3H13 were significantly downregulated in the KIRP tissues compared to the normal kidney tissue samples (Figure 1B). The expression of WTAP, RBM15, YTHDC2, FTO, METTL3 and YTHDC1 was similar in KIRP and normal kidney tissue samples (Figure 1A).

PPI network and correlation analysis between m6A RNA methylation regulators

We then downloaded the data for the 14 differentially expressed m6A RNA methylation regulators from the STRING database and constructed a protein–protein interaction (PPI) network using Cytoscape. PPI network analysis showed that KIAA1429, HNRNPC, METTL14, HNRNPA2B1 and ALKBH5 were the hub genes (Figure 2A). The interaction between KIAA1429 and YTHDF3 (r = 0.79) was most significant among all the m6A RNA methylation regulators (Figure 2B).

Identification of a prognostic risk signature based on three m6A RNA methylation regulators in the training cohort of KIRP patients

Univariate Cox regression analysis of the transcriptome data from the TCGA-KIRP dataset showed that the expression of four m6A RNA methylation regulators (IGF2BP3, KIAA1429, YTHDF3 and HNRNPC) was significantly associated with the overall survival (OS) of KIRP patients (P < 0.05; Figure 3A). Furthermore, based on the results of the least absolute shrinkage and selection operator (LASSO) Cox regression analysis, we constructed a prognostic risk signature with three genes, including IGF2BP3, KIAA1429 and HNRNPC (Figure 3B, 3C).

Next, we analyzed the efficacy of the risk signature in predicting the prognosis of KIRP patients in the training cohort. We calculated the risk score for each patient in the training cohort and divided them into high-risk (n=71) and low-risk (n=71) groups according to the median risk score. Kaplan-Meier survival curve analysis showed that overall survival time was significantly shorter for KIRP patients in the high-risk group compared to the KIRP patients in the low-risk group (log-rank test P < 0.001; Figure 3D). The 3-yr OS rates for the high- and low-risk group were 81.3% and 96.3%, respectively. The 5-yr OS rates for the high- and low-risk group were 62.8%and 96.3%, respectively. The ROC curve analysis demonstrated significant prognostic predictive value for the risk signature in determining the 1-yr, 3-yr and 5-yr survival rates of KIRP patients (1-year AUC = 0.864, 3-year AUC = 0.903, 5-year AUC = 0.714; Figure 3E). The risk score distribution of the high- and low-risk group patients in the training cohort showed that the survival rates were significantly higher for the low-risk group compared to the high-risk group (Figure 3F, 3G). The heatmap showed higher expression levels of the three risk-related m6A RNA methylation regulators (IGF2BP3, KIAA1429 and HNRNPC) in the high- risk KIRP patients compared to the low-risk KIRP patient group (Figure 3H).

Validation of the prognostic risk signature in the testing and entire TCGA-KIRP cohorts

We validated the accuracy and robustness of prognostic risk signature in the testing cohort and the entire TCGA-KIRP cohort. We calculated the prognostic risk scores of patients in both the testing (n=94) and the entire TCGA-KIRP (n=237) cohorts based on the prognostic risk signature and stratified the patients into high-risk (testing cohort: 47; entire TCGA cohort: 119) and low-risk (testing cohort: 47; entire TCGA cohort: 118) groups based on the median cut-off value. The detailed clinicopathological features of all the TCGA-KIRP patients are listed in Table 1.

Kaplan-Meier survival curve analysis showed that the overall survival of high-risk patients were significantly shorter compared to the low-risk patients in both the testing cohort (P < 0.05; Figure 4A) and the entire TCGA-KIRP cohort (P < 0.001; Figure 4B). In the testing cohort, the 3-year and 5-year survival rates were shorter for the high-risk patients compared to the low-risk groups (3-year: 73.3% vs. 97.7%; 5-year: 61.9% vs. 82.5%). Similarly, in the entire TCGA-KIRP cohort, the 3-year and 5-year survival rates in high-risk group were lower than those in the low-risk group (3-year: 77.7 vs. 97.0%; 5-year: 62.0 vs. 88.6%). ROC curve analysis showed that the AUC values for the 1-, 3- and 5-year survival in the testing cohort were 0.988, 0.853 and 0.712, respectively (Figure 4C). The AUC values for the 1-, 3- and 5-year survival in the entire TCGA-KIRP cohort were 0.925, 0.869 and 0.708, respectively (Figure 4D). The risk score distribution, survival status and the risk gene expression in the testing cohort and the entire TCGA-KIRP cohort are shown in Figures 4E and 4F. Overall, our results showed that the prognostic risk signature accurately predicted the survival outcomes of KIRP patients.

Independent prognostic value of the risk signature

To determine whether the risk signature can be used as an independent prognostic factor, Next, we performed univariate and multivariate Cox regression analyses to determine the independent prognostic significance of the risk score and the relevant clinicopathological factors, including age, gender, AJCC stage, T stage, N stage, and M stage. In the TCGA training cohort, univariate analyses showed that AJCC stage, T stage, N stage, M stage and risk score were significantly associated with OS. Subsequently, multivariate analyses showed that AJCC stage, T stage, M stage and risk score were significantly associated with OS (Table 2). Similar results were obtained for both the testing and the entire TCGA-KIRP cohorts (Table 2). These results demonstrate that the risk score calculated based on the prognostic risk signature is an independent prognostic factor in KIRP patients.

Consensus clustering of KIRP patients based on the expression of the three prognostic risk-related m6A RNA methylation regulators

We then performed consensus clustering based on the expression levels of the three m6A RNA methylation regulators in the TCGA-KIRP dataset. We chose K = 2 as the most optimal clustering of the TCGA-KIRP patients because the clustering was suboptimal when divided into more than 2 clusters (Figure 5A–5D). Principal component analysis (PCA) also divided the TCGA-KIRP patients into two clusters (clusters 1 and 2) based on their transcriptional profiles (Figure 5E). Subsequently, Kaplan-Meier survival curve analysis showed that the OS was significantly shorter for the KIRP patients in cluster 2 compared to those in cluster 1 (Figure 5F). We then analyzed the correlations between the two clusters and their corresponding clinicopathological features. KIRP patients in clusters 1 and 2 showed significant differences in AJCC stage (P < 0.05), N stage (P < 0.01) and survival status (P < 0.01), but did not show any significant differences in age, gender, T stage and M stage (Figure 5G). Moreover, the expression of m6A RNA methylation modulators was significantly higher in the cluster 2 KIRP patients compared to the cluster 1 KIRP patients.

Relationship between prognostic risk signature and clinicopathological parameters of KIRP patients

Next, we analyzed the association between the prognostic risk signature and the clinicopathological parameters. The heatmap showed that the expression levels of the three risk-related m6A RNA methylation regulators correlated with the clinicopathological variables in the high- and low-risk groups. We observed significant differences between the high- and low-risk groups in regard to T stage (P < 0.01), AJCC stage (P < 0.01), gender (P < 0.01) and survival status (P < 0.001) (Figure 6A). Moreover, advanced-stage tumors significantly associated with the high-risk group, whereas, the early-stage tumors correlated with the low-risk group (Figure 6B–6E).

We further analyzed the prognostic value of our risk score model in different subgroups of KIRP patients that were stratified based on clinicopathological parameters. We observed significantly shorter overall survival rate in the male patients (P =2.776e-04) and those with age ≤ 65 (P = 5.846e-04), AJCC stage III/IV (P = 1.775e-02), T1-2 stage (P = 2.829e-02), T3-4 stage (P= 2.591e-02), N0 stage (P = 8.657e-03), N1-2 stage (P=2.839e-02) and M0 stage (P =4.862e-04) in the high-risk group compared to those in the low-risk group (Figure 7). However, OS rates were similar in female patients (P = 2.314e-1), and those with age > 65 (P = 1.494e-1), AJCC stage I/II (P = 1.324e-1), and M1 stage (P = 8.386e-1) in both the high- and low-risk groups.

Validation of the three prognostic signature-related genes

Subsequently, we analyzed the correlation between the expression of each of the three prognostic risk signature genes and the clinicopathological features of KIRP patients in the TCGA cohort. We observed differential expression of the 3 risk signature-related genes across various clinicopathological parameters (Table 3). Higher expression of IGF2BP3, KIAA1429 and HNRNPC correlated with higher AJCC, T, N and M stages in the KIRP patients (Figure 8). Moreover, we observed age-dependent differences in the expression of IGF2BP3, and gender-related differences in the expression of KIAA1429 (Table 3).

We further analyzed the relationship between the risk signature-related genes and the OS and DFS rates of KIRP patients in the GEPIA database. Kaplan-Meier survival curves and log-rank test showed that the OS rate of KIRP patients with higher expression of IGF2BP3 (P = 6.4e-05), KIAA1429 (P = 0.05) and HNRNPC (P = 7.2e-04) was significantly shorter compared to those with lower expression of the three risk signature-related genes (Figure 9A). Moreover, higher expression of IGF2BP3 (P = 0.018) and KIAA1429 (P = 9.1e-04) correlated with significantly shorter DFS rate (Figure 9B). These data are consistent with our previous findings in the TCGA cohort of KIRP patients. However, we did not find any significant differences in the DFS rates of KIRP patients with differential expression (high or low) of HNRNPC (Figure 9B). Overall, our results demonstrate that the three m6A RNA methylation regulators are potential prognostic biomarkers that can accurately predict survival outcomes of KIRP patients.

Patient data

In this study, we obtained RNA-sequencing data of 289 KIRP and 32 normal kidney samples from the TCGA database (https://cancergenome.nih.gov/). We downloaded the corresponding clinical data regarding gender, age, clinicopathological parameters, and survival using the GDC data transfer tool (https://portal.gdc.cancer.gov/). We excluded 14 KIRP patients from the study because their survival time was less than 30 days. We also excluded 38 KIRP patients because of incomplete data. Ethics approval and informed consent were not required for this study because we used publicly available data from TCGA.

Identification of differentially expressed m6A RNA methylation regulators

We systematically analyzed the expression of 20 m6A-related genes, including seven m6A writers (METTL3, METTL14, WTAP, KIAA1429, RBM15, RBM15B and ZC3H13), two m6A erasers (FTO and ALKBH5), and eleven m6A readers (IGF2BP1, IGF2BP2, IGF2BP3, YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2, HNRNPC, HNRNPA2B1 and RBMX) in KIRP and non-cancerous kidney samples using the LIMMA package (version 3.30.3; http://bioconductor.org/packages/release/bioc/html/limma.html) from the R software (version 3.6.2; https://cran.r-project.org/bin/windows/base/) and identified differentially expressed genes using P-value < 0.05 as a threshold parameter.

PPI network and Pearson correlation analysis

We constructed a PPI network between the 14 differentially expressed m6A RNA methylation regulatory genes by downloading the data from the Search Tool for the Retrieval of Interacting Genes (STRING; version 11.0, http://string-db.org) online database [17] and constructing the network with the Cytoscape software (version 3.7.2) between the differentially expressed m6A RNA methylation regulators [35]. Pearson correlation analysis was performed to determine positive or negative association between different m6A RNA methylation regulators based on the Pearson’s correlation coefficient (r) value between −1 and +1.

Construction of the prognostic risk model

We randomly assigned the 237 KIRP patients to a training cohort (n = 143) and a testing cohort (n = 94). We then performed univariate Cox regression analysis to determine the correlation between differentially expressed m6A RNA methylation regulatory genes and OS of KIRP patients in the training set. We selected the genes showing significant correlation with OS (P < 0.05) as prognosis-related genes. Subsequently, we performed the least absolute shrinkage and selection operator (LASSO) Cox regression analysis and identified 3 potential genes to develop the prognostic risk signature. The risk score (RS) of the patients was estimated using the following formula:

$\text{RiskScore}={\displaystyle {\sum }_{i=1}^{n}Coef(i)\times x(i)}$

wherein, Coef(i) and x(i) represent the estimated regression coefficient and the expression value of each target gene by LASSO analysis, respectively. The risk score formula for the 3-gene prognostic risk signature was as follows: risk score = (0.6715 × expression value of IGF2BP3) + (0.2675 × expression value of KIAA1429) + (0.0109 × expression value of HNRNPC). We then divided the KIRP patients in the training cohort based on their risk scores into high- and low-risk groups by using the median risk score as the cut-off value. We then analyzed the survival parameters of the two groups using the Kaplan-Meier survival curve and two-sided log-rank test. We also used receiver operating characteristic (ROC) curves and area under the ROC curve (AUC) values to evaluate the prediction accuracy of our prognostic model. We also performed univariate Cox regression analysis to evaluate the prognostic prediction ability of other clinicopathological factors such as age, gender, AJCC (The American Joint Committee on Cancer) stages, and TNM (Tumor, Node, Metastasis) stages. Multivariate Cox regression analysis was used to assess if the risk score was an independent prognostic factor. The robustness and reliability of the prognostic risk score was then evaluated in the testing cohort (n=94) and the entire TCGA-KIRP cohort (n=237).

Consensus clustering analysis

We used the Consensus Clusterplus R package to identify different KIRP patient subgroups based on the cumulative distribution function (CDF), delta area, tracking plot and consensus matrix. We evaluated k=2 to 9 potential subgroups for class discovery and clustering validation (50 iterations and 80% resampling rate). Principal component analysis (PCA) was used to assess the signature-related genes expression differences in the two clusters. Kaplan-Meier survival curves and log-ranktest were used to determine the OS rates of KIRP patients in the two clusters. Chi-square test was used to evaluate the differences of clinicopathological characteristics between the two clusters.

Validation of the clinical utility of the prognostic risk model

Using the Kaplan-Meier (K-M) survival curves, we tested the ability of the prognostic risk signature to predict the survival outcomes of TCGA-KIRP patients stratified by various clinicopathological characteristics, including age (≤ 65 and > 65), gender (female and male), AJCC stage(I/II and III/IV), T stage (T1-2 and T3-4), N stage (N0 and N1-2) and M stage (M0 and M1). The relationship between the expression of the three individual risk signature genes and the clinicopathological characteristics was evaluated using the Student’s t-test. We also analyzed the relationship between the risk signature-related genes and the survival parameters (OS and DFS) of KIRP patients in the GEPIA database (http://gepia.cancer-pku.cn/).

Statistical analysis

All statistical analyses were performed using the R software (version 3.6.2). Student’st-test was applied to examine the differences among variables. Pearson’s or Spearman correlation analyses were used to determine the association between various parameters. Kaplan-Meier curve analysis with log-rank test was used to analyze survival rates between different patient subgroups. Chi-square tests were performed to compare categorical variables. Univariate and multivariate Cox regression analyses identify prognostic significance of various clinicopathological characteristics and the prognostic risk score. ROC curve analysis was used to evaluate the efficacy of the prognostic risk model to discriminate between high- and low-risk KIRP patients as well as cluster 1 and cluster 2 KIRP patients based on the AUC values. An AUC value of 1.0 denotes perfect prognostic prediction, whereas an AUC value below 0.5 denotes poor prediction. P < 0.05 was considered statistically significant.