Identification of 27 differentially expressed prognostic PCDs and PRLs in KIRC patients

The screening criteria led to the selection of 50 differentially expressed genes, 27 of which were associated with prognosis (Fig. 2A and D). It is shown in Fig. 2B that 27 intersection genes have different expression levels between KIRC tissues and normal kidney tissues. The interactive information of 27 intersection genes is displayed in Fig. 2C. We then screened 659 lncRNAs that were closely related to the 27 PCD genes (Fig. 3A, Supplementary Table S3), and 491 LncRNAs with considerably differential expression levels were identified (Supplementary Fig. 1, Table S4).

Fig. 2

Identification of the prognostic PCD-related DEGs in the KIRC cohort. A Venn diagram showing the overlapping genes between PCD-related DEGs and OS-related genes; B heatmap showing the candidate genes expression; C the correlation network of candidate genes; D forest plots showing the candidate genes identified by univariate Cox regression analysis

Fig. 3

Identification of a prognostic PRL signature in KIRC. A Sankey diagram of the relationship between PCD and PRLs; B LASSO coefficient distribution of PRLs; C variable selection cross-validation in the LASSO analysis; D coefficients of the seven chosen PRLs

Developing and validating a PRL prognostic model for KIRC patients

Univariate Cox proportional hazards regression analyses was employed to determine 147 lncRNAs that were considerably related to the OS (p < 0.05), and then multivariate Cox analysis and LASSO analysis were conducted to screen the PRLs. Seven PRLs were used to establish a prognostic risk model: LINC02747, AP001636.3, AC022126.1, LINC02657, long intergenic non-protein coding RNA 2609 (LINC02609), lncRNA associated with SART3 regulation of splicing (LINC02154), and ZNF706 neighboring transcript 1 (ZNNT1). The risk scores of the seven lncRNAs were calculated using their expression levels and regression analysis coefficients (Fig. 3B–D). Risk was calculated as follows: Risk score = LINC02747*(-0.160234181791124) + AP001636.3*(-2.21600604144298) + AC022126.1*(1.07480838933924) + LINC02657*(0.806615089929543) + LINC02609*(0.34964876422504) + LINC02154*(0.478170693602118) + ZNNT1*(0.447721371873328). Low-risk patients in the training set had a notably lower mortality rate than high-risk patients (Fig. 4A and D). Figure 4G illustrates the expression patterns of the seven PRLs in the training set. In the training set, Fig. 4J reveals that the low-risk patients are more likely to survive than high-risk patients. The prognostic risk model was then validated in the testing set (Fig. 4B, E, H, and K), overall set (Fig. 4C, F, I, and L), ICGC set (Fig. 5A), and results from them were highly consistent with those from the training set.

Fig. 4

Prognostic analysis of the risk model in different sets. A–C The overall survival risk score distribution; D–F the distribution of survival time and survival status; G–I heatmaps of the expressions of the seven PRLs; J–L survival outcomes of different sets

Fig. 5

Accuracy assessment of the risk characteristics. A Validation of the prognostic signature for KIRC patients in ICGC cohorts; B validation of prognostic model effectiveness in ICGC cohorts; C evaluation of prognostic model effectiveness in training; D testing, E overall; F–H comparison of single prognostic factor and nomogram prognostic efficacy

Analysis of the clinical characteristics of the PRL prognostic risk model

In order to evaluate the model’s accuracy, ROC curves and AUC values were calculated. As of the training set, AUC values were 0.750, 0.770, and 0.811 for 1-, 3-, and 5-year OS (Fig. 5C), which indicates that the prognostic risk model we constructed had a promising predictive ability in the training set. Moreover, ICGC data were used to validate the model’s predictive ability (Fig. 5B), the testing dataset (Fig. 5D) and the overall dataset (Fig. 5E). This study also examined whether a prognostic signature could be more accurate at predicting prognosis than common clinicopathological characteristics through multivariate ROC analysis. Based on the training set, AUC value for risk score was 0.751, which was higher than that of age (0.672), gender (0.499), and grade (0.725), but lower than stage (0.805) (Fig. 5F). AUC values for the testing and the overall set are shown in Fig. 5G and H, respectively, which were similar to the training set, indicating that the combination of the prognostic signature and TNM stage may be more valuable for diagnosis. Additionally, as shown in Fig. 6A and B, each of the p values for age, stage, grade, and risk score was less than 0.05, indicating that they could be considered independent variables. In order to facilitate clinical decision-making, a nomogram was constructed based on clinical characteristics and the risk score. Clinicopathological features including age, metastasis, gender, as well as AJCC T stage, grade, stage, and the risk score are all shown in Fig. 6C. In ccRCC patients, the nomogram predicting OS was highly accurate according to the calibration curve.

