Identification of a three-m6A related gene risk score model as a potential prognostic biomarker in clear cell renal cell carcinoma

Collection of primary data

The TCGA-KIRC mRNA FPKM transcriptome profiling data (with 539 tumor samples and 72 normal samples) and associated clinical information (N = 491) were downloaded from TCGA (http://cancergenome.nih.gov/), we then extracted the expression data of these seventeen m6A genes for following operations.

Consensus clustering of m6A genes

Identification of distinct clusters of ccRCC patients using consensus clustering was performed primarily, by which we decided to divide patients into k (k = 2–9) clusters according to the expression of m6A genes, moreover, we did survival analyses of these k clusters to see whether the clustering is reasonable.

Data processing and DEGs analysis

DEGs were filtered by comparing the expression level of ccRCC samples and normal samples in TCGA dataset utilizing the wilcox test in Limma package in R software with the cutoff criteria of p-value < 0.05. A violin plot was depicted to clearly display the differential expression of m6A genes.

Selection of potential prognostic-related genes

We did univariate cox analysis of all DEGs by survival package in R software. The DEGs with a p-value less than 0.05 were considered as prognostic-associated genes and selected as candidate genes for further procedure.

Construction of RS

In order to gain better prediction of m6A genes and ccRCC, we compared three analysis models (Least Absolute Shrinkage and Selection Operator (LASSO), Ridge and Elastic Net analysis) to construct better RS. Finally, we implemented LASSO COX regression analysis to candidate genes from TCGA (which was used as a training set) utilizing the survival and glmnet packages in R. Finally, genes and their coefficients were determined. The RS was calculated by the formula: $$\mathrm{RS}={\displaystyle \sum _{i=1}^n\mathrm{coefi}}\times x_i$$

In which the coefi is the coefficient, and x_i is the expression value of each selected gene. This formula was used to calculate a RS for each patient in both the training (TCGA) dataset.

OS analysis of RS

To investigate the relationship between the RS and 5 year OS, we used survival package in R to map a Kaplan–Meier plot (K–M plot) by TCGA data, patients were divided into high and low risk groups by the median value of RS, p < 0.05 was applied as cut-off standard. In addition, we plotted a ROC curve by survivalROC package, moreover, the AUC, was calculated to predict the accuracy of the survival analysis.

Cox regression analysis of RS

We ran univariate and multivariate cox regression analyses by survival and forestplot packages to inspect the independency of RS as a prognostic factor. Except for RS, we analyzed the impact of age (which was classified into two groups: >60 and ≤60), gender, T stage, M stage, grade and stage. All the ambiguous value are excluded, including Mx, Tx N/A and unknown (we deleted all the N stage data about TCGA patients because the large amount of Nx). Each factor with a p value < 0.05 in univariate analysis was then analyzed by multivariate cox regression, those with a p value < 0.05 and HR > 1 were viewed as independent prognostic factors.

Validation of the RS

For the purpose of RS validation, we downloaded another dataset, E-MTAB-1980 (N = 101) from ArrayExpress (https://www.ebi.ac.uk/arrayexpress/). However, the expression data from this validation set was log2 processed, thereby, we reversely normalized the data so that it is in consistent with that of TCGA. The RS formula was then applied to each sample in validation set, likewise, the patients were separated to high and low risk groups in term of median RS value. Eventually, a K–M plot and an ROC curve were constructed for verification.

Validation by quantitative real-time PCR from clinical specimens of ccRCC

To further validate our finding from TCGA and ArrayExpress database, we performed quantitative real-time PCR (qRT-PCR) to detect expression of HNRNPA2B1, METTL3 and METTL14 in 16 pairs of ccRCC samples and matched adjacent normal kidney tissues (n = 32). All patients received radical nephrectomy from January 2017 to March 2018 in Shengjing Hospital of China Medical University and were pathologically diagnosed as ccRCC. This study was approved by Ethics Committee of Shengjing Hospital (No. 2019PS050J), all patients provided written scientific ethics consents.

RNA was extracted from frozen samples using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instructions. RNA purity and concentrations were identified, then reverse transcription and qRT-PCR were performed. Detailed procedures had been described in our previous paper (Chen et al., 2016). U6 acted as internal control to normalize the results. The relative quantification equation (RQ = 2^−ΔΔCt) was used to calculate the relative expression of RNA.

Clustering of m6A genes

The overall workflow of the process that we used to develop and validate the RS to predict prognostic outcomes was depicted in Fig. 1. Based on the expression similarity of m6A genes, k = 2 seemed to be an adequate selection with clustering stability increasing from k = 2 to 9 in the TCGA datasets, which we named as cluster1 and cluster2 (Figs. 2A–2C). Significant difference of OS was noticed between these two clusters (p = 0.033) (Fig. 2D), we believed that the results may reveal a potential correlation of m6A genes and ccRCC patients, which led to the following exploration.

DEGs and prognostic-related genes screening

From 17 m6A genes, a total of 12 DEGs with a p value < 0.05 were identified. Including YTHDF2, KIAA1429, ALKBH5, HNRNPA2B1, WTAP, METTL3, YTHDF3, FTO, RBM15, YTHDC2, RBM15B and METTL14. The expression level of all 17 genes were illustrated in Fig. 3A. Additionally, univariate cox regression analysis of the DEGs was operated, seven genes (YTHDF2, KIAA1429, HNRNPA2B1, YTHDF3, METTL3, FTO and METTL14) were considered as prognostic-related genes, the results was shown in Fig. 3B.

Formula and OS analysis of the RS

After running the LASSO COX regression analysis of the seven genes, three genes (METTL14, HNRNPA2B1 and METTL3) were filtered and selected to build the prediction model, as is demonstrated in (Figs. 4A and 4B). The total RS was imputed as follows: (−0.413491830046862 × expression level of METTL14) + (0.0149601455134957 × expression level of HNRNPA2B1) + (0.0260335348056223 × expression level of METTL3). The 5-year OS and ROC curve were presented in Figs. 4C and 4D while we found significant differences between the high and low risk groups (p = 3.18E−11), the time-dependent AUC further proved the accuracy of the result (AUC = 0.712).

RS and other clinicopathological characteristics

Age, grade, stage, T stage, M stage and RS showed statistical significance after the univariate analysis (all p value < 0.001), which is illustrated in Fig. 5A, these six factors were then involved in multivariate analysis when we observed age (p = 0.006), grade (p = 0.010), stage (p = 0.018) and RS (p < 0.001) remaining to be significant (Fig. 5B). Besides, all four factors were with an HR that is bigger than one.

Performance of RS in validation set

OS and ROC curve constructed by E-MTAB-1980 dataset (Figs. 6A and 6B) displayed similar results as that in the training set. The p value = 8.659E−03 and AUC = 0.7 all meet the standard of significance. The result confirmed that the RS attained in our study may be a potential prognostic biomarker.

Comparison of LASSO, ridge and elastic net model

In addition to LASSO COX regression analysis, Ridge and Elastic Net analysis were also applied two different RS were obtained. We then did survival and cox analysis via these two RSs. We got meaningful results by all three models in training cohort (Figs. S1–S3). Nevertheless, the Ridge model showed no significance for prognosis in the validation cohort (Fig. S4C). Despite Elastic Net model generated smaller p value to LASSO model in the training cohort (Figs. S1A and S1B), but got higher p value than LASSO model in the validation cohort (Figs. S4A and S4B). Moreover, its AUC of the ROC curve in the validation dataset was 0.686, which did not meet the cut-off standard that AUC ≥ 0.7(Fig. S5). Therefore, we finally chose lasso model to construct RS.

Unicox survival analysis of m6A genes and RS in validation cohort

Forest plot of m6A related genes and RS was mapped (Fig. S6), although all RS related genes showed weak significance in the validation cohort (METTL14, HR = 0.986 (0.9775–0.9945); METTL3, HR = 0.9993 (0.9980–1.0007); HNRNPA2B1, HR = 1.0006 (1.0004–1.0008)), the RS we obtained still had comparative advantage over other genes (RS p < 0.001, HR = 1.0236 (1.0118–1.0355)).

Validation based on qRT-PCR

To further validate the expression pattern of the three differentially expressed genes (METTL14, HNRNPA2B1 and METTL3) from our RS model. The qRT-PCR was performed on 16 paired ccRCC samples and matched adjacent normal kidney tissues. The expression level of METTL3 was significantly higher in ccRCC specimens than those in normal tissues (p = 0.0095) (Fig. 7A). Inversely, the expression of METTL14 and HNRNPA2B1 were relatively lower in tumor tissues (p = 0.0016, p = 0.0163 respectively) (Figs. 7B and 7C). These results were consistent with our finding based on TCGA database.