Predictive value of m1A-related genes in kidney renal clear cell carcinoma

Abstract

Background：Kidney renal clear cell carcinoma (KIRC) is a prevalent type of renal malignancy characterized by high mortality rates and poor response to treatment. N1-methyladenosine (m1A) is a type of RNA methylation modification that has received considerable attention due to its crucial roles in various biological functions. With the advancement of genomics and molecular biology, m1A-related genes (m1A RGs) have been confirmed to be intimately connected with the development and occurrence of various tumors. Nevertheless, the role of m1A RGs in KIRC remains poorly understood.

Methods：This study aims to investigate the prognostic significance of 10 major m1A RGs in KIRC patients, utilizing data from The Cancer Genome Atlas (TCGA) dataset. A prognostic model was constructed using Lasso regression analysis, and risk scores were calculated. KIRC patients were classified into high- and low-risk groups based on the median of the average risk score. The prognostic value of the model was evaluated using two independent datasets, GSE537574 and GSE265745, by assessing the sensitivity and specificity using Kaplan-Meier survival analysis and receiver operating characteristic curves. Additionally, gene set enrichment analysis was conducted to explore the possible biological behavior and pathways of m1A RGs. Ultimately, 5 m1A RGs were identified to construct the prognostic model. Furthermore, nomogram and decision curve analyses were performed to evaluate the model's predictive performance and clinical application value.

Results：Our study demonstrates that the expression of m1A RGs might serve as a prognostic biomarker for KIRC patients and provides a new perspective for cancer prognosis screening in clinical practice.

Background

Kidney cancer is among the most prevalent malignancies affecting the urinary system. Its incidence has increased in recent years, primarily due to an aging population and smoking[1, 2]. In 2022, there were 71,676 new cases reported in the United States, representing 3.02% of all cancer cases. It ranks seventh in the incidence of cancer in men and tenth in women, following prostate and bladder cancer [3]. Kidney renal clear cell carcinoma (KIRC) stands as the predominant subtype of kidney cancer, constituting an estimated 80% of all cases. Regrettably, KIRC is also among the most aggressive subtypes with the poorest prognosis [4]. Although KIRC could be detected early and treated successfully with surgery, early symptoms are not always apparent. Most KIRC patients die within 1-2 years of diagnosis, 30% have distant metastases, and KIRC is insensitive to radiotherapy and chemotherapy [5]. Although immunotherapy for different targets offers hope for these patients, the effective response rate to treatment is limited. Thus, there is an urgent need for additional research on relevant prognostic markers and potential treatment targets to improve patient outcomes in KIRC.

With the development of advanced sequencing technologies, RNA modifications in the progression of various types of cancer have received considerable attention in recent years. RNA methylation is currently the most extensively researched among the more than 150 known RNA modifications [6, 7]. RNA methylation is critical for controlling gene expression and is engaged in carcinogenesis through regulating cancer cell proliferation, invasion, metastasis, and treatment resistance, mainly at the post-transcriptional level [8]. N1-methyladenosine (m1A) is an essential post-transcriptional RNA modification discovered in the 1960s [9]. The regulation of m1A modification status is predominantly achieved by three categories of genes, namely "writers," "readers," and "erasers" [10-12]. The "writers" refer to methyltransferases, including TRMT10C, TRMT61B, and TRMT6/61A, which regulate the level of m1A and influence translation efficiency [13]. Conversely, the "readers," such as YTHDF1, YTHDF2, YTHDF3, and YTHDC1, recognize the methylation mark and engage in RNA translation and degradation processes [10]. Demethylases, including ALKBH1 and ALKBH3, also known as "erasers," remove the m1A modification from single-stranded DNA and RNA [11]. The function of m1A RGs in the occurrence and development of tumors is becoming progressively clearer. TRMT10C silencing inhibits the malignant characteristics of ovarian and cervical cancer cells [14], whereas TRMT6 overexpression is related to an unfavorable prognosis in hepatocellular carcinoma patients [15]. Additionally, ALKBH1 regulates the mTOR and ErbB signaling pathways to promote pancreatic cancer progression [16]. However, the role of m1A RGs in the development of KIRC is poorly understood. Therefore, additional research is required to explore the possible role of m1A in KIRC.

In this study, we utilized the TCGA-KIRC dataset to investigate the connection between m1A RGs and the prognosis of KIRC patients. Two expression profile datasets (GSE537574 and GSE265745) were acquired from the Gene Expression Omnibus (GEO) database for further validation. Subsequently, we established 5 prognostic models for m1A-related genes and constructed a nomogram to assist clinicians in predicting clinical outcomes and tailoring personalized treatments. Finally, we analyzed the expression of m1A RGs and their potential roles in cancer prognosis in KIRC through immunohistochemical (IHC) analysis using the Human Protein Atlas (HPA) database. 

Results

Technology Roadmap (Figure 1)

Table 1: List of Kidney Renal Clear Cell Carcinoma Information.

Development of a prognostic model for m1A RGs and evaluation of clinical relevance

