Identification and Validation of a Prognostic Model Based on Tumour Necrosis Factor‐Related mRNAs for Kidney Renal Clear Cell Carcinoma

TL;DR

This study identifies a new model using tumor necrosis factor-related mRNAs to predict survival and guide treatment in kidney cancer.

Contribution

A novel TRmRNA-based prognostic model for kidney renal clear cell carcinoma is developed and validated.

Findings

A nine TRmRNA signature was identified, with high-risk patients showing significantly shorter survival.

The model's nomogram achieved 76.6% accuracy in predicting patient outcomes.

High-risk patients exhibited higher tumor mutational burden and distinct enriched biological pathways.

Abstract

Tumour necrosis factor (TNF) plays a critical role in tumour progression, but the specific involvement of mRNA in this process, particularly in kidney renal clear cell carcinoma (KIRC) remains insufficiently understood. Our study aims to develop a TNF‐related mRNA (TRmRNA) model to predict prognosis and inform treatment strategies in KIRC. KIRC expression data from The Cancer Genome Atlas (TCGA) and TNF‐related genes (TRGs) from the Genecards database were used to construct and validate a TRmRNA prognostic model. A nomogram integrating clinical features with the risk model was also developed to enhance prognostic accuracy. Enrichment analysis, drug sensitivity analysis and RT‐qPCR validation were performed to further explore the biological mechanisms and clinical applicability of the model. A prognostic signature consisting of nine TRmRNAs was identified. Kaplan–Meier analysis showed that the high‐risk (HRK) group had significantly shorter overall survival (OS) compared to the low‐risk (LRG) group (p < 0.001). The nomogram, incorporating the risk model, yielded an area under the curve (AUC) of 0.766, indicating robust prognostic accuracy. Enrichment analysis identified solute sodium symporter and proximal tubule transport pathways enriched in the LRG group, whereas the HRK group exhibited enrichment in CD22‐mediated BCR regulation and immunoglobulin complex pathways. The HRK group also showed a higher tumour mutational burden (TMB), correlating with a poorer prognosis. RT‐qPCR confirmed the differential expression of mRNAs in KIRC cells. The TRmRNA‐based prognostic model holds significant promise for predicting patient outcomes and guiding personalised treatment strategies in KIRC.

Figures

Flow Chart.

PPI network (A); Filter for TNF‐related mRNAs. Volcano plot of differentially expressed TNF‐associated mRNAs (B); LASSO regression analysis (C, D).

The model prediction effect is validated by the train group, test group, and entire group. K‐M analysis (A–C) and time‐dependent ROC curves (D–F) to compare the survival of the high‐risk group and low‐risk group.

PCA and independent prognostic analysis of the signature. PCA based on all genes (A), all mRNAs (B), TNF‐related mRNAs (C), and risk signature (D); univariate (E) and multivariate (F) independent prognostic analysis.

Further validation of model effects. Survival curves of patients in different clinical states (A–F).

Nomogram predicts patient prognosis. Decision curve to test for forecast value (A) ROC curves containing different clinical information (B) A clinical prognosis nomogram is constructed by age, gender, risk, and stage together (C). Nomogram with (D) and without (E) risk model.

Tumour mutation burden in different risk groups. Percentage bar graph showing TMB for different risk subgroups (A) High‐risk group waterfall chart (B) Low‐risk group waterfall chart (C).

Analysis of tumour immune microenvironment. Violin plots of differences in ESTIMATE scores (A), stromal scores (B), immune scores (C) and tumour purity (D) for different risk subgroups; Bubble plots of correlations between immune cells and risk scores under six algorithms (E); Proportions of 22 immune cells in two subgroups under the CIBERSORT algorithm (F); single sample gene set enrichment analysis (G). *p < 0.05, p < 0.01, *p < 0.001.

Tables

TABLE 1. Clinical information of the patients in the testing and training groups.

Characteristics	Train cohort (n = 264)		Test cohort (n = 262)		Entire cohort (n = 526)
	n	%	n	%	n	%
Age
≤ 65	173	65.53	174	66.41	347	65.97
> 65	91	34.47	88	33.59	179	34.03
Status
Alive	178	67.42	177	67.56	355	67.49
Dead	86	32.58	85	32.44	171	32.51
Gender
Female	94	35.61	89	33.97	183	34.79
Male	170	64.39	173	66.03	343	65.21
Stage
Stage I	146	55.30	115	43.89	261	49.62
Stage II	22	8.33	35	13.36	57	10.84
Stage III	53	20.08	70	26.71	123	23.38
Stage IV	42	15.91	40	15.27	82	15.59
Unknown	1	0.38	2	0.76	3	0.57
T stage
T1	150	56.82	117	44.66	267	50.76
T2	28	10.61	41	15.65	69	13.12
T3	81	30.68	98	37.40	179	34.03
T4	5	1.89	6	2.29	11	2.09
M stage
M0	205	77.65	213	81.30	418	79.47
M1	58	21.97	38	14.50	96	18.25
Unknown	1	0.38	11	4.20	12	2.28
N stage
N0	102	38.64	136	51.91	238	45.25
N1	9	3.41	7	2.67	16	3.04
Unknown	153	57.95	119	45.42	272	51.71
Race
White	226	85.61	231	88.17	457	86.88
Black or African American	29	10.98	25	9.54	54	10.27
Asian	5	1.89	3	1.15	8	1.52
Unknown	4	1.52	3	1.15	7	1.33

TABLE 2. Clinical information for 526 patients in different risk categories.

Characteristics	High‐risk group (n = 263)		Low‐risk group (n = 263)
	n	%	n	%
Age
≤ 65	166	63.12	181	68.82
> 65	97	36.88	82	31.18
Status
Alive	132	50.95	223	84.79
Dead	131	49.05	40	15.21
Gender
Female	94	35.74	89	33.84
Male	169	64.26	174	66.16
Stage
Stage I	93	35.36	168	63.88
Stage II	33	12.55	24	9.13
Stage III	71	27.00	52	19.77
Stage IV	64	24.33	18	6.84
Unknown	2	0.76	1	0.38
T stage
T1	98	37.26	169	64.26
T2	42	15.97	27	10.27
T3	112	42.59	67	25.48
T4	11	4.18	0	0
M stage
M0	184	69.96	234	88.97
M1	60	22.81	18	6.84
Unknown	19	7.22	11	4.18
N stage
N0	129	49.05	109	41.44
N1	12	4.56	4	1.52
Unknown	122	46.39	150	57.03
Race
White	223	85.93	234	88.97
Black or African American	34	12.93	20	7.60
Asian	3	1.14	5	1.90
Unknown	3	1.14	4	1.52

Full text

Introduction

1

Kidney renal clear cell carcinoma (KIRC) represents one of the most prevalent urogenital cancers, with the incidence of this cancer rising by 2%–3% each year [1]. The tumour‐node‐metastasis (TNM) classification by the American Joint Committee on Cancer (AJCC) is a common method to evaluate disease progression for prognosis [2]. However, this staging system has notable limitations in accurately predicting patient outcomes, reducing its utility in optimising treatment strategies [3]. Consequently, the demand for accurate and dependable prognostic models to support clinical decisions and enhance KIRC patient management has become increasingly urgent. In recent years, biomarker‐based models have gained attention as a promising tool, offering improved accuracy and precision in cancer prognosis assessment.

The tumour necrosis factor (TNF) family comprises 19 ligands and 29 receptors, both of which are significant in cancer progression [4]. Previous research has demonstrated that signatures based on the TNF family can serve as effective prognostic tools and predict the response to immunotherapy and chemotherapy in various cancers, such as lung adenocarcinoma (LUAD) and colorectal cancer [5, 6]. Moreover, TNF and its receptors are emerging as promising therapeutic targets for treating rheumatic and other inflammatory diseases [7]. Therefore, the TNF family represents a promising avenue for enhancing cancer prognosis and developing novel immunotherapeutic strategies. At the same time, messenger RNA (mRNA) has become an essential component in cancer prognosis, with numerous mRNA‐based prognostic models successfully developed for different types of cancer. For example, Zou et al. established an mRNA‐based prognostic model for cervical squamous carcinoma and adenocarcinoma, which was validated for identifying precise and effective biomarkers [8]. Likewise, Jiang et al. proposed an mRNA‐based model for stomach adenocarcinoma (STAD) that provided reliable predictions of clinical outcomes and targeted therapy responses [9]. However, to date, no research has focused on developing a prognostic model for KIRC based on TNF‐related mRNA (TRmRNA).

Building upon this foundation, a predictive model using TRmRNAs was developed and validated to explore its underlying mechanisms and evaluate its clinical relevance. This was accomplished by performing a range of analyses, such as enrichment analysis, assessment of tumour mutational burden (TMB), profiling of the tumour microenvironment (TME) and predictions of drug sensitivity.

Materials and Methods

2

Data Acquisition and Organisation

2.1

Data for our research were primarily collected from the University of California, Santa Cruz (UCSC) Xena database (https://xena.ucsc.edu/, available until July 1st, 2024) [10]. We collected KIRC transcriptional data for 607 cases and clinical data for 979 cases from The Cancer Genome Atlas (TCGA) database via Xena. Our study included only patients who had complete clinical and transcriptional data available from datasets. In total, 598 samples were used for the further analysis. To reduce potential technical variability in sequencing depth and platform artefacts within TCGA data, we used FPKM‐normalised transcriptomic data downloaded via the UCSC Xena platform, which had undergone uniform preprocessing.

Select TNF‐Related Genes (TRGs) and TRmRNAs

2.2

TRGs were sourced from the GeneCards database, which were subsequently selected for further analysis [11]. To investigate the relationship between TRGs and all mRNAs, Pearson correlation analysis was conducted (|correlation coefficient| > 0.3 and p < 0.05), aiming to identify TRmRNAs. To further refine the results, variance analysis was conducted through the “DESeq2” package to filter mRNAs (p < 0.05 and |log2‐fold change| > 1.2) [12].

Subsequently, we split the KIRC patient data into train and test groups at a 1:1 ratio. TRmRNAs linked to survival were detected through univariate Cox (uCox) analysis. TRmRNAs strongly correlated with prognosis were pinpointed through the least absolute shrinkage and selection operator (LASSO) and multivariate Cox (mCox) analysis.

Establish TRmRNA Predictive Model

2.3

The findings were used to develop a predictive model, with the risk score determined by the following way: risk score = (coefficient of mRNA1 × expression level of mRNA1) + … + (coefficient of mRNAn × expression level of mRNAn). Using the midpoint of the score, patients were divided into two groups: high‐risk (HR) and low‐risk (LR) groups. Kaplan–Meier (K‐M) analysis for these cohorts was formulated through the “survival” package in R to analyse overall survival (OS) differences between the HR and LR groups. The “timeROC” package in R was employed to evaluate the model's predictive ability over 1, 3 and 5 years by generating Receiver operating characteristic (ROC) curves for the cohorts. Confidence intervals (95% CI) for AUC values were estimated via bootstrapping with 1000 iterations. The statistical significance of AUCs was assessed using permutation‐based p‐values.

Principal component analysis (PCA) was applied to examine the expression profiles of all genes, mRNAs, TRmRNAs and those included in the prognostic models, aiming to evaluate subgroup classification effectiveness. The uCox and mCox analyses were conducted to assess the independent prognostic capacity for our model, integrating risk scores with clinical variables. To further confirm the model's performance, K‐M survival analysis assessed its prognostic accuracy across multiple clinical subgroups.

Permutation Testing

2.4

To evaluate whether the predictive performance of the TRmRNA‐based prognostic model could have arisen by chance, we conducted permutation testing. We performed 1000 iterations, in each of which the overall survival status was randomly permuted across patients while maintaining the original TRmRNA expression matrix. For each permuted dataset, we applied the same LASSO and multivariate Cox regression pipeline to reconstruct the risk model and calculate its C‐index and AUC values. The distribution of performance metrics from the permuted models was then compared to that of the original model. A P‐value was calculated based on the proportion of permuted models whose performance exceeded that of the original model.

Creation of Nomogram to Predict OS

2.5

A nomogram was developed for a prognostic model, integrating clinical factors such as age, gender and stage to assess OS at 1, 3 and 5 years [13]. Additionally, we created a calibration plot that included the risk model and assessed the impact of the prognostic model on the consistency between survival prediction and actual outcomes by comparing its C‐index with that of the plot without the risk model. To compare the predictive power of the nomogram with that of the TNM stage alone, we performed statistical model comparisons using both the concordance index (C‐index) and the likelihood ratio test (LRT). Cox proportional hazards models were constructed using the coxph () and cph () functions in R, and the LRT was performed using the anova () function. A P‐value < 0.05 was considered statistically significant.

Enrichment Analysis

2.6

Based on outcomes derived from the expression of 9 TRmRNAs, the functional enrichment analysis was conducted to investigate functions linked to the identified genetic set, utilising Gene Set Enrichment Analysis (GSEA) 4.3.2 software. The P‐value was < 0.05 and the false discovery rate (FDR) was < 0.25 for screening suitable pathways [14].

Immune Checkpoint Expression

2.7

Using data from previous studies, we extracted information on immune checkpoints and analysed their expression levels in both HR and LR groups. The analysis was conducted using the “reshape2” package for data manipulation and processing [15].

Tumour Mutational Burden

2.8

Tumour somatic mutation data were collected from the TCGA and subsequently processed using the “TCGA biolinks” package. Waterfall plots were created in R utilising the “maftools” package, and the TMB was calculated [16, 17].

Evaluation of the TME and Immune Cell Distribution

2.9

To assess differences in stromal score, ESTIMATE score, immune score and tumour purity between the HR and LR groups, we utilised the “ESTIMATE” R package. Immune‐related data were further downloaded from the Tumour Immune Estimation Resource (TIMER) 2.0 website. A variety of computational tools were employed to refine the connection between risk scores and specific immune cell types [18]. For a more detailed analysis of immune cell composition, we performed immune infiltration analysis on 22 immune cell types. Furthermore, we used the single‐sample gene set enrichment analysis (ssGSEA) to compute function scores and compare these scores between the HR and LR groups.

Chemotherapy Assessment of TNF‐Related Signature

2.10

To explore tumour immune dysfunction and rejection (TIDE) in KIRC patients across both groups, we obtained TIDE scores and associated information from the TIDE website (https://tide.dfci.harvard.edu/) [19]. Furthermore, we calculated half‐maximal inhibitory concentrations (IC50) for standard chemotherapeutic drugs with the “oncoPredict” R package [20]. To evaluate differences in IC50 values among different groups, we employed the Wilcoxon signed‐rank test. Lastly, we compared drug responses between HR and LR groups based on the IC50 values.

Validation by Quantitative Reverse Transcription Polymerase Chain Reaction (RT‐qPCR) Three Tumour Cell Lines and One Normal Cell

2.11

KIRC cell lines, including 786‐O, A‐498, ACHN and the normal HK‐2 cell line, were obtained from Procell Life Science & Technology Co. Ltd. (Wuhan, China). RNA extraction was carried out using TRIzol reagent (Invitrogen, USA) and an RNA extraction kit. The RNA was reverse transcribed into cDNA using the PrimeScript RT Reagent Kit (Takara, Japan). TRmRNA expression levels, as detailed in Table S1, were quantified by RT‐qPCR, with relative RNA levels calculated using the 2^−ΔΔCT^ method [21].

Construct an External Validation Model of TRmRNAs

2.12

The E‐MTAB‐1980 data was retrieved from the ArrayExpress database, and we created ROC curves integrating clinical factors. Survival analyses were then conducted across patient subgroups to assess the prognosis of patients. Independent prognostic factors were identified through u/mCox. To estimate OS at 1, 3 and 5 years for KIRC, we utilised the “rms”, “timeROC” and “tidyverse” R packages, ensuring robust survival predictions [22].

Results

3

Screening for TRmRNAs in KIRC Patients

3.1

The rationale for our study is presented in a flowchart (Figure 1). First, we selected TRGs and developed protein–protein interaction (PPI) networks for 15 of these genes via the STRING database (https://string‐db.org/) (Figure 2A). Next, Pearson correlation analysis was conducted on TRGs, revealing significant associations between the expression of 17,846 mRNAs and TRGs (|Pearson R| > 0.3, p < 0.05). After refining these to 341 mRNAs, a differential expression analysis was performed, highlighting a significant difference in tumour tissue expression (Figure 2B).

Flow Chart.

PPI network (A); Filter for TNF‐related mRNAs. Volcano plot of differentially expressed TNF‐associated mRNAs (B); LASSO regression analysis (C, D).

Following the differential expression analysis, we randomly assigned 526 patients to test and train groups (Table 1). The uCox regression signature was employed to identify 78 mRNAs in the train cohort that were significantly related to differential prognosis (p < 0.05) (Table S2). Next, LASSO analysis was conducted, which led to the identification of 20 TRmRNAs linked to KIRC (Figure 2C,D) (Table S3). Subsequently, mCox regression analysis further reduced the list to 9 mRNAs that were significantly correlated with OS (Figure S1A). The interrelationships between these 9 mRNAs were explored (Figure S1B), along with their associations with TRGs (Figure S1C). Importantly, these mRNAs showed significant upregulation and downregulation in tumour tissues (Figure S1D).

Development and Verification of the Predictive Model

3.2

Using the nine TRmRNAs, we established a formula to calculate risk scores for each patient: Risk score = coefficient (SPTBN2) × expression (SPTBN2) + coefficient (SCNN1B) × expression (SCNN1B) + coefficient (GPC3) × expression (GPC3) + coefficient (DPEP1) × expression (DPEP1) + coefficient (FGF1) × expression (FGF1) + coefficient (MPP7) × expression (MPP7) + coefficient (PROX1) × expression (PROX1) + coefficient (SIM2) × expression (SIM2) + coefficient (ODF3B) × expression (ODF3B). Table 2 provides an overview of clinical characteristics for individuals in both HR and LR groups. In datasets for train, test and the total population, the HR group showed significantly poorer survival than the LR group (p < 0.001) (Figure 3A–C). The time‐dependent ROC curves verified strong model predictions, with Area Under the Curve (AUC) of 0.872 (95% CI: 0.813–0.932, p < 0.001), 0.823 (95% CI: 0.756–0.889, p < 0.001) and 0.838 (95% CI: 0.774–0.902, p < 0.001) for one, three and five years in the training group; 0.686 (95% CI: 0.590–0.781, p < 0.001), 0.678 (95% CI: 0.599–0.756, p < 0.001) and 0.718 (95% CI: 0.634–0.802, p < 0.001) in the test group; and 0.773 (95% CI: 0.718–0.828, p < 0.001), 0.726 (95% CI: 0.673–0.779, p < 0.001) and 0.765 (95% CI: 0.710–0.821, p < 0.001) in the overall group (Figure 3D–F). Additional visualisations, including the expression profiles of the 9 mRNAs in patients with different risk levels (Figure S2A–C), risk curves (Figure S2D–F), risk distribution plots (Figure S2G–I) and scatter plots (Figure S2J–L), all support the model's prognostic value.

The model prediction effect is validated by the train group, test group, and entire group. K‐M analysis (A–C) and time‐dependent ROC curves (D–F) to compare the survival of the high‐risk group and low‐risk group.

PCA analysis (Figure 4A–D) revealed clear separation between HR and LR groups. Using uCox regression, age, tumour stage and risk score were determined to be standalone prognostic variables (HR = 1.02, 1.98 and 2.27; 95% CI amount = 1.00–1.04, 1.65–2.37 and 1.82–2.82; p < 0.05, p < 0.001 and p < 0.001) (Figure 4E). Moreover, mCox regression validated that the stage and risk score continued to be significant independent factors of OS (HR = 1.85 and 1.83; 95% CI amount = 1.53–2.24 and 1.49–2.26; Both P values < 0.001) (Figure 4F).

PCA and independent prognostic analysis of the signature. PCA based on all genes (A), all mRNAs (B), TNF‐related mRNAs (C), and risk signature (D); univariate (E) and multivariate (F) independent prognostic analysis.

The heatmap, which incorporates various clinical features, risk groups and mRNAs associated with the prognostic model (Figure S3A), along with survival curves for patients across different clinical subgroups (Figure 5), demonstrates the robust predictive power of the model. Significant differences in survival were observed among various clinical categories. To confirm the robustness of our model and exclude the possibility that its performance arose by random chance, we conducted a permutation test with 1000 iterations. The results revealed that the original model's C‐index and AUC values were significantly higher than those obtained from the permuted datasets (p < 0.001), suggesting that the predictive ability of the model is unlikely to be due to chance. A visual comparison of the performance distributions is provided in Figure S11.

Further validation of model effects. Survival curves of patients in different clinical states (A–F).

Nomogram for Predicting Clinical OS

3.3

Evaluation of the link between risk score and clinical factors indicated that age (p = 0.51) showed no significant relationship with the risk score, whereas clinical stage (p < 0.001) exhibited a strong link (Figure S3B,C). Decision curve analysis (DCA) and ROC curve analysis (Figure 6A,B) showed that the risk score offered superior predictive value compared to other clinical factors. A prospective estimator for KIRC patients, which integrated clinical factors, was developed (Figure 6C); the AUC for the risk score in predicting 3‐year overall survival (OS) was 0.766, which was comparable to TNM stage (AUC = 0.812) and superior to age (AUC = 0.649) and sex (AUC = 0.492), indicating that the risk score provides robust prognostic value, and its validation using patient data yielded promising results (Figure S4). To formally assess the additive prognostic value of the nomogram, we compared the combined model (TNM stage + TRmRNA‐based risk score) with a model based solely on the TNM stage. The likelihood ratio test confirmed that the addition of the risk score significantly improved model fit (χ ^2^ = 37.15, df = 1, p = 1.10 × 10^−9^), suggesting that the TRmRNA risk score provides meaningful prognostic information beyond conventional clinical staging (Figure 6D,E).

Nomogram predicts patient prognosis. Decision curve to test for forecast value (A) ROC curves containing different clinical information (B) A clinical prognosis nomogram is constructed by age, gender, risk, and stage together (C). Nomogram with (D) and without (E) risk model.

Enrichment Analysis

3.4

The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed functions of model‐related genes, such as neuroactive ligand–receptor interaction, cell adhesion molecules, collecting duct acid secretion, viral protein interaction with cytokines and cytokine receptors, cytokine–cytokine receptor interaction, protein digestion and absorption and phagosome formation (Figure S5). GSEA further revealed distinct pathway enrichments across various gene sets. In the HR group, the GOBP library showed significant enrichment in functions related to immunoglobulin production, B cell‐mediated immunity, lymphocyte‐mediated immunity, molecular mediator production for immune response, immunoglobulin complex formation and antigen binding (Figure S6). A detailed list of the enriched pathways can be found in Table S4.

Tumour Mutation Burden

3.5

TMB was evaluated in patients, revealing a significantly higher TMB in the HR (Figure 7A). To further investigate the genetic landscape, a waterfall plot visualised the 20 most mutated genes across subgroups (Figure 7B,C). Among the HR group, the most commonly mutated genes were VHL (40%), PBRM1 (32%), BAP1 (19%), TTN (19%) and SETD2 (18%). In contrast, LR patients showed the highest mutation frequencies in VHL (51%), PBRM1 (49%), TTN (14%), SETD2 (7%) and MUC16 (8%).

Tumour mutation burden in different risk groups. Percentage bar graph showing TMB for different risk subgroups (A) High‐risk group waterfall chart (B) Low‐risk group waterfall chart (C).

Tumour Immune Infiltration Status Analysis

3.6

TME analysis identified elevated ESTIMATE scores (p < 0.05; Figure 8A) and Immune scores (p < 0.001; Figure 8C) in the HR group, while the Stromal score showed no significant variation (p > 0.05; Figure 8B). Also, the HR group displayed reduced tumour purity (p < 0.05; Figure 8D). These observations highlight that HR is linked to an immune‐suppressive tumour microenvironment, potentially impairing the antitumor immune response. We also examined the connection between immune cell composition and risk scores using multiple algorithms (Figure 8E). Figure 8F exhibited that in the HR group, there is a higher proportion of T cells CD4 memory resting (p < 0.001). In contrast, the LR group showed an increased proportion of NK cells activated (p < 0.001) and Macrophages M2 (p < 0.001). The connection between immune cells and risk scores is illustrated (Figure S7). Functional pathway analysis revealed a reduction in type II IFN response (p < 0.001) and APC coinhibition (p < 0.001) in HR patients, whereas pathways such as T cell stimulation (p < 0.001) and type I IFN response (p < 0.05) were upregulated (Figure 8G).

*Analysis of tumour immune microenvironment. Violin plots of differences in ESTIMATE scores (A), stromal scores (B), immune scores (C) and tumour purity (D) for different risk subgroups; Bubble plots of correlations between immune cells and risk scores under six algorithms (E); Proportions of 22 immune cells in two subgroups under the CIBERSORT algorithm (F); single sample gene set enrichment analysis (G). *p < 0.05, **p < 0.01, **p < 0.001.

Benefits of Predictive Models for KIRC Treatment

3.7

TIDE revealed that both T cell exclusion and T cell dysfunction were significantly more pronounced in the HR group with TNF‐associated signatures (Figure S8A–C). This suggests that the HR subgroup has a higher potential for immune evasion, which could complicate therapeutic interventions. Further research identified 19 differentially expressed immune checkpoints including PPP6R1, IL15RA and DDX47 that were upregulated in the HR group (Figure S8D), which indicates that the HR group may benefit from immune therapies targeting these specific immune checkpoints. Additionally, we assessed the sensitivity of tumours in both the LR and HR groups to various anti‐tumour drugs. Drugs exhibiting differential efficacy between the subgroups are listed in Table S5 (LR) and S6 (HR), while those with no notable differences are shown in Table S7. The HR patients notably revealed heightened sensitivity to therapies aimed at the PI3K/MTOR signalling pathway, which may guide the decision of targeted therapies.

Risk Model Validation In Vitro

3.8

Immunohistochemical images, obtained from the HPA website (Figure S9A), illustrated the differences in protein expression of TRGs in KIRC and nearby normal tissues. RT‐qPCR analysis demonstrated that, of nine TRmRNAs, DPEP1, FGF1 and MPP7 were largely expressed in normal tissues, whereas the other six mRNAs exhibited elevated expression levels in tumour tissues (Figure S9B).

Validation of the Prognostic Traits Externally

3.9

To assess the predictive value of the signature, we developed an external validation using the E‐MTAB‐1980 cohort [23]. K‐M analysis found that patients in the HR group experienced worse outcomes than those in the LR group (p = 0.003; Figure S10A). The model's robustness was further confirmed through ROC curve analysis, which produced AUC of 0.871, 0.776 and 0.732 for 1, 3 and 5 years (Figure S10B). Importantly, the AUC for the risk score (0.896) was notably higher than those for age (0.573) and gender (0.627), highlighting the excellent prognostic capability of the risk score model (Figure S10C). To validate the TRmRNA signature's predictive value across various clinical subgroups, uCox regression models were performed, with stage (HR = 2.28, 95% CI amount 1.68–3.10, p < 0.001) and risk score (HR = 2.17, 95% CI amount 1.66–2.84, p < 0.001) (Figure S10D). The mCox analysis had similar results: stage (HR = 2.03, 95% CI amount 1.42–2.90, p < 0.001) and risk score (HR = 1.80, 95% CI amount 1.31–2.45, p < 0.001) (Figure S10E). The heatmap demonstrated a clear survival difference between the LR and HR groups, with the LR group showing significantly longer OS (Figure S10F). Despite potential cross‐platform variability between RNA‐seq (TCGA) and microarray (E‐MTAB‐1980), our model exhibited consistently high prognostic performance, suggesting that the core biological signals captured by TRmRNAs are stable and generalisable across platforms.

Discussion

4

Over the past few years, KIRC incidence has steadily increased, posing significant challenges in prognosis prediction and therapeutic decision‐making. Despite advancements in surgical and systemic treatments, the heterogeneity in patient outcomes highlights the limitations of current staging systems, such as the AJCC TNM classification, which fail to capture the molecular complexity of tumours. In response to these challenges, biomarker‐based prognostic models have emerged as promising tools to enhance prediction accuracy [3]. While TNF‐related gene signatures have been applied in malignancies such as cervical squamous cell carcinoma and bladder cancer, their applicability to kidney renal clear cell carcinoma (KIRC) remains limited [24]. KIRC exhibits distinct biological and immunological hallmarks—particularly immune evasion and pronounced vascularisation—that necessitate cancer‐specific prognostic modelling. In our research, we developed and validated a predictive signature utilising nine TRmRNAs, which demonstrated extremely high accuracy in assessing prognosis in HR and LR groups. Additionally, we integrated this molecular model with clinical factors into a nomogram, providing a more comprehensive approach for predicting survival outcomes. By incorporating enrichment analysis, TMB evaluation and drug sensitivity predictions, our study provides important perspectives for the development of personalised treatment strategies in KIRC.

In this study, we propose a novel prognostic model for kidney renal clear cell carcinoma (KIRC) based on a refined subset of TNF family‐related mRNAs. Unlike previous models—including those by Huang et al. and Zhang et al. [25, 26]. That employed broad TNF gene panels in lung adenocarcinoma—our approach tailors gene selection to the distinct immune and molecular landscape of KIRC. This cancer‐specific focus enhances both prognostic precision and clinical relevance. This model exhibited high predictive accuracy across multiple evaluation metrics, including the ROC curve, C‐index and calibration curve. Notably, our risk model demonstrated prognostic power comparable to that of the TNM staging system and outperformed age and sex in ROC analysis. When integrated into a nomogram with standard clinical variables, the inclusion of the TRmRNA risk score significantly improved predictive accuracy, as reflected by a higher C‐index. These results notably surpassed the performance of the TNF‐related model proposed by Ma et al. [27]. Our approach, utilising TNF‐related genes to construct prognostic signatures, aligns with similar studies in cervical and bladder cancers, further validating the robustness of this method [28]. Several key mRNAs in our model, such as PROX1 and ODF3B, are well‐established in tumour biology, further reinforcing the biological relevance of our prediction framework. Specifically, Lv et al. have shown that PROX1 promotes KIRC cell proliferation [29]. Also, ODF3B, a gene closely linked to glioma growth and apoptosis, has been reported to exhibit elevated expression levels in gliomas. Notably, higher expression of ODF3B is correlated with a poorer prediction in glioma patients [30]. SCNN1B is downregulated in KIRC due to promoter hypermethylation and may preserve epithelial integrity [31]. GPC3 and MPP7 inhibit proliferation and invasion through modulation of growth factor signalling and maintenance of epithelial polarity, respectively [32, 33]. DPEP1 is involved in immune‐metabolic regulation, and its low expression is linked to poor prognosis. Conversely, some TRmRNAs appear oncogenic [34]. FGF1 promotes RCC proliferation via PI3K/Akt and MAPK pathways [35]. SIM2 is linked to increased invasion and Wnt pathway activation [29, 36]. Additionally, DCA indicated that our risk score outperforms traditional prognostic predictors, offering a superior net benefit in risk assessment [37]. Together, the TNF‐related prediction model demonstrates potential as an effective and dependable approach to enhance prognostic assessment for cancer patients.

To further investigate the prognostic differences and molecular alterations between HR and LR groups, we conducted pathway enrichment and immune function analysis in the tumour microenvironment using our predictive model. Functional enrichment analysis showed that the HR group displayed increased expression of pathways related to immunity, indicating an activated but dysregulated immune response [38]. Among these pathways, the enrichment of CD22‐mediated BCR regulation in the high‐risk group may represent an immunosuppressive mechanism mediated by dysfunctional or regulatory B‐cell populations. CD22 negatively regulates BCR signalling and contributes to immune tolerance. Its increased activity in high‐risk tumours could suppress antigen presentation and cytotoxic immune responses, supporting immune escape and tumour progression [39]. This immune dysregulation may contribute to the more aggressive tumour behaviour observed in the HR group. Conversely, the LR group showed enrichment in solute carrier and proximal tubule transport pathways; these pathways are characteristic of functional proximal tubular epithelial cells and suggest that tumours in this group may retain more segment‐specific functions. The preservation of such renal epithelial programs may reflect a less aggressive tumour biology and explain the improved outcomes observed in this subgroup [40]. In addition, immune function analysis detected notable differences between the two risk groups in processes like antigen‐presenting cell (APC) co‐suppression, parainflammatory responses and type II interferon signalling. Notably, IFN‐γ has emerged as a pivotal mediator in the immune microenvironment of tumours. It is identified as a significant marker of immunotherapy effectiveness and a key factor in regulating tumour progression [41]. The high TMB observed in the HR group further supports the potential for immunotherapy, as TMB is typically associated with an increased burden of neoantigens. However, it also correlates with poorer prognosis in KIRC patients [42, 43]. Additionally, the inactivation of tumour suppressor genes such as VHL has been implicated in cancer progression. Loss of VHL function facilitates histone lactylation, which activates the expression of platelet‐derived growth factor receptor β (PDGFRβ) at the transcriptional level, a key factor in KIRC development [44]. Overall, the poorer prognosis observed in the HR group is closely linked to elevated TMB, dysregulated immune‐related pathways and inactivation of tumour suppressor genes. These findings underscore important mechanisms driving KIRC progression and provide important guidance for identifying potential therapeutic targets and research directions.

As previously discussed, TMB is a critical determinant of immunotherapy efficacy in cancer treatment. Studies suggest that patients with high TMB are prone to respond well to immunotherapy [45]. In our study, the HR group exhibited elevated TMB levels, indicating that these patients may be more responsive to immunotherapeutic interventions. To further investigate the potential for personalised treatment approaches, we compared the immunotherapy responses between the HR and LR groups; although these drug sensitivity predictions were based on in silico IC50 estimation from transcriptomic data, rather than real‐world clinical responses, they provide a preliminary pharmacogenomic framework that warrants further experimental and clinical validation. Notably, HR patients exhibited heightened responsiveness to therapies aimed at the PI3K/MTOR signalling pathway, such as Dactolisib and Buparlisib, consistent with the drug sensitivity findings, suggesting elevated pathway activity that may contribute to both poor prognosis and enhanced sensitivity to pathway‐specific inhibitors. Notably, Buparlisib has been shown to downregulate epithelial‐mesenchymal transition pathways and matrix metalloproteinases, thereby limiting metastasis in urinary tract cancers, indicating promising therapeutic potential [46, 47]. These findings align with the recognised role of PI3K/MTOR signalling in developing tumour growth and conferring resistance to conventional therapies [48, 49], likely due to aberrant activation of this pathway. Ye et al. demonstrated that GNE‐477 effectively suppressed RCC growth by blocking PI3K/AKT/mTOR signalling [50]. The enrichment of pro‐survival and immunosuppressive features in high‐risk patients suggests greater pathway dependency, supporting their potential responsiveness to targeted inhibition. In contrast, LR patients demonstrated higher responsiveness to agents such as Cediranib and Osimertinib, which specifically target the VEGFR and EGFR pathways, respectively [51]. This variation in drug sensitivity between the two groups reflects distinct molecular and immune characteristics associated with each risk profile. Additionally, the elevated expression of immune checkpoints, including IL15RA and DDX47, in the HR group suggests that these individuals might benefit from immune checkpoint inhibitors (ICIs). DDX47, an RNA‐binding protein, has been closely linked to KIRC progression [52, 53]. Recent studies have highlighted the effectiveness of ICIs in treating KIRC, particularly in patients with high TMB or increased expression of immune checkpoints [54, 55]. These observations emphasise the potential of integrating risk‐based stratification with immunotherapy strategies, which could significantly improve therapeutic outcomes for KIRC patients.

Our study makes several key contributions to the field of KIRC research. Firstly, our bioinformatics‐based pipeline is in line with recent prognostic studies in oncology, such as those reported by Dang et al. and Tran et al., supporting the methodological soundness and translational potential of our approach [56, 57]. We present a prognostic model that demonstrates high predictive accuracy, surpassing similar models, thus establishing TRmRNAs as promising biomarkers for KIRC prognosis. Secondly, we strengthen the validity of our findings by performing external validation using independent datasets, which enhances both the biological relevance and reliability of our results [58]. While the final model includes only nine TRmRNAs, we acknowledge the potential risk of overfitting due to the limited sample size. To reduce this risk, we applied stringent feature selection (univariate Cox, LASSO and ML‐based filtering), and more importantly, validated the model in an external cohort (E‐MTAB‐1980), where it achieved robust predictive performance. These results support that the model captures reproducible prognostic signals rather than overfitting noise and demonstrate its generalisability across independent datasets. Thirdly, our study is the first to investigate the link between TRmRNAs and KIRC prognosis, highlighting these mRNAs as potential biomarkers for the disease. Although our TRmRNA‐based prognostic model was validated in two independent public datasets, the lack of prospective clinical validation remains a key limitation. As this study is retrospective in nature, future work will involve prospective, multicentre cohorts to evaluate the model's clinical utility in risk stratification and treatment guidance. While no direct functional experiments were performed, our results provide a hypothesis‐generating framework. Follow‐up studies will include genetic perturbation and phenotypic assays in cellular and animal models to elucidate the biological roles of key TRmRNAs. Notably, external validation using the E‐MTAB‐1980 dataset supports the model's robustness. Nevertheless, residual batch effects cannot be fully excluded, and overfitting from random train‐test splits may still occur. Although feature selection and model training were strictly limited to the training cohort, future validation strategies will consider cohort‐level or subtype‐aware designs to minimise data leakage. In conclusion, our model shows strong prognostic potential, but further experimental and clinical validation is essential to confirm its applicability in KIRC.

Conclusion

5

Our risk model, developed based on 9 TRmRNAs, demonstrates strong predictive accuracy for patient outcomes. When integrated with clinical factors, the accompanying nomogram further enhances its predictive power, offering exceptional effectiveness in forecasting survival. The poor prognosis observed in HR patients is strongly correlated with factors such as elevated TMB, enrichment of CD22‐mediated BCR regulation and immunoglobulin complex pathways, etc. The results of drug susceptibility analysis provide valuable information that can help in the creation of personalised treatments. However, additional clinical trial validation is essential to fully realise its applicability and impact in clinical practice.

Author Contributions

Zijian Hu: conceptualization (equal), data curation (equal), formal analysis (equal), investigation (equal), methodology (equal), project administration (equal), resources (equal), software (equal), supervision (equal), validation (equal), visualization (equal), writing – original draft (lead), writing – review and editing (equal). Yajie Zhou: conceptualization (equal), data curation (equal), formal analysis (equal), investigation (equal). Lei Xie: conceptualization (equal), data curation (equal), formal analysis (equal), investigation (equal), methodology (equal). Shuwen Zhang: conceptualization (equal), data curation (equal), formal analysis (equal), investigation (equal), methodology (equal), writing – original draft (equal). Yijiang Liu: conceptualization (equal), data curation (equal), formal analysis (equal), investigation (equal), methodology (equal), writing – original draft (equal). Wenxiong Zhang: conceptualization (equal), data curation (equal), formal analysis (equal), funding acquisition (lead), investigation (equal), methodology (equal), project administration (equal), resources (equal), software (equal), supervision (equal), validation (equal), visualization (equal), writing – original draft (equal), writing – review and editing (lead). Ting Huang: conceptualization (equal), data curation (equal), formal analysis (equal), investigation (equal), methodology (equal), project administration (equal), resources (equal), software (equal), supervision (equal), validation (equal), visualization (equal), writing – original draft (equal), writing – review and editing (equal).

Ethics Statement

This article does not contain any studies with human participants or animals performed by any of the authors.

Consent

For this type of study formal consent is not required.

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Figure S1. Construct a prognostic model. Forest plots associated with multifactor regression (A); Correlation circle plot of 9 TNF‐related mRNAs (B) Mulberry plots of 9 TNF‐related mRNAs with TNF‐related genes (C) Deviation plots indicating upregulation and downregulation changes of 9 TNF‐related mRNAs (D).

Figure S2. The model prediction effect is validated by the train group, test group and the entire group. Heat map of 9 TNF‐related mRNAs expressions (A‐C). Risk curve for risk scores (D‐F) and Scatterplot (G‐I) for the survival status of each patient; Scatterplot of risk scores of patients with different survival statuses (J‐L).

Figure S3. Further validation of model effects. mRNA signature expression heat map with clinical information (A); Correlation analysis of risk signature with age (B) and stage (C).

Figure S4. Nomogram with patient verification.

Figure S5. KEGG enrichment analysis. Bar chart (A); Circle chart (B).

Figure S6. GSEA enrichment analysis. Enrichment pathways in entire group in different gene sets (A‐F).

Figure S7. Correlation plot of risk scores with immune cells. Scatter plot of the correlation between risk scores and cancer associated fibroblast MCPCOUNTER (A), cancer associated fibroblast EPIC (B), T cell NK (C), T cell regulatory (Tregs) (D), Neutrophil (E), T cell (F), Myeloid dendritic cell (G), T cell follicular helper (H), T cell CD4+ Th1 (I), Macrophage M1 (J), Mast cell activated (K), NK cell (L), Mast cell activated (M), T cell follicular helper (N).

Figure S8. The TNF‐related mRNAs risk pattern in tumour therapy. Differences in dysfunction (A), immune exclusion (B), TIDE (C) and immune checkpoints (D) in high and low‐risk groups. * p < 0.05, ** p < 0.01, *** p < 0.001.

Figure S9. In vitro experimental validation of the risk model. Immunohistochemical staining images of partial TNF‐related gene proteins in KIRC tissue and normal tissue (A); Relative expression of 9 TNF‐related mRNAs in different risk subgroups (B). * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001.

Figure S10. External validation of prognostic characterisation. K‐M analysis (A) and Time‐dependent ROC curves (B) to compare the survival of the high‐risk group and low‐risk group. ROC curves containing different clinical information (C). Univariate (D) and multivariate (E) independent prognostic analysis. The heatmap between high‐risk and low‐risk group (F).

Figure S11. Results of permutation testing (1000 iterations). The histogram shows the distribution of C‐index values from permuted models. The red dashed line indicates the C‐index of the original model, which lies far in the tail of the null distribution (p < 0.001), confirming that the model’s performance is unlikely to occur by chance.

Table S1. Primer sequences for 9 TNF‐related mRNAs. mRNA signature expression heatmap (C).

Table S2. 78 differentially expressed TNF‐related mRNAs.

Table S3. 20 differentially expressed TNF‐related mRNAs.

Table S4. GSEA pathways for different risk groups.

Table S5. Antineoplastic drug sensitivity information (sensitive group: low).

Table S6. Antineoplastic drug sensitivity information (sensitive group: high).

Table S7. Antineoplastic drug sensitivity information (no obviously sensitive group).