Fig. 6

Development of a nomogram by integrating the risk score and clinicopathological features in the KIRC cohort. A Univariate Cox analysis; B multivariate Cox analysis; C nomogram for clinical prognosis assessment (1-year, 3-year, and 5-year); D calibration curve to assess nomogram accuracy; *p < 0.05, **p < 0.01, ***p < 0.001, ns, non-significant

Subgroup analysis of prognosis-related clinical features

To clarify whether PRLs retain their prognostic ability in different subgroups and to assess their prognostic ability, in order to stratify patients according to their clinical characteristics, we conducted a stratification analysis. In Fig. 7A–N, p < 0.05, nearly all low-risk patients had significantly longer OS compared to high-risk patients with varying clinical characteristics, except in Fig. 7H (p = 0.863), which might be caused by an insufficient sample size.

Fig. 7

Kaplan–Meier curve analysis for the high- and low-risk groups stratified by clinical features, including A, B gender; C, D grade; E, F M; G, H N; I, J stage; K, L T; M, N; D age

Analysis of biological functions and immune cell infiltration levels

Aims to identify the potential causes of prognostic differences between two risk subgroups, GSEA enrichment analysis, GO, and KEGG were conducted on the DEGs within the subgroups. As shown in Supplementary Fig. 2A, ERBB, TGF-BETA, JAK-STAT, MTOR, VEGF, MAPK, and WNT signaling pathways were notably enriched in low-risk patients. According to the KEGG analysis, DEGs from the entire population mainly control the interaction between cytokine receptors, complement, and coagulation cascades, and the signaling pathway for IL-17 (Supplementary Fig. 2B). Based on the GO analysis, DEGs had a very high immune pathway enrichment (Supplementary Fig. 2C). Figure 8A depicts the intrinsic link between immune cells infiltration levels and different risk scores. There was a significant positive relationship between the high-risk subgroup and B cells, M0 macrophages, cancer-associated fibroblasts and T cells, as well as a negative relationship between the high-risk subgroup and endothelial and NK cells. Based on Fig. 8B and C, the immune and ESTIMATE scores of the high-risk patients were significantly higher than those of the low-risk patients, however, the stromal score did not show a significant difference (Fig. 8D), indicating that immune molecules were more enriched in high-risk patients. Compared to the low-risk subgroup, high-risk patients exhibited greater immune cell infiltration and immunity (Fig. 8E and F).

Fig. 8

Analysis of immune infiltration in different risk subgroups. A Immune cell bubbles of the different groups; B ESTIMATE scores of the low- and high-risk subgroups; C immune scores of the low- and high-risk subgroups; D stromal score of the low- and high-risk subgroups; E the ssGSEA scores of immune cells in different risk groups; F immune function scores of different risk groups. *p < 0.05, **p < 0.01, ***p < 0.001, ns, non-significant

Analysis of immunotherapy and drug sensitivity

Given the special status of immune checkpoints in tumor immunotherapy, both risk subgroups were compared regarding their expression levels (Fig. 9A). Compared with the low-risk subgroup, high-risk patients had higher levels of immune checkpoint expression, signifying that immunotherapy might be more effective in high-risk patients. Moreover, as shown in Fig. 9B and C, we investigated whether two different risk subgroups had different somatic mutations and found that the mutation rates of PBRM1, SETD2, and BAP1 varied significantly; in particular, the mutation rate of SETD2 varied by up to 14%, and high-risk patients had a higher frequency of genetic mutations, which might lead to their poor prognoses. To determine whether TMB score correlates with KIRC patient prognosis, a survival analysis was conducted. As shown in Fig. 9D, those with low TMB scores had a better prognosis than those with high scores. As illustrated in Fig. 9E, low PRL risk score groups with low TMB scores had a palpable survival benefit. Meanwhile, the drug sensitivity of the various groups was determined using TIDE analysis, and it seems that immunotherapy may be less effective in the high-risk subgroup as a result of lower scores on TIDE and dysfunction (Fig. 9F and H). Finally, the results of the drug sensitivity analysis showed that patients in the low-risk group were generally sensitive to most drugs, such as sorafenib, paclitaxel, sunitinib, vinblastine, and temsirolimus (Supplementary Figs. 3A–Y).