To evaluate the prognostic significance of 10 m1A RGs in the TCGA-KIRC dataset, we constructed a prognostic model for these genes employing Lasso regression analysis (Figure 2A). To visualize the regression results, we also created a Lasso variable trajectory graph (Figure 2B, Table 2). The trajectory graph indicates that our m1A RGs prognostic model was effective. We used a risk factor map to visualize the sample groupings in the prognostic model (Figure 2C). The risk factor plot includes risk grouping, survival outcome, and a heat map. Risk grouping involves grouping risk scores obtained from predicting the model for the dataset samples by the median. The survival outcome presents clinical samples' survival time and outcome in the TCGA-KIRC dataset through point plots. Finally, a heat map is employed to visualize the expression of m1A RGs for the samples in the prognostic model. To further elucidate the disparities among the high- and low-risk groups within the TCGA-KIRC dataset, we conducted a Mann-Whitney U test for the 10 m1A RGs and presented the results in subgroup comparison plots (Figure 2D). The analysis results indicate that the expression levels of the 7 m1A RGs (ALKBH1, ALKBH3, TRMT10C, TRMT6, TRMT61B, YTHDC1, and YTHDF2) exhibit significant variation among diverse groups. Specifically, TRMT6 (P-value < 0.001) expression is notably increased in the high-risk group, while TRMT10C (P-value < 0.001), TRMT61B (P-value < 0.001), ALKBH1 (P-value < 0.001), YTHDF2 (P-value < 0.001), and YTHDC1 (P-value < 0.001) are significantly decreased. The expression level of ALKBH3 (P-value > 0.05) shows no statistically significant variation between the two groups.

Table 2: Expression Genes of m1A between high- and low-risk groups.

Functional enrichment analysis(GO) of m1A RGs

To explore the interaction between BP, MF, and CC for the 7 m1A RGs in the TCGA-KIRC dataset, we performed a GO analysis on these genes (Table 3). Statistical significance was considered for entries with an enrichment P-value < 0.05 and FDR q-value < 0.05. The results revealed that the 7 m1A RGs in KIRC patient samples were significantly enriched in various biological processes, such as RNA modification, mRNA methylation, tRNA modification, mRNA modification, and tRNA processing, as well as cellular components such as methyltransferase complex, ribonuclease P complex, nuclear euchromatin, mitochondrial matrix, and endoribonuclease complex. Additionally，the molecular functions enriched by these RGs included catalytic activities acting on tRNA and RNA, as well as tRNA, RNA, and mRNA methyltransferase activities. The results were presented using bubble plots (Figure 3A) and circular network diagrams for the BP pathway, CC pathway, and MF pathway (Figure 3B-D).

Table 3: GO enrichment analysis results of m1A RGs.

GSEA of DEGs between high- and low-risk groups

To investigate the effect of gene expression levels on KIRC development in distinct groups of the m1A RGs prognosis model in the TCGA-KIRC dataset, we performed a GSEA on the expression of all genes and their participation in BP, CC, and MF between the two groups. We used a significant enrichment screening criterion of P-value < 0.05 and FDR q-value < 0.05 (Figure 4A). The results reveal that genes in both groups of m1A RGs prognosis model in the TCGA-KIRC dataset were significantly enriched in pathways such as IL-5 pathway (Figure4B), IL-6-7 pathway (Figure 4C), IL-27 pathway (Figure 4D), P130Cas linkage to MAPK signaling for integrins (Figure 4E), and Grb2-Sos provides linkage to MAPK signaling for integrins (Figure 4F) (Figure 4B-F, Table 4).

Table 4: GSEA results of m1A RGs Prognosis Model high- and low-risk groups genes.

GSVA of DEGs between high- and low-risk groups

We used GSVA on the DEGs in the TCGA-KIRC dataset to investigate the variation in hallmark gene sets between the two groups (Table 5). Our analysis revealed 32 hallmark gene sets that exhibited significant differences (P-value < 0.05) between the two groups. Of these sets, 17 pathways exhibited significantly higher enrichment scores in the high-risk group (P-value < 0.05), while 15 pathways exhibited significantly higher enrichment scores in the low-risk group (P-value < 0.05). Based on the GSVA results, we analyzed the differential expression of the 32 Hallmark pathways between different groups in the TCGA-KIRC dataset and used the R package "pheatmap" to display the specific differential analysis results in a heatmap (Figure 5A). We also analyzed the degree of group differences of the 32 Hallmark pathways between different groups in the TCGA-KIRC dataset by using the Mann-Whitney U test and presented the results using a group comparison plot (Figure 5B), which demonstrated that the differential expression of the 32 Hallmark pathways between different groups in the TCGA-KIRC dataset had statistically significant differences (P-value < 0.05).

Table 5: GSVA results of m1A RGs Prognosis Model high- and low-risk groups genes.

Comparison of subgroups of m1A RGs prognostic model in patients with KIRC, PPI network and clinical analysis

To analyze the expression of 7 m1A RGs (ALKBH1, ALKBH3, TRMT10C, TRMT6, TRMT61B, YTHDC1, YTHDF2) in KIRC patients with KIRC tumor tissue (subgroup: KIRC) and healthy kidney tissue (subgroup: Normal), we utilized the Mann-Whitney U test to examine the expression of these genes in subgroup samples from the TCGA-KIRC dataset, GSE53757 dataset, and GSE26574 dataset. The results were presented using group comparison plots (Figure 6A-C). Among the 7 m1A RGs in FPKM format from the TCGA-KIRC dataset, YCHDC1(P-value > 0.05) expression did not exhibit any statistically significant differences between the KIRC and Normal groups. The KIRC group showed significant upregulation in ALKBH1 (P-value < 0.01) and ALKBH3 (P-value < 0.001) expression while TRMT10C (P-value < 0.001), TRMT6 (P-value < 0.001), TRMT61B (P-value < 0.01), and YTHDF2 (P-value < 0.05) expression was significantly downregulated (Figure 6A). In the GSE53757 dataset, TRMT6 and YTHDF2 expression among the 7 m1A RGs did not show any statistically significant differences among different subgroups (P-value > 0.05). The KIRC group exhibited significant upregulation in ALKBH3 (P-value < 0.001) and YTHDC1 (P-value < 0.001) expression while ALKBH1 (P-value < 0.001), TRMT10C (P-value < 0.001), and TRMT61B (P-value < 0.001) expression was significantly downregulated (Figure 6B). In the GSE26574 dataset, TRMT61B, YTHDC1, and YTHDF2 expression did not exhibit any statistically significant differences (P-value > 0.05), whereas TRMT6 (P-value < 0.05) expression was significantly upregulated in the KIRC group and ALKBH1 (P-value < 0.05) and ALKBH3 (P-value < 0.05) expression was significantly downregulated (Figure 6C).

To further investigate the protein-protein interactions of the 7 m1A RGs in the TCGA-KIRC dataset, we constructed a PPI network using the STRING database and visualized it with Cytoscape software (Figure 6D).

Finally, we plotted ROC curves of the 7 m1A RGs in the Tumor and Normal groups of the TCGA-KIRC dataset and presented the results (Figure 6E-K). The expression of ALBKH1 (AUC = 0.610), ALBKH3 (AUC = 0.627), TRMT61B (AUC = 0.605), YTHDC1 (AUC = 0.518), and YTHDF2 (AUC = 0.592) in the TCGA-KIRC dataset in the Tumor and Normal groups had low accuracy. In contrast, TRMT6 (AUC = 0.857) and TRMT10C (AUC = 0.808) showed relatively high accuracy in distinguishing between Tumor and Normal groups in the TCGA-KIRC dataset.

Prognostic analysis of m1A-related genes

We performed a prognostic analysis of 7 m1A RGs, including ALKBH1, ALKBH3, TRMT10C, TRMT61B, TRMT6, YTHDC1, and YTHDF2, using the constructed Lasso regression model, and considered associated molecules to be statistically significant at P-value <0.05. We plotted prognostic survival KM curves for each of the 7 m1A RGs using P-value <0.05 as the screening criterion, and obtained 5 eligible m1A RGs (ALBKH1, TRMT10C, TRMT61B, YTHDC1, YTHDF2) (Figure 7A-E). The results showed that ALBKH1 (P-value < 0.001), TRMT10C (P-value = 0.008), TRMT61B (P-value < 0.001), YTHDC1 (P-value = 0.002), and YTHDF2 (P-value < 0.031) were all statistically significant.

Subsequently, we also plotted the ROC curves of 5 m1A RGs based on disease-specific survival (DSS) events in the TCGA-KIRC dataset, and presented the results with AUC>0.6 (Figure 7F-I, Table 6). As shown in Figure 5, the expression of m1A RGs ALBKH1 (AUC = 0.641, Figure 7F), YTHDC1 (AUC = 0.684, Figure 7H), and YTHDF2 (AUC = 0.605, Figure 7I) only had low accuracy in predicting survival events based on clinical variables in the TCGA-KIRC dataset. TRMT61B (AUC = 0.803, Figure 7G) had moderate accuracy in predicting survival events based on clinical variables in the TCGA-KIRC dataset.

Table 6: Patient Characteristics of KIRC patients in the TCGA datasets. 

Construction of a prognostic model for DEGs in the high- and low-risk groups of the TCGA-KIRC dataset

We conducted a univariate Cox regression analysis to evaluate the prognostic relationship between each variable and m1A RGs. Next, we included all variables that were significantly associated with m1A RGs in the univariate analysis, and built a prognostic model to distinguish high-risk and low-risk groups using multivariate Cox regression analysis. Additionally, we used a forest plot to display the results of the univariate Cox regression analysis, which helped readers better understand the prognostic relationship between each variable and m1A RGs (Figure 8A,Table 7).

Subsequently, we utilized nomogram analysis to assess the prognostic ability of the risk model and constructed a nomogram (Figure 8B) to facilitate visualization. Based on the chart, we observed that the ALKBH1 and TRMT6B1 genes exhibit higher utility compared to other variables. Additionally, we performed a prognostic calibration analysis for the 5 m1A RGs using the nomogram at 1-year, 3-year, and 5-year (Figure 8C-E) intervals and constructed calibration curves (Figures 8F-H). The x-axis of the calibration curve represents the forecasted probability of survival, while the y-axis corresponds to the actual survival probability. The degree of proximity between the calibration curve and the diagonal line serves as an indicator of the accuracy of the model's predictions.

Finally, we utilized DCA to assess the clinical utility of the prognostic model at 1year, 3years, and 5 years (Figure 8F-H). The threshold probability is depicted on the x-axis of the DCA, while the y-axis represents the net benefit. The stable range of x values above the lines for "all positive" and "all negative" could be used to judge the results, with a more extensive x-value range indicating better model performance. The results show that our constructed multivariate Cox regression model has a clinical predictive effect of 5 years > 3 years > 1 year.

Table 7: COX regression to hub genes associated with OS in TCGA-KIRC.

Immunohistochemical analysis of prognostic m1A RGs

We conducted IHC analysis using gene-specific antibodies against prognostic 5 m1A RGs in normal kidney and KIRC tissues, as provided by the HPA database. The IHC analysis revealed more intense ALKBH1 and TRMT61B staining in normal kidney tissue (Figure 9A, 9E) than in KIRC tumor tissue (Figure 9B, 9F), indicating lower expression levels in KIRC tumor tissue compared to normal kidney tissue. Conversely, YTHDC1 and YTHDF2 exhibited weaker staining in normal kidney tissue (Figure 9G, 9I) and more vital staining in KIRC tumor tissue (Figure 9H, 9J), suggesting significantly higher expression levels of these genes in KIRC tumor tissue than in normal kidney tissue. TRMT10C did not exhibit significant differences in staining intensity between normal kidney tissue (Figure 9C) and KIRC tumor tissue (Figure 9D), indicating that its expression levels did not significantly change in KIRC tumor tissue compared to normal kidney tissue.

Discussion

KIRC patients exhibit a higher incidence of tumor recurrence and metastasis compared to other subtypes of kidney cancer. In clinical practice, most KIRC patients experience significant side effects after undergoing chemotherapy, radiation therapy, or targeted therapy, and surgical resection remains the primary treatment method [17]. Despite recent advancements in diverse treatment approaches, the prognosis for most patients with advanced KIRC remains dismal. Therefore, accurate prognostication of KIRC patients is crucial for selecting treatment options and directly affects treatment outcomes. We urgently need novel biomarkers to enhance the accuracy of early screening and enable more precise monitoring of tumor progression. With the advancement of genomics and molecular biology, numerous studies have shown that m1A RGs play a pivotal role in the occurrence and progression of various tumors. For instance, in pancreatic cancer, m1A-mediated regulation of gene expression could serve as a prognostic indicator. ALKBH1 exerts a crucial promoting effect on the progression of pancreatic cancer by activating the mTOR and ERBB signaling pathways. In gastrointestinal tumors, studies have shown that m1A modification promotes tumor progression by regulating the PI3K/AKT/mTOR and ErbB signaling pathways [18]. In summary, m1A RGs hold tremendous potential in tumor diagnostic, therapeutic, and prognosis evaluation. However, the extent of our comprehension regarding the regulatory function of m1A RGs in KIRC remains exceedingly restricted.

Our study initially performed Lasso regression analysis to identify 7 m1A RGs (ALKBH1, ALKBH3, TRMT10C, TRMT6, TRMT61B, YTHDC1, and YTHDF2) for constructing a prognostic risk model. Subsequently, we validated the expression differences of these 7 m1A RGs in the TCGA-KIRC dataset and the GSE53757 and GSE26574 datasets. By evaluating the accuracy of the m1A RGs in kidney cancer diagnosis using ROC curves, significant accuracy in kidney cancer diagnosis was found in the expression of TRMT6 (AUC = 0.857) and TRMT10C (AUC = 0.808). This discovery could offer new clues for early identification and treatment of KIRC. Subsequently, we conducted a comprehensive analysis of the possible biological roles and signaling pathways within these m1A RGs through GO, GSEA, and GSVA analyses and constructed a PPI network. After performing Kaplan-Meier survival analysis on the 7 m1A RGs, We selected 5 m1A RGs (ALKBH1, TRMT10C, TRMT61B, YTHDC1, and YTHDF2) to construct prognostic risk models. We then plotted ROC curves based on DSS events and found that TRMT61B (AUC = 0.803) might have clinical value in predicting survival events in KIRC patients. We also constructed a nomogram to improve the accuracy of 1-year, 3-year, and 5-year OS predictions and plotted calibration curves. Our results showed that ALKBH1 and TRMT6B1 genes exhibit higher utility compared to other variables. Subsequently, we utilized DCA to evaluate the clinical utility of our prognostic model, revealing that the clinical predictive performance of the model was superior at 5-year. Finally, we conducted immunohistochemistry analysis based on gene-antibody staining and found that ALKBH1 and TRMT61B expression levels were lower in KIRC tissues. Our findings suggest that the risk model is capable of effective predicting the prognosis of KIRC patients. ALKBH1 and TRMT61B might play a critical role in the prognosis of KIRC, and their aberrant expression could potentially impact essential biological processes. Therefore, future research on these genes carries paramount clinical significance.

ALKBH1 is a tRNA demethylase that mediates the demethylation of m1A in tRNA. This demethylation leads to weakened translation initiation and reduced involvement of tRNA in protein synthesis [11]. This unique mechanism might provide a new perspective for studying the progression of various cancers. For example, studies have found that low expression of ALKBH1 is correlated with an unfavorable prognosis in individuals with pancreatic cancer [16]. However, in most tumors, ALKBH1 acts as an oncogene. For example, in gastric cancer, overexpression of ALKBH1 significantly enhances tumors' survival and migration ability [19]. In colorectal cancer, elevated ALKBH1 expression has been linked to increased cell migration and invasion, as well as poor clinical outcomes [20]. TRMT61B is an mitochondrial RNA methyltransferase that could manipulate m1A levels to interfere with translation and catalyze the formation of m1A modifications [21]. Studies have revealed that TRMT61B is a novel susceptibility locus that might be associated with estrogen receptor-negative breast cancer, and as an overexpressed transcript contributes to the development of gastric cancer and might predict unfavorable clinical outcomes in patients with gastric cancer [19, 22]. The correlation analysis between the two in KIRC has not been specifically studied. Despite recent findings that suggest high expression of ALKBH1 and TRMT61B in KIRC patients might predict better survival rates, the current study lacks specificity and does not adequately explore the systemic relationship between KIRC and m1A RGs [23]. In summary, the expression of ALKBH1 and TRMT61B is related to the development and prognosis of various cancers. Therefore, further expansion of the correlation analysis of m1A RGs with KIRC, especially the in-depth study of the two regulatory genes of m1A, ALKBH1 and TRMT61B, would facilitate a more profound comprehension of their contributions to the progression and prognosis of KIRC, enabling accurate clinical outcome predictions and personalized treatment strategy development.

Utilizing GO, GSEA, and GSVA analysis, this study demonstrates the enrichment of m1A RGs in biological processes and pathways, revealing the involvement of multiple biological processes and associated pathways in the progression of KIRC. The enriched biological processes and pathways of m1A RGs in this study include mRNA methylation, RNA modification, tRNA modification, and mRNA modification, which have been reported to influence the pathogenesis of various diseases [24]. Furthermore, the functional relevance of m1A RGs enrichment in the cytoplasmic ribosome and mitochondrial matrix is noteworthy. Ribosome assumes a pivotal role in protein synthesis. It exerts a profound influence on cancer development by modulating mRNA translation. Mitochondria serve as the principal generators of intracellular energy. Damaged mitochondrial function leads to metabolic changes in tumor cells, promoting cell proliferation and metastasis [25, 26]. In terms of molecular function, significantly differentially expressed m1A RGs are significantly enriched in the catalytic activity of tRNA and RNA and the activity of methyltransferases in tRNA, mRNA, and RNA. GSEA indicated that m1A RGs might regulate the occurrence and development of KIRC by participating in pathways related to Interleukin Cytokines signaling, such as IL-5, IL-6, IL-7, and IL-27. Interleukins are a class of cytokines that regulate the development, proliferation, differentiation, and function of immune cells by triggering a series of signaling pathways through binding to surface receptors on immune cells [27]. Meanwhile, the P130Cas and Grb2-Sos linkage to MAPK signaling for integrins is also significantly enriched in m1A RGs. The interaction between the MAPK signaling pathway and integrin could regulate biological processes such as cell adhesion, movement, and proliferation. This signaling mechanism plays an important role in the invasion and metastasis of tumor cells, which might affect the occurrence and progression of KIRC [28, 29]. GSVA analysis revealed significant differences in the 32 hallmark gene sets enrichment between both groups in the TCGA-KIRC dataset. In summary, m1A RGs could affect the survival and growth of KIRC tumor cells by regulating various biological processes.

The primary innovative aspect of our study lies in the systematic analysis of patient cohorts within the TCGA database, which enabled us to establish an effective prognostic model based on 5 m1A RGs with potential clinical significance. According to our knowledge, the correlation analysis between m1A RGs and KIRC has not been previously investigated, highlighting the significance of our study. Additionally, our prognostic model demonstrates high precision in predicting survival, which might assist clinicians in making more accurate individualized survival predictions.

Undoubtedly, it is essential to consider limitations in our work. Firstly, the sample size of the external validation set is relatively small. Secondly, the TCGA database was deficient in essential clinical features of KIRC patients, such as their history of smoking and alcohol consumption, family medical history, and lifestyle, among other factors. Incorporating these variables could enhance the model's effectiveness and reliability. Finally, although the model's AUC demonstrates acceptable discriminatory ability, its performance still requires further improvement. Finally, the bioinformatics analysis proposed in this article should undergo validation through experimental studies.

Conclusions

In summary, our study utilizes the TCGA database through high-throughput bioinformatics analysis to assess the prognostic value of m1A RGs in patients with KIRC. We constructed a prognostic risk model with 5 m1A RGs that could predict the survival time of KIRC patients independently and effectively. Additionally, we developed a customized nomogram for predicting survival that might be particularly advantageous in assisting with clinical decision-making. In summary, our study might assist in identifying new prognostic biomarkers and potential therapeutic targets for KIRC patients.

Materials And Methods

Data Acquisition

We downloaded the KIRC dataset from TCGA using the “TCGAbiolinks” package (https://portal.gdc.cancer.gov/) [30]. By removing samples lacking clinical information, we obtained count sequencing data for 611 KIRC samples, including 72 normal samples and 539 tumor samples. The count sequencing data for KIRC was subsequently transformed into Fragments Per Kilobase per Million formats, and we obtained the corresponding clinical data from the “UCSC Xena” database (http://genome.ucsc.edu).

Furthermore, we retrieved expression profiling datasets for KIRC patients from the GEO database (https://www.ncbi.nlm.nih.gov/geo) [31] by utilizing the “GEOquery R” package [32]. Specifically, we obtained the GSE537574 [33] and GSE265745 [34] datasets (GPL570 and GPL11433 platforms, respectively), both annotated with the Affymetrix Human Genome U133 Plus 2.0 Array (HG U133 Plus 2.0) for probe name identification. The GSE53757 dataset contained 144 samples, with 72 KIRC patient tumor samples and 72 corresponding normal kidney samples. The GSE26574 dataset consisted of 16 samples, including 8 KIRC patient tumor samples and 8 corresponding normal kidney samples (Table 1). These datasets were utilized as validation sets for subsequent analyses.

Prognostic model for m1A-related differentially expressed genes

To construct a predictive model for m1A-related differentially expressed genes in KIRC, we utilized Least Absolute Shrinkage and Selection Operator (Lasso) regression analysis with the TCGA-KIRC dataset, a 10-fold cross-validation approach, and a significance threshold of P-value < 0.05. We ran the Lasso regression for 1,000 iterations to prevent overfitting. Lasso regression is a commonly used method for developing prognostic models that build on linear regression by adding a penalty term (lambda x the absolute value of the slope) to reduce overfitting and improve model generalization. We visualized the results of the Lasso regression and used a risk factor plot to explain the molecular expression patterns of m1A-related prognostic genes in each risk score group. Additionally, we examined the survival outcomes for each sample in the Lasso regression prognostic model to further analyze the possible effect of these genes on patient prognosis.

Risk score stratification in the model

To detect differentially expressed genes (DEGs) associated with high- and low-risk scores in the Lasso regression prognostic model, we utilized the “limma” package [35] to conduct differential expression analysis on expression profile data from the TCGA-KIRC dataset. We selected DEGs in KIRC patients by comparing the high- and low-risk score groups with a threshold of |logFC| > 1 and P-value < 0.05.

Gene Ontology (GO) Analysis

GO analysis is a highly utilized approach for conducting comprehensive functional enrichment analyses. GO analysis is designed to annotate the functional information of genes and their products and encompasses three domains: molecular functions (MF), biological processes (BP), and cellular components (CC) [36]. In this study, we conducted GO analysis on m1A RGs obtained from the TCGA-KIRC dataset using the “clusterProfiler R” [36]package. To ascertain the statistical significance of the results, we set the threshold at P-value < 0.05 and false discovery rate (FDR) q-value < 0.05 and performed P-value correction via the Benjamini-Hochberg (BH) method.

Gene Set Enrichment Analysis (GSEA) and Gene Set Variation Analysis (GSVA) between high- and low-risk Subgroups

GSEA is utilized to detect gene distribution patterns within a pre-defined gene set in a gene expression matrix sorted based on their correlation to a particular phenotype to ascertain their potential contribution towards the phenotype [37]. For this study, we first divided all the genes within the m1A RGs prognostic model of the TCGA-KIRC dataset into two distinct groups based on their logFC values, arranged in descending order. Next, we utilized the “clusterProfiler” package to conduct enrichment analysis of DEGs within the different groups of m1A RGs expression prognostic models using the TCGA-KIRC dataset. We set the following parameters for GSEA: the seed was 2022, 1000 permutations were performed, and gene sets included 10 genes and a maximum of 500 genes. The P-value correction method used was BH. The gene set "c2.all.v7.2.symbols.gmt" from the Molecular Signatures Database 3.0 (MSigDB) [38] was employed as the reference set, and significant enrichment was determined based on a P-value < 0.05 and q-value < 0.05.

GSVA is an unsupervised analysis method that determines the enrichment of various pathways across samples by converting gene expression matrices into gene set expression matrices to assess the gene set enrichment of the transcriptome [39]. We performed GSVA using the “h.all.v7.4.symbols.gmt” reference gene set from the MSigDB database (https://www.gsea-msigdb.org/) on the gene expression matrices of the different groups of the TCGA-KIRC dataset. The enriched pathways between two subgroups of KIRC were compared using functional differences calculated from the GSVA results, with statistical significance defined as P-value < 0.05.

Protein-protein Interaction Network (Ppi Network)

PPI network is established through interactions between multiple proteins, allowing for their involvement in various biological processes within organisms, such as signal transduction, energy metabolism, and gene regulation. The STRING database (http://string.embl.de/) [40] is commonly employed for searching identified proteins and predicting their interactions. In the current study, we constructed a PPI network for 7 m1A RGs through the STRING database and rendered it via Cytoscape. The highly interconnected nodes within the PPI network could potentially signify localized regions where distinct biological processes of molecular complexes are executed with specificity.

Prognostic Clinical Correlation Analysis

Initially, we performed a univariate Cox regression analysis on key genes to develop prognostic models for the DEGs in the high- and low-risk groups of the TCGA-KIRC dataset. Next, we included all variables for multivariate Cox regression analysis. Using the “rms” R package, we built column line plots, or nomogram plots, based on the multivariate Cox regression analysis results. Nomogram plots are graphical representations of multivariable regression models that predict the probability of a specific event by assigning values of the different independent variables to a numerical scale [41].

To evaluate the precision and discriminatory power of the nomogram plots, we employed calibration curves [42], which compare actual and predicted probabilities for different scenarios. Meanwhile, we used decision curve analysis (DCA) [43] to assess the accuracy and clinical usefulness of the DEG prognostic models, which compares the clinical outcomes of using the models with those of treating all or no patients. The R package “ggDCA” was used to generate DCA plots, which depict the net benefits of different prediction models across a range of decision threshold probabilities. We evaluated the models' net benefit across various decision threshold probabilities by evaluating their effectiveness in clinical settings.

Declarations

Ethics approval and consent to participate 

The public database mentioned in this study is publicly available for re-analyzing, and no ethical approval was required by the local ethics committees, so that this study does not require the ethics approval.

Consent for publication 

Not applicable.

Ethical Guidelines 

Not applicable.

Availability of data and materials 

All sequencing and clinical data of KIRC patients in this study were obtained from TCGA (https://portal.gdc.cancer.gov/) and GEO (https://www.ncbi.nlm.nih.gov/geo) databases. After organizing all the raw data, they were uploaded to a public database with the DOI: 10.6084/m9.figshare.22262221.

Funding Statement

This work was supported by Key Projects of University of South China (grant number：USC [2019] No. 02）

Competing interests 

The authors declare that they have no competing interests.

Author contributions

ZTF was responsible for data collection and initial analysis, YN designed the project and reviewed the manuscript, LX, JT and DBT were responsible for guiding the statistical methods and revising the discussion section. All authors read and approved the final manuscript.

Acknowledgements

The authors would like to express their gratitude to the staff members of public databases, such as TCGA, GEO, and HPA.

References

1. Bai X, Yi M, Dong B, Zheng X, Wu K: The global, regional, and national burden of kidney cancer and attributable risk factor analysis from 1990 to 2017. Exp Hematol Oncol 2020, 9:27.
2. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A: Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2018, 68(6):394-424.
3. Xia C, Dong X, Li H, Cao M, Sun D, He S, Yang F, Yan X, Zhang S, Li N et al: Cancer statistics in China and United States, 2022: profiles, trends, and determinants. Chin Med J (Engl) 2022, 135(5):584-590.
4. Moch H, Cubilla AL, Humphrey PA, Reuter VE, Ulbright TM: The 2016 WHO Classification of Tumours of the Urinary System and Male Genital Organs-Part A: Renal, Penile, and Testicular Tumours. Eur Urol 2016, 70(1):93-105.
5. Jonasch E, Walker CL, Rathmell WK: Clear cell renal cell carcinoma ontogeny and mechanisms of lethality. Nat Rev Nephrol 2021, 17(4):245-261.
6. Boccaletto P, Stefaniak F, Ray A, Cappannini A, Mukherjee S, Purta E, Kurkowska M, Shirvanizadeh N, Destefanis E, Groza P et al: MODOMICS: a database of RNA modification pathways. 2021 update. Nucleic Acids Res 2022, 50(D1):D231-D235.
7. Roundtree IA, Evans ME, Pan T, He C: Dynamic RNA Modifications in Gene Expression Regulation. Cell 2017, 169(7):1187-1200.
8. Dominissini D, Nachtergaele S, Moshitch-Moshkovitz S, Peer E, Kol N, Ben-Haim MS, Dai Q, Di Segni A, Salmon-Divon M, Clark WC et al: The dynamic N(1)-methyladenosine methylome in eukaryotic messenger RNA. Nature 2016, 530(7591):441-446.
9. RajBhandary UL, Stuart A, Faulkner RD, Chang SH, Khorana HG: Nucleotide sequence studies on yeast phenylalanine sRNA. Cold Spring Harb Symp Quant Biol 1966, 31:425-434.
10. Dai X, Wang T, Gonzalez G, Wang Y: Identification of YTH Domain-Containing Proteins as the Readers for N1-Methyladenosine in RNA. Anal Chem 2018, 90(11):6380-6384.
11. Liu F, Clark W, Luo G, Wang X, Fu Y, Wei J, Wang X, Hao Z, Dai Q, Zheng G et al: ALKBH1-Mediated tRNA Demethylation Regulates Translation. Cell 2016, 167(3):816-828 e816.
12. Safra M, Sas-Chen A, Nir R, Winkler R, Nachshon A, Bar-Yaacov D, Erlacher M, Rossmanith W, Stern-Ginossar N, Schwartz S: The m1A landscape on cytosolic and mitochondrial mRNA at single-base resolution. Nature 2017, 551(7679):251-255.
13. Vilardo E, Nachbagauer C, Buzet A, Taschner A, Holzmann J, Rossmanith W: A subcomplex of human mitochondrial RNase P is a bifunctional methyltransferase - extensive moonlighting in mitochondrial tRNA biogenesis. Nucleic Acids Res 2018, 46(20):11126-11127.
14. Wang Q, Zhang Q, Huang Y, Zhang J: m(1)A Regulator TRMT10C Predicts Poorer Survival and Contributes to Malignant Behavior in Gynecological Cancers. DNA Cell Biol 2020, 39(10):1767-1778.
15. Wang Y, Huang Q, Deng T, Li BH, Ren XQ: Clinical Significance of TRMT6 in Hepatocellular Carcinoma: A Bioinformatics-Based Study. Med Sci Monit 2019, 25:3894-3901.
16. Zheng Q, Yu X, Zhang Q, He Y, Guo W: Genetic characteristics and prognostic implications of m1A regulators in pancreatic cancer. Biosci Rep 2021, 41(4).
17. Motzer RJ, Jonasch E, Agarwal N, Alva A, Baine M, Beckermann K, Carlo MI, Choueiri TK, Costello BA, Derweesh IH et al: Kidney Cancer, Version 3.2022, NCCN Clinical Practice Guidelines in Oncology. J Natl Compr Canc Netw 2022, 20(1):71-90.
18. Zhao Y, Zhao Q, Kaboli PJ, Shen J, Li M, Wu X, Yin J, Zhang H, Wu Y, Lin L et al: m1A Regulated Genes Modulate PI3K/AKT/mTOR and ErbB Pathways in Gastrointestinal Cancer. Transl Oncol 2019, 12(10):1323-1333.
19. Li J, Zuo Z, Lai S, Zheng Z, Liu B, Wei Y, Han T: Differential analysis of RNA methylation regulators in gastric cancer based on TCGA data set and construction of a prognostic model. J Gastrointest Oncol 2021, 12(4):1384-1397.
20. Chen W, Wang H, Mi S, Shao L, Xu Z, Xue M: ALKBH1-mediated m(1) A demethylation of METTL3 mRNA promotes the metastasis of colorectal cancer by downregulating SMAD7 expression. Mol Oncol 2023, 17(2):344-364.
21. Chujo T, Suzuki T: Trmt61B is a methyltransferase responsible for 1-methyladenosine at position 58 of human mitochondrial tRNAs. RNA 2012, 18(12):2269-2276.
22. Kim J, Kwon J, Kim M, Do J, Lee D, Han H: Low-dielectric-constant polyimide aerogel composite films with low water uptake. Polymer Journal 2016, 48(7):829-834.
23. Song R, J-L Wong J: Altered expression of m1A regulatory genes is associated with oncogenic pathways, overall survival, and infiltration of immune cells in diverse human cancers. Genes & Diseases 2022.
24. Barbieri I, Kouzarides T: Role of RNA modifications in cancer. Nat Rev Cancer 2020, 20(6):303-322.
25. Goard CA, Schimmer AD: Mitochondrial matrix proteases as novel therapeutic targets in malignancy. Oncogene 2014, 33(21):2690-2699.
26. Goudarzi KM, Lindstrom MS: Role of ribosomal protein mutations in tumor development (Review). Int J Oncol 2016, 48(4):1313-1324.
27. Briukhovetska D, Dorr J, Endres S, Libby P, Dinarello CA, Kobold S: Interleukins in cancer: from biology to therapy. Nat Rev Cancer 2021, 21(8):481-499.
28. Cabodi S, del Pilar Camacho-Leal M, Di Stefano P, Defilippi P: Integrin signalling adaptors: not only figurants in the cancer story. Nat Rev Cancer 2010, 10(12):858-870.
29. Hamidi H, Ivaska J: Every step of the way: integrins in cancer progression and metastasis. Nat Rev Cancer 2018, 18(9):533-548.
30. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I et al: TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res 2016, 44(8):e71.
31. Barrett T, Troup DB, Wilhite SE, Ledoux P, Rudnev D, Evangelista C, Kim IF, Soboleva A, Tomashevsky M, Edgar R: NCBI GEO: mining tens of millions of expression profiles--database and tools update. Nucleic Acids Res 2007, 35(Database issue):D760-765.
32. Davis S, Meltzer PS: GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics 2007, 23(14):1846-1847.
33. von Roemeling CA, Radisky DC, Marlow LA, Cooper SJ, Grebe SK, Anastasiadis PZ, Tun HW, Copland JA: Neuronal pentraxin 2 supports clear cell renal cell carcinoma by activating the AMPA-selective glutamate receptor-4. Cancer Res 2014, 74(17):4796-4810.
34. Ooi A, Wong JC, Petillo D, Roossien D, Perrier-Trudova V, Whitten D, Min BW, Tan MH, Zhang Z, Yang XJ et al: An antioxidant response phenotype shared between hereditary and sporadic type 2 papillary renal cell carcinoma. Cancer Cell 2011, 20(4):511-523.
35. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK: limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research 2015, 43(7):e47-e47.
36. Yu G, Wang LG, Han Y, He QY: clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012, 16(5):284-287.
37. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES et al: Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005, 102(43):15545-15550.
38. Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P: The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst 2015, 1(6):417-425.
39. Wang HQ, Yu XD, Liu ZH, Cheng X, Samartzis D, Jia LT, Wu SX, Huang J, Chen J, Luo ZJ: Deregulated miR-155 promotes Fas-mediated apoptosis in human intervertebral disc degeneration by targeting FADD and caspase-3. J Pathol 2011, 225(2):232-242.
40. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P et al: STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res 2019, 47(D1):D607-D613.
41. Park SY: Nomogram: An analogue tool to deliver digital knowledge. J Thorac Cardiovasc Surg 2018, 155(4):1793.
42. Perkins NJ, Weck J, Mumford SL, Sjaarda LA, Mitchell EM, Pollack AZ, Schisterman EF: Combining Biomarker Calibration Data to Reduce Measurement Error. Epidemiology 2019, 30 Suppl 2(Suppl 2):S3-S9.
43. Van Calster B, Wynants L, Verbeek JFM, Verbakel JY, Christodoulou E, Vickers AJ, Roobol MJ, Steyerberg EW: Reporting and Interpreting Decision Curve Analysis: A Guide for Investigators. Eur Urol 2018, 74(6):796-804.
44. Magaki S, Hojat SA, Wei B, So A, Yong WH: An Introduction to the Performance of Immunohistochemistry. Methods Mol Biol 2019, 1897:289-298.
45. Thul PJ, Lindskog C: The human protein atlas: A spatial map of the human proteome. Protein Sci 2018, 27(1):233-244.
46. Engebretsen S, Bohlin J: Statistical predictions with glmnet. Clin Epigenetics 2019, 11(1):123.