Data for transcriptome RNA sequencing and clinical information for 598 ccRCC cases (72 normal samples and 526 tumor samples) were sourced from the TCGA (The Cancer Genome Atlas) database (https://xenabrowser.net/datapages/) (17). Cases with clinical data and the overall survival (OS) missing were excluded. Meanwhile, two independent databases: ICGC (International Cancer Genome Consortium) database (https://icgc.org/) containing 91 ccRCC patients; and E-MTAB-1980 database (https://www.ebi.ac.uk/arrayexpress/) (18) comprising 101 ccRCC patients, were retrieved to validate our findings.

A group of 69 genes linked with NETs were extracted from prior studies (11, 15, 16, 19, 20) and shown in Table S1 . The differently expressed genes (DEGs) associated with NETs between tumor and normal tissues were identified using the “DESeq2” package in R. Candidate NETs genes were obtained after setting the fold change (|fold change| ≥ 1) and adj. P < 0.05. The “pheatmap” R package was utilized to generate heatmap of the DEGs. Then, STRING (Search Tool for the Retrieval of Interacting Genes; Version 11.0, https://string-db.org/) was used to acquire a protein −protein interaction (PPI) network.

Initially, we assessed the prognostic significance of the NETs-associated genes using univariate Cox regression analysis, with a significant filter set at 0.05, employing the “survival” R package. Subsequently, based on the results obtained from the aforementioned analysis, we utilized the “glmnet” R package to develop a prognosis model, utilizing LASSO (least absolute shrinkage and selection operator). The following formula was used to calculate the risk score: risk score = β (1) gene (1) × expression of gene (1) + β (2) gene (2)× expression of gene (2)+…+ β(n) gene(n) × expression of gene(n). Patients were classified into high-risk (>median number) and low-risk (≤median number) groups based on the median value of risk score. The “survival” and “survminer” R package were used to conduct the Kaplan–Meier analysis (21, 22). The R package “survival ROC” was used to examine the predictive accuracy at 1-, 3-, and 5-year intervals through a time-dependent receiver operating characteristic (ROC) curve analysis.

Clinical parameters (age, gender, and tumor stage) of the patients in the TCGA cohort were extracted. Together with the risk score, these clinical variables were analyzed by univariate and multivariate Cox proportional hazards regression. The Schoenfeld residuals were calculated for each variable to assess if each variable independently satisfied the assumptions of the Cox model. By using the “rms” and “survival” R packages, a nomogram was further depicted to predict OS in TCGA patients (1 year, 3 years, and 5 years) based on the previously mentioned NETs signature and clinical prognostic factors. Besides, we have performed the decision curve analysis to compare the prognostic value of models including and not including risk score.

The TCGA cohort of ccRCC patients were divided into two subgroups, characterized as per their risk scores. DEGs were separated using FDR<0.01 and |log2FC| ≥ 1.5 filtering criteria between the high- and low-risk group. Subsequently, we executed enrichment analyses (GO and KEGG analysis) with the aid of the “ClusterProfiler” R package to identify potential correlation between the enrichment terms.

The TCGA patients were divided into two groups according to their risk scores. The assessment of tumor microenvironment for both groups along with stromal score, immune score and ESTIMATE score was conducted using the “ESTIMATE” R package (23). To quantify the immunologic cell abundances in the immune microenvironment, the infiltration levels of 28 immune cell types were evaluated by the enrichment score, which was calculated via the single-sample gene set enrichment analysis (ssGSEA) method from “GSVA” R package (24, 25). Furthermore, we gathered a group of key immune checkpoints and compared their different expressions between two NET risk groups.

The GDSC (Genomics of Drug Sensitivity in Cancer) database provides information regarding the sensitivity of cancer cell lines to diverse anticancer drugs (26), which was used to predict drug susceptibility for TCGA-KIRC patients. Drug sensitivity were calculated with the “oncoPredict” R package (27). Then we compared the therapeutic ability of several common targeted therapies in KIRC between two NET risk groups.

R software (version 4.3.0) and SPSS software (version 23.0) were employed for all statistical analyses. Box plot analyses were conducted using the Wilcoxon rank-sum test. Kaplan-Meier method was utilized to generate survival curves for high-risk and low-risk groups of patients. The significance of differential groups was assessed using the log-rank test. All hypothetical analyses were two-sided, where a P value < 0.05 was deemed significant.

Expression levels of 69 genes related to NETs were compared between 526 tumors and 72 normal tissues in the database TCGA-KIRC, leading to the identification of 31 DEGs. Six genes (DNASE1, CYP4F3, MME, KCNJ15, MTOR, SELP) were down-regulated and 25 other genes (SLC22A4, SLC25A37, CD93, FCAR, PADI4, FPR2, FCGR3B, S100A12, PTAFR, DYSF, CREB5, FPR1, TLR2, SIGLEC14, CEACAM3, CYBB, ITGAM, SELPLG, TLR7, TLR8, LILRB2, ITGB2, VNN3, CSF3R, MMP9) were up-regulated in tumor tissue samples. RNA expression of these genes was visualized through a heatmap ( Figure 1A , blue: low level of expression; red: high level of expression). To explore potential interactions among these NETs-related genes, PPI analysis was conducted by applying the STRING platform, and the results were shown in Figure 1B .

To identify the significant NETs-related genes for TCGA patients, univariate Cox regression analysis was contacted to evaluate their prognostic value. 17 genes were selected after filtered by P<0.05. The forest plot represented the HR of each single gene ( Figure 2A ). LASSO regression analysis was utilized to construct a prognostic model, and nine genes along with their coefficients were ultimately included in the NETs-related gene prognostic signature, based on the optimal λ value. ( Figures 2B, C ). The risk score was calculated as following, NET risk score = 0.003 × MMP9 + 0.08 × VNN3 + 0.386 × SLC25A37 - 0.154 × CD93 -0.201 × KCNJ15 - 0.121 × SLC22A4 + 0.152 × LILRB2 + 0.049 × TLR2 + 0.020 × FPR2. Then, the 519 TCGA patients were classified into high- and low- risk groups based on the median value of risk score ( Figure 2D ). The Kaplan-Meier curves demonstrated that patients in the high-risk group had a significantly shorter OS compared to those in the low-risk group ( Figure 2E , P<0.0001), and the AUCs for 1-, 3- and 5-year OS rates were found to be 0.71, 0.71 and 0.75, respectively ( Figure 2F ).

Subsequently, we used ICGC and E-MTAB-1980 cohorts as two independent validation datasets and observed a similar outcome in the high risk group of ccRCC patients ( Figures 2G–J ).

To assess the independence of NET risk score, the univariate and multivariate Cox regression analyses for the NETs-related signature and other clinical pathological factors associated with OS were performed. Results were presented in Table 1 , revealing that the risk score was deemed an independent risk factor for the prognosis of TCGA patients with ccRCC (HR: 2.53, 95% CI: 1.77–3.61, P <0.001). Besides, results from the Schoenfeld residual revealed that each variable in the multiple Cox model independently satisfied the assumptions of the Cox model ( Figure S1 ).

By comparing the AUC among different models, we can determine whether the NET risk signature in this study improves the performance in terms of OS prediction. ROC analysis showed that the model comprised of NET risk signature and clinical variables (age and stage) have the higher AUC value at 1-, 3- and 5-year than the other model only including clinical variables, confirming the valuable contribution of NET risk signature ( Figure S2 ).

In addition, Kaplan–Meier survival curves for nine NETs genes and correlation between their expression and clinical characteristics were presented in Figures S3 , S4 , respectively, to shed light on the role of prognostic signature genes in the development of ccRCC. In brief, except KCNJ15 was down regulated in tumor, the other eight genes were all up regulated. KCNJ15 together with CD93 and SLC22A4 were decreasingly expressed in the higher grade of tumor. Patients with low expression of KCNJ15, CD93 and SLC22A4 had shorter OS. All together, KCNJ15 may act as a tumor suppressor in ccRCC progression.

To enhance the predictive capacity of the signature for ccRCC patients, a nomogram developed by incorporating available clinical pathological parameters and risk score. Through the multivariate Cox regression analysis, age, stage, and risk score were considered as independent factors and were included in the nomogram ( Figure 3A ). Moreover, calibration curves were plotted for 1-, 3-, and 5-year survival rates to assess the accuracy of the nomogram ( Figure 3B ). The conspicuous agreement emerged for the OS predicted by the nomogram and the reliable values across the following period. The stability and accuracy of the nomogram including our NETs-related genes signature as well as age and tumor stage can predict the outcome of individual patients. In addition, the decision curve model showed that the prognostic model with risk score had a good performance ( Figure 3C ). The association between signature risk scores and clinical pathological characteristics was presented in the form of a Sankey diagram ( Figure 3D ). A higher signature risk score was significantly correlated with higher tumor stage (p < 0.001).

GO and KEGG analyses were conducted to assess the biological involvement of the DEGs. As presented in Figure 4A , the top GO terms comprised acute inflammatory response, lipid catabolic process, collagen-containing extracellular matrix, etc. Furthermore, KEGG analysis indicated that the DEGs were primarily involved in pathways such as complement and coagulation cascades, staphylococcus aureus infection and estrogen signaling pathway, as shown in Figure 4B .

The association between NET risk classification and the immune status in patients with KIRC was assessed. We found that KIRC patients in low risk group had notably elevated ESTIMATE scores and immune scores when compared to those in high risk group ( Figure 5A ). These results suggest an inverse relationship between NET scores and immune status in KIRC patients. Then, we compared the immune cell infiltration in the high versus low NET risk groups based on”ssGSEA” algorithm. The findings revealed that the immunosuppressive cells, such as myeloid-derived suppressor cells (MDSC) and regulatory T cells (Tregs), were significantly higher in the NET high risk group ( Figure 5B ). This phenomenon may imply the immunosuppressive microenvironment in the NET high risk group. In addition, the role played by the 9 NETs genes in immune cell infiltration were displayed separately ( Figure S5 ).

We figured out the disparate expression of immune checkpoints between low and high NET risk groups. Except TIM-3, the expression of majority selected immunosuppresive checkpoints PD-1, CTLA4, LAG3, A2BR, and B7-H3 were all significantly upregulated in the high risk group ( Figure 5C ). Moreover, a significant association between NETs risk group and the efficacy of targeted therapies such as gefitinib, axitinib, afatinib, erlotinib, saptinib, and ibrutinib was detected. Specifically, patients with lower NETs-scores had a better response to these targeted therapies ( Figure 5D ).