Fig. 9

Analysis of immunotherapy. A Expression of immune checkpoints in different risk groups; B waterfall plot of tumor somatic mutation in the high-risk subgroup; C waterfall plot of tumor somatic mutation in the low-risk subgroup, D Kaplan–Meier curve analysis of high- and low-TMB groups in TCGA database; E Kaplan–Meier curve analysis of TCGA KIRC data stratified by TMB and risk score; F–H TIDE, IFNG, and dysfunction scores of high- and low-risk subgroups. *p < 0.05, **p < 0.01, ***p < 0.001, ns, non-significant

Immunotherapy and PRL subgroup generation in KIRC

To explore the relationship between PRL expression levels and KIRC subtypes, KIRC patients were classified using the ConsensusClusterPlus package, and two groups were found to be the optimal clustering method (Fig. 10A). In order to identify clinical differences between the two clusters, as a first step, we compared survival curves of two groups and found significant differences in prognosis (Fig. 10B). Furthermore, using PCA and t-SNE, patients with different risks within clusters were well separated in two directions (Fig. 10C–F). Then, in order to obtain a better understanding of the clinical implications of the seven PRLs in ccRCC, we correlated clinicopathological variables with their expression levels (Fig. 10G). Risk scores varied significantly based on lymph node metastasis, tumor metastasis, cluster stratification, age, gender, AJCC T stage, grade, stage, and immune core. Additionally, tumors in cluster 2 were significantly more likely to be infiltrated by immune cells than tumors in cluster 1 (Fig. 11A). In cluster 1, immune checkpoint expression was higher than in cluster 2, which may have contributed to the lower overall survival of cluster 2 (Fig. 11B). Subsequently, according to our analysis, patients in cluster 1 displayed greater drug sensitivity, such as paclitaxel, vinblastine, sunitinib, sorafenib, and temsirolimus (Fig. 12A–Y). Based on these findings, cluster 1 and low-risk groups should have better outcomes and respond better to immunotherapy.

Fig. 10

Analysis of KIRC subtypes. A KIRC divided into two clusters; B Kaplan–Meier curve analysis of OS in clusters; C PCA analysis of the two risk groups; D t-SNE analysis of the two risk groups; E PCA analysis of the two clusters; F t-SNE analysis of the two clusters; G heatmap of correlations among clinical features, immune scores, and risk scores. *p < 0.05, **p < 0.01, ***p < 0.001, ns, non-significant

Fig. 11

Analysis of immune infiltration in the different clusters. A Heatmap of immune cells in the different clusters; B expression of immune checkpoints in the different clusters. *p < 0.05, **p < 0.01, ***p < 0.001, ns, non-significant

Fig. 12

Drug sensitivity analysis in the different clusters (A–Y)

Seven PRL expression levels in KIRC cells

Seven PRLs were determined by qRT-PCR in the 293 T and 4 ccRCC cell lines, 769-P, ACHN, 786-O, and CAKI-1. When compared to 293 T, the expression of LINC02657, ZNNT1, and LINC02747 was considerably higher in the four KIRC cell lines (Fig. 13A, B and D). In contrast, LINC02154 and AC022126.1 were significantly downregulated in KIRC cell lines (Fig. 13C and E). AP001636.3 has three subtypes: AP001636.3–1, AP001636.3–2, and AP001636.3–3, all of which had higher expression in tumor cell lines than in normal kidney cells (Fig. 13F–H). Similarly, LINC02609 has two subtypes, and the expression trends of both subtypes were consistently downregulated in tumor cells (Fig. 13I–J).

Fig. 13

qRT-PCR analysis of the seven PRLs. A LINC02657; B ZNNT1; C LINC02154; D LINC02747; E AC022126.1; F–H AP001636.3; I, J LINC02609. *p < 0.05, **p < 0.01, ***p < 0.001, ns, non-significant