Data on relevant clinical information, nucleotide variations and mRNA expressions were obtained from the TCGA public database. The dataset GSE29609 (20) was sourced from the Gene Expression Omnibus (GEO) database (ncbi.nlm.nih.gov/). Leveraging the capabilities of the ‘GEOquery’ package (bioconductor.org/packages/GEOquery/,2.72.0), both gene expression and clinical datasets were extracted. After integrating the clinical data with the GEO transcriptome data by sample name, the final sample size consisted of 39 cases. A total of 39 samples derived from GSE29609 were employed as a validation cohort to assess the precision of the prognostic model. Additionally, 100 TMRGs were identified using the Molecular Signatures Database (MSigDB; https://www.gsea-msigdb.org/gsea/msigdb/human/genesets.jsp; Table SI). First, on the MSigDB website, ‘tRNA’ was searched, and ‘GOBP_TRNA_MODIFICATION’ was selected. Next, the gene set tab-separated values metadata was downloaded. Finally, a comparative analysis of the expression levels of tRNA-modified proteins within ccRCC tissues was performed utilizing the comprehensive resources obtained from the Human Protein Atlas (HPA) database (proteinatlas.org/).

Univariate Cox regression analysis was utilized to identify tRNA modification genes associated with prognostic outcomes. To investigate the correlation among TMRGs, a protein-protein interaction (PPI) network was created utilizing the capabilities of the STRING database (string-db.org/) as follows: i) Network type, full STRING network active; ii) interaction source: Text mining, experiments, databases, co-expression, neighborhood, gene fusion, co-occurrence; and iii) minimum required interaction score, medium confidence 0.400. Subsequently, the Cytoscape Cytohubba (cytoscape.org/,v3.10.2) plug-in identified 26 hub genes and modules in the PPI network related to tRNA modification genes. Cluster analysis was performed using the package ‘ConsensusClusterPlus’ (bioconductor.org/packages/ConsensusClusterPlus/1.66.0) program to identify molecular subtypes associated with tRNA modification; parameters included maxK, 9; reps, 10; and pItem, 0.8. Kaplan-Meier analysis was used to analyze the differences between the two groups. χ² analysis was conducted concurrently to generate heat maps to systematically illustrate the associations between distinct clusters and their corresponding clinical features.

The least absolute shrinkage and selection operator (LASSO) and multivariable Cox regression analyses were performed to identify genes associated with ccRCC and establish prognostic signatures. The GraphPad software (version 9.5; Dotmatics) was employed to display the coefficients in the selected genes. The risk score was calculated as follows: ∑ⁿ_n coef(i) × Expr(i) (where i represents genes). Kaplan-Meier analysis and receiver operating characteristic (ROC) curves were generated to assess the predictive value of each attribute. Cox regression analyses were applied to verify whether a signature was an independent risk factor. Correlation analysis, stratification analysis and the establishment of a nomogram, based on risk scores and pertinent clinical attributes, were performed according to clinicopathological standards. Calibration plots were created for the 1-, 3- and 5-year survival rates to evaluate the alignment between the predicted probabilities and the actual survival results.

Gene Set Enrichment Analysis (GSEA; gsea-msigdb.org/gsea/index.jsp) was performed to examine pathway enrichment in high-risk groups. The internal reference gene sets included three categories: C2 KEGG, HALLMARK and C5GO. Significant results were indicated by normalized enrichment scores (NES) >1, a nominal P<0.05 and a false discovery rate (FDR) q-value <0.25.

A quartet of immune-related algorithms, including single sample (ss) GSEA, ESTIMATE, TIDE and CIBERSORT algorithms (RStudio IDE; Posit Software, PBC) were applied to assess and compare the immunological profiles of high-risk and low-risk cohorts. The ssGSEA methodology was utilized to elucidate the activity of immune cells, immune functions and the pathways associated with immunity for each sample. Marker genes for various immune cells were identified from previous studies, such as FOXP3/CTLA4, CD56, which are the marker genes of Tregs and Natural killer cells (21,22). The ESTIMATE method was utilized to calculate immune, stromal, estimate score and tumor purity by analyzing the ratio of immune cells to stromal cells. The composition of immune cell populations infiltrating each tumor specimen was predicted utilizing the CIBERSORT algorithm. Following cluster analysis and subsequent characterization, a comparative examination of the expression of the major histocompatibility complex (MHC) and immune checkpoint molecules was performed. Higher Tumor Immune Dysfunction and Exclusion (TIDE) scores are associated with longer survival and poorer response to checkpoint blockade therapy. Employing the TIDE database (tide.dfci.harvard.edu/), the TIDE scores for TCGA-kidney RCC (KIRC) cohort were calculated to predict immunotherapy responses in both subpopulations.

The R package ssGSEA (parameters included: Method, ‘ssgsea’; kcdf, ‘Gaussian’; abs.ranking, TRUE) was applied to assess the angiogenic activity, tumorigenic cytokine score, mesenchymal EMT and stemness score for individual samples within the TCGA-KIRC cohort.

Based on TCGA-KIRC somatic mutation data, the ‘maftools’ (bioconductor.org/packages/maftools/,2.18.0) package was applied to analyze gene mutations. Tumor mutational burden (TMB) was calculated for each patient and compared between the high-risk and low-risk cohorts. The TMB scores were subsequently utilized to perform a comprehensive survival analysis. The cBioPortal database (cbioportal.org/v4.0) indicated that some genes selected by the signature have somatic mutations.

The Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) was used to evaluate the impact of predictive signatures on treatment response in ccRCC. The ‘oncoPredict’ (https://bioconductor.org/packages/oncoPredict, 1.0.1) package was employed to retrieve the gene expression profiles from the GDSC dataset with the corresponding drug response data. Sensitivity scores were applied to predict the maximum IC₅₀ of all drugs in patients with KIRC.

The primer design process in the present study typically adhered to the following: i) The nucleotide sequence of the target gene was retrieved from gene databases, such as National Center for Biotechnology Information GenBank (ncbi.nlm.nih.gov/genbank/), ii) based on the sequence of the target gene, an oligonucleotide sequence was designed that can anneal and bind effectively; the optimal length for primers generally ranges from 18 to 25 nucleotides, which ensures high specificity and appropriate annealing temperatures, iii) calculation of primer melting temperature (Tm) was performed, aiming for a Tm between 55–65°C as recommended; primers were calculated using online tools such as OligoCalculator (idtdna.com/calc/analyzer) and iv a design software such as Primer3 (tmcalculator.neb.com; v2.6.1) minimized non-specific binding between primer sequences and target sequences. Accession numbers are shown in Table SII.

HK2 and ccRCC lines (786-O, 769-P, ACHN, Caki-1 and OS-RC-2) were purchased from the Shanghai Institutes for Biological Sciences. HK2 and ACHN cells were cultured in MEM (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (FBS; Gibco; Thermo Fisher Scientific, Inc.), while 786-O, 769-P and OS-RC-2 cells were cultured in RPMI-1640 (Gibco; Thermo Fisher Scientific, Inc.) medium supplemented with 10% FBS. Caki-1 cells were maintained in a McCoy's 5A medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% FBS. All cell lines were maintained in a sterile incubator (Thermo Fisher Scientific, Inc.) at 37°C with 5% CO₂. Total RNA was extracted from the tissues or cell lysates of RCC patients utilizing a TRIzol™ reagent kit (Qiagen, Inc.). Subsequent cDNA synthesis employed a Thermo-script RT kit (Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. The 2(−ΔΔ C(T)) Method to analysis the relative gene expression (23). qPCR was conducted utilizing the CFX96™ Real-Time System (Bio-Rad Laboratories, Inc.) with SYBR Green PCR reagent (Takara Bio, Inc.). The thermocycling conditions were as follows: 95°C for 5 min, followed by 40 cycles of 95°C for 15 sec, 60°C for 25 sec and 72°C for 30 sec. GADPH served as a standardized internal reference for normalization. The primer sequences utilized for qPCR are listed in Table SIII.

All statistical analyses were performed using R software (version 4.3.1; Posit Software, PBC) and GraphPad Prism (version 9.5; Dotmatics). The data are presented as mean ± standard deviation (SD) or mean ± standard error of the mean (SEM), as specified in the figure legends. P<0.05 was considered to indicate a statistically significant difference. For survival analysis, Kaplan-Meier curves were generated and the log-rank test was applied to evaluate differences between groups. Univariate and multivariate Cox proportional hazards regression models were used to determine independent prognostic factors. The LASSO regression and multivariate Cox regression were performed to construct the prognostic model. Comparisons between two groups were conducted using the Wilcoxon rank-sum test or Student's t-test, depending on the normality of the data distribution assessed by the Shapiro-Wilk test. For comparisons among multiple groups, one-way ANOVA followed by Tukey's post hoc test was applied when the data were normally distributed, while the Kruskal-Wallis test followed by Dunn's post hoc test was used for non-normally distributed data. χ² analysis was used to examine associations between categorical variables. GSEA was conducted using the clusterProfiler package (bioconductor.org/packages/clusterProfiler, 4.12.0). P<0.05 and FDR <0.25 were considered to indicate a statistically significant difference in terms of enrichment. The TIDE score differences between groups were assessed using the Mann-Whitney U test. The correlation between continuous variables was evaluated using Spearman's or Pearson's correlation analysis, depending on the distribution of the data. The ROC curve and area under the curve (AUC) were used to assess the predictive accuracy of the prognostic model. Differences in TMB scores were analyzed using the Wilcoxon test. Drug sensitivity prediction was performed using the oncoPredict package and IC₅₀ values were compared using the Wilcoxon test. All experiments were performed with at least three biological replicates. Statistical significance for each analysis is indicated in the corresponding figure legends.

A total of 100 genes (Table SI) associated with tRNA modification were identified in the MSigDB. Univariate Cox analysis identified 24 genes significantly associated with ccRCC prognosis, including six protective factors, leucine carboxyl methyltransferase 2 (LCMT2), SEPSECS, TRNA methyltransferase O (TRMO), TRNA methyltransferase 1L (TRMT1L), TRMT5, TRMT61B and 18 risk factors (Fig. 1A). The expression levels of the 24 genes related to prognosis in ccRCC tissues were compared with that of normal kidney tissues (Fig. 1B). The PPI network and correlation analysis illustrated the complex associations among the 24 prognosis-related genes (Fig. 1C and D). Clustering studies using these 24 genes (Table SIV) indicated that the optimal classification was to stratify patients with ccRCC into two categories (k=2) based on the expression patterns of 24 tRNA modification-associated genes (Fig. 2A-C). The cases represented by clusters 1 and 2 are shown in Table SV. Survival analysis of the two subgroups demonstrated that patients classified within cluster 1 exhibited a significantly more favorable prognosis compared with that of the patients in cluster 2 (Fig. 2D) The heat map significant associations between the clusters and clinicopathological parameters including clinical stage, T, M stage and grade, advanced stage (III/IV) or T4/N1/M1are enriched in Cluster2 (Fig. 2E).

In addition, cluster 2 exhibited significant association with decreased expression levels of various MHC molecules (Fig. 2F). Considering the disparities in immune infiltration between the two cohorts, the relationship with immune checkpoints was examined. Cluster 2 demonstrated significantly increased expression levels of immune checkpoints such as TNF receptor superfamily member 7, programmed cell death protein 1 (PD-1), T-cell immunoreceptor with Ig and ITIM domains (TIGIT), lymphocyte-activation gene 3 (LAG3), CD40 and cytotoxic t-lymphocyte antigen 4 (CTLA-4; Fig. 2G). CIBERSORT algorithm showed that cluster 2 exhibited a markedly increased infiltration of immune cells, encompassing CD8⁺ T lymphocytes, activated memory CD4⁺ T lymphocytes, follicular helper T cells and regulatory T cells (Tregs; Fig. 2H) Using the ESTIMATE methodology, cluster 2 demonstrated significantly increased tumor purity, accompanied by decreased stromal and ESTIMATE scores (Fig. 2I). Prolonged survival and diminished efficacy of checkpoint blockade therapy are linked to elevated TIDE scores (24). Cluster 2 exhibited significantly increased TIDE scores compared with that of cluster 1 (Fig. 2J).

LASSO analysis in conjunction with multivariate Cox regression analysis was applied to refine the gene set of the present model, which resulted in the inclusion of six genes, FtsJ RNA 2′-O-methyltransferase 1 (FTSJ1), LCMT2, METTL6, pseudouridine synthase 1 (PUS1), TRMO and TRMT5, within the signature (Fig. 3A). The coefficients of the six genes in the signature are shown Fig. 3B. The relationships between risk score and six genes are illustrated in Fig. 3C. The risk score was calculated according to the formula: ∑ⁿ_n coef(i) × Expr(i) and the model was then used to classify patients with ccRCC into low- and high-risk groups using optimized cut-off values (1.053). The risk score was correlated with clinicopathological parameters such as grade and clinical stage (Fig. 3D). At the 1-, 3- and 5-year AUCs (0.751, 0.726 and 0.739), patients that exhibited increased risk scores demonstrated significantly poorer prognosis compared with those with low-risk scores (Fig. 3E-G). The results of the survival analysis were validated through the utilization of the GSE29609 dataset (Fig. 3H-J). Additionally, changes in risk assessments among subgroups associated with numerous clinicopathological parameters were observed. Patients within grades 3 and 4, and stages T3-T4 and III–IV exhibited significantly increased risk scores, which suggested that more advanced tumors were associated with higher risk scores (Fig. 4A-C). Finally, Cox regression analyses demonstrated that the signature constituted an independent risk factor (Fig. 4D and E). From the Cox regression analysis results, age, score, clinical stage and signature were included when constructing the nomogram, with the gene signature as the key component of significance (Fig. 4F). According to the calibration plot, the survival times at 1-, 3- and 5-years exhibited alignment with the projected survival estimates (Fig. 4G).

GSEA was employed to examine pathways that govern tumorigenesis within the high-risk cohort. Numerous tumor-associated pathways were significantly enriched within the high-risk cohort (Fig. 5B). The findings demonstrated that the high-risk cohort exhibited a significant enrichment in the EMT pathway [NES=1.72; nominal (NOM) P=0.029; FDR q-value=0.12], IL6-JAK-STAT3 signaling pathway (NES=1.75; NOM P=0.027; FDR q-value=0.18), P53 pathway (NES=1.80; NOM P=0.018; FDR q-value=0.14) and the interactions of cytokine-cytokine receptors (NES=1.81; NOM P=0.014; FDR q-value=0.17; Fig. 5A). The scores of angiogenic activity, mesenchymal EMT, tumorigenic cytokines and stemness were assessed in patients with ccRCC. The scores for EMT were significantly decreased, while the assessments of tumorigenic cytokines, angiogenic activity and stemness exhibited a significant increase within the high-risk cohort compared with that of the low-risk cohort (Fig. 5C). The risk score demonstrated a positive correlation with the tumorigenic cytokines scores (R=0.22; P=8.4×10⁻⁰⁷), angiogenic activity scores (R=0.25; P=2.4×10⁻⁰⁸) and stemness scores (R=0.12; P=0.0092), whereas the risk score negatively correlated with mesenchymal EMT scores (R=−0.11; P=0.015; Fig. 5D).

Previous studies showed that the TME serves a pivotal role in tumorigenesis, and GSEA conducted in the present study indicated that several immune-related pathways were associated with the high-risk group. Therefore, the association with the tumor immune microenvironment was investigated. The ssGSEA algorithm demonstrated that the high-risk group showed significantly enhanced immune-related activities and pathways, as well as a more pronounced infiltration of immune cells in comparison with the low-risk group (Fig. 6A and D). Additionally, MHC expression levels were significantly increased in the high-risk group, compared with that of the low-risk group (Fig. 6B). The high-risk cohort exhibited significantly increased expression levels of immune checkpoint inhibitors, including TNFRSF9, PDCD1, TIGIT, LAG3, CD40 and CTLA-4 (Fig. 6C). The CIBERSORT algorithm indicated that plasma cells, CD8⁺ T cells, activated memory CD4⁺ T cells, Tregs, follicular helper T cells and M0 macrophages emerged as the predominant immune cells infiltrating the high-risk cohort (Fig. 6E and G). High TIDE scores in high-risk subtypes suggested that patients with high-risk ccRCC may exhibit a diminished response to immunotherapy (Fig. 6F). The ESTIMATE algorithm demonstrated that the high-risk cohort exhibited decreased tumor purity, and increased immunological and estimated scores compared with that of the low-risk group (Fig. 6H).

Primary nucleotide variation data pertaining to ccRCC were acquired from the TCGA database to investigate the genomic mutation disparities between high-risk and low-risk cohorts. Notably, the five genes exhibiting the highest mutation frequency within the low-risk group were identified as von Hippel-Lindau tumor suppressor (VHL; 47%), polybromo 1 (PBRM1; 43%), titin (TTN; 19%), mTOR (7%) and SET domain containing 2 histone lysine methyltransferase (SETD2; 7%; Fig. 7A). In the high-risk group, the top five genes with the highest mutation frequency were VHL (46%), PBRM1 (39%), SETD2 (20%), TTN (20%) and BRCA1-associated protein 1 (15%; Fig. 7B). Using a cut-off value is 0.925). The cases stratified into a low-TMB cohort exhibited a significantly prolonged survival duration compared with of the high-TMB cohort (Fig. 7C). The present model demonstrated that the high-risk and high-TMB cohort exhibited a markedly poorer prognosis when compared with that of the low-risk and low TMB group (Fig. 7D). Investigation into the mutation rates of marker genes demonstrated that FTSJ1 exhibited amplification mutations, while METTL6 was characterized by a predominance of profound deletion mutations. By contrast, LCMT2 displayed a higher frequency of missense mutations (Fig. 7E).

GDSC was used to predict the therapeutic response of high- and low-risk patient cohorts to widely utilized chemotherapeutic agents. Significant differences were observed in the responsiveness of the high- and low-risk cohorts to an array of chemotherapy agents, including tyrosine kinase (Fig. 8A), mTOR (Fig. 8B), AKT and Erk inhibitors (Fig. 8C). Additionally, the three-dimensional structures of potential drugs (sorafenib, axitinib, vistusertib, dactolisib, taselisib and afuresertib) were displayed using the PubChem database (Fig. 8D).

The protein expression levels of the six signature genes in renal cancer were derived from immunohistochemical staining data from the HPA database (Fig. 9A). To validate the results, total RNA was extracted from an array of ccRCC cell lines (786-O, 769-P, ACHN, Caki-1 and OSRC-2) and HK2 cells. Subsequently, the mRNA expression levels FTSJ1, LCMT2, METTL6, PUS1, TRMO and TRMT5 of mRNA were assessed. The RT-qPCR results showed that the mRNA expression levels of FTSJ1, METTL6 and PUS1 were markedly increased in ccRCC cells, while the mRNA expression levels of LCMT2, TRMO and TRMT5 were significantly decreased ccRCC cells compared with that of HK2 cells (Fig. 9B).