The transcriptome and survival data of the ccRCC patients in this study came from The Cancer Genome Atlas portal website (TCGA). Also, the AS data of TCGA came from SpliceSeq, and the SF expression data were obtained from the SpliceAid 2 database (www.introni.it/spliceaid.html). All analyses strictly followed TCGA’s published guidelines, and the detailed analysis flowchart can be found in Figure 1 .

When setting the PSI value above 0.75 as the point for filtration, samples were partitioned. Using the UpSetR software package, the UpSet plot was drawn and seven subtypes of AS events were found. We named AS by splicing types, ID numbers in splicing sequences, and corresponding parental gene names. There was a case that C4orf19 was the corresponding parent gene name, 69001 was the ID number in SpliceSeq, and AT was the splicing type in “C4orf19|69001|AT”.

When we detected the PSI standard deviation less than 0.01, the data of AS events were deleted. The connection between the overall survival (OS) and AS events was found in the univariate Cox regression analysis ( Additional file 1: Table 1 ), which was exhibited in the UpSet map and volcano map. In addition, each bubble chart of the seven subtypes summarized the 20 most important AS events.

Firstly, candidate models for each splicing pattern were determined by least absolute shrinkage and selection operator (LASSO) regression analysis, in which way we could also avoid model overfitting. Next, multivariate Cox regression analysis was applied to screen prognostic predictors from the identified AS events. Because the pattern of AS events in post-transcriptional modification was independent of each subtype, the AS events identified in each of the splicing subtypes described above were integrated and then another prognostic feature was generated. Afterward, risk scores were calculated according to the formula: risk score = βAS event1 × PSIAS event1 + ⋯ + βAS eventn × PSIAS eventn. The specific formulas for each prognostic signature can be found in Additional file 1: Table 2 . Consequently, the low-risk group and the high-risk group were born by the calculated median risk scores. The “survival” R package was employed to analyze K–M survival curves. The predictive value of this prognostic signature was validated by using time-dependent receiver operating characteristic (ROC) curves. Then, univariate and multivariate Cox regression analyses were exploited to ascertain whether this signature could be used as an independent prognostic factor. In addition, stratified survival analysis further verified whether prognostic performance in patients was independent of clinical data including age; sex; pathological grade; T, N, and M categories; and tumor stage. Then, we calculate the AUC from the ROC curve to systematically measure the value of the accuracy of the model for 1-, 2-, and 3-year OS. Finally, to accurately calculate the OS of ccRCC patients, we established prognostic nomograms to obtain the survival probabilities of 1, 2, and 3 years. Then, there was a calibration curve showing the prognostic value of the AS-constructed nomogram. It should be noted that the model was highly predictive when the calibration curve was close to 45°.

Information on immune infiltrates such as B cells from each specimen was downloaded from TIMER. The ssGSEA algorithm of the R package “GSEAbase” was performed to elucidate the enrichment of two different risk subgroups in 29 gene sets related to immune function. Subsequently, we calculated the purity of the tumor and the degree of cell invasion (stromal and immune cells) using the R package “ESTIMATE” to validate the significantly different characteristics of the TIME between the low-risk and high-risk groups. The proportion of 22 immune cell types in the tumor sample was recognized by assessing the relative subsections of RNA transcripts from CIBERSORT.

According to existing research, the expression level of key genes associated with immune checkpoint blockade might have a close relationship with the clinical results of ICB treatment (14, 15). Six key genes (PD‐L1, IDO1, PD‐L2, PD‐1, CTLA‐4, and TIM‐3) of immune checkpoint blockade therapy in ccRCC (16, 17) were obtained. Afterward, to investigate the potential role of risk score in immune checkpoint blockade therapy of ccRCC, AS-based prognostic characteristics were significantly related to the expression levels of four key genes for immune checkpoint blockade. At last, the expression levels of 47 immune checkpoint genes (i.e., CTLA4, BTLA, etc.) were compared in low-risk and high-risk patients.

A total of 404 SFs derived from a previous study (18) are exhibited in Additional file 1: Table 3 , and the RNA-seq profiles of SFs can be found in the TCGA database. In addition, we conducted a Spearman correlation analysis to assess the connection between SFs and survival-related AS events. p <0.001 and correlation coefficient >0.6 were the cutoff values. In the end, Cytoscape (version 3.9.1) was applied to build an underlying SF-AS regulatory network.

In this study, for comparisons between two different groups, we used the Wilcoxon test, and for comparisons between more than two groups, we used the Kruskal–Wallis test. OS was the time between diagnosis and death. The K–M log-rank test was employed to plot the survival curse. Moreover, the Pearson correlation test was applied to explore the correlation between risk score, clinical variables, and degree of immune cell infiltration and immune checkpoint. When the result of the CIBERSORT algorithm p ≥ 0.05, further study was abandoned. Then, in order to verify the independent prognostic prediction abilities of risk signatures, univariate and multivariate analyses were carried out by the Cox regression model. For 1-, 2-, and 3-year OS, we used ROC curves to evaluate their prognostic value. p <0.05 was regarded as statistically significant. All statistical analyses were performed using R software in version 4.1.2.

Five hundred thirty-seven patients with ccRCC were obtained from the TCGA database, and 11 patients without complete information were rejected. Therefore, 526 patients in total were included. Table 2 presents the basic clinical data of all ccRCC patients. In addition, by using the UpSet plot ( Figure 2A ), we analyzed the AS event profiles comprehensively and displayed gene intersections among the seven subtypes of AS events. It could be seen that ES was the most frequent splicing pattern, while ME was the least frequent.

Univariate Cox regression analysis showed that 9,347 AS events were significantly associated with survival (p < 0.05). Furthermore, a detailed record of the data can be found in Additional file 1: Table 1 . In Figure 2B , the gene interactions among the seven types of survival-related AS events are shown. Moreover, ES was still the main splicing pattern. On the other hand, the volcano map was designed to present the distribution of AS events ( Figure 3A ), and the top 20 AS events with significant survival correlation from seven subtypes were summarized by using the bubble graphs ( Figures 3B–H ).

In this study, the prognostic abilities of the survival-related AS events found in the previous step were evaluated by using the stepwise LASSO algorithm and multivariate Cox regression analysis. Moreover, the LASSO regression analysis results of ALL AS events and seven AS event subtypes are exhibited in Figures 4A, B , 5A–G , 6A–G . Next, the best survival-related AS events, which were determined by multivariate Cox analysis, were performed to build eight AS prognostic signatures, namely, AA, AD, AP, AT, ES, ME, RI, and ALL. The formulas for each prognostic signature are detailed in Additional file 1: Table 2 . Using the median risk score as a standard for further study, ccRCC patients were ranked into low- and high-risk groups. The distribution of eight different AS events (AA, AP, AT, AD, ME, RI, ES, and ALL) and their PSI values in the two subgroups and patients was exhibited in the heatmap ( Figures 4C , 7A, D , 8A, D , 9A, D , 10A ). In the same way, the distribution of risk score ( Figures 4D , 7B, E , 8B, E , 9B, E , 10B ) and the dot plot of survival status ( Figures 4E , 7C, F , 8C, F , 9C, F , 10C ) indicated a lower overall survival in the higher-risk patients. Furthermore, the Kaplan–Meier curve also confirmed that in the low-risk subgroup, patients had a significantly better prognosis than those in the high-risk subgroup ( Figures 4A , 11A, C, E, G , 12A, C, E, G ; all P < 0.05). The results showed that the areas under the risk score curves of 1-, 2- and 3-year survival were all greater than 0.70, indicating that the established prognostic signature had highly sensitive and specific survival prediction ability ( Figures 4G , 11B, D, F, G , 12B, D, F, G ). Moreover, the risk score might become an independent prognostic signature of the ccRCC (univariate Cox model in Figures 4H and 13A, C, E, G, I, K, M and multivariate Cox regression analysis in Figures 4I and 13B, D, F, H, J, L, N ).

According to the difference in the risk score in different subtypes of clinical variables, clinical significance was explored. With the progression of tumor grade (most p < 0.05, Figure 14A ); clinicopathological stage (most p < 0.05, Figure 14B ); and T, M, and N stages (most p < 0.05, Figures 14C–E ), the risk score significantly rose, suggesting that prognostic risk score had a positive correlation with tumor progression. Next, the prognostic nomogram established for forecasting the prognosis of ccRCC patients is exhibited in Figure 14F . It was well known that there was a great prognostic capability of 1-, 2-, and 3-year OS in the nomogram plot when the calibration curve was close to the diagonal ( Figures 14G–I ).

In order to further investigate the possibility of using risk score as an immune indicator, we performed correlation analyses between risk score and immune score (from the ESTIMATE algorithm), ssGSEA characteristics, and TIC subtypes and levels (from the CIBERSORT method). The high-risk patients achieved a higher immune score and ESTIMATE score and lower tumor purity ( Figures 15A–C ), which suggested higher immune infiltration. However, there was no significant difference in stromal score ( Figure 16 ). Then, in Figures 15D, E , immune-related signatures were shown to differ between the two subgroups, where immunological scores corresponding to immune-related signatures were exhibited for each patient in the low-/high-risk group. The results revealed that the infiltration of immune cells such as CD8+ T cells, macrophages, T helper cells, Tfh, Th1 cells, Th2 cells, and TIL and the immune signatures such as parainflammation, T-cell co-inhibition, T-cell co-stimulation, checkpoint, inflammation-promoting, and cytolytic activity were significantly increased with increased risk score ( Figure 15F ). On the contrary, iDCs, mast cells, and type IIIFN response were significantly decreased with increased risk score ( Figure 15F ). The CIBERSORT algorithm results showed that the proportion of CD8+ T cells, activated CD4 memory T cells, follicular helper T cells, Tregs, and M0 macrophages was positively associated with risk score, and the proportion of naive B cells, memory B cells, M1 macrophages, M2 macrophages, resting dendritic cells, and resting mast cells was negatively associated with risk score ( Figure 15G ). In conclusion, the ALL prognostic signature could be a kind of new method to clarify the ccRCC immunoregulatory network.

With the increasing attention paid to ICB therapy in clinical work, immune checkpoint inhibitors have greatly changed the clinical decision-making of cancer oncology (19, 20). We screened out six key immune checkpoint inhibitor genes (PDCD1, CD274, PDCD1LG2, CTLA‐4, HAVCR2, and IDO1) (21, 22) for further analysis. Then, to uncover the potential role of risk signature in ICB therapy for ccRCC, we comprehensively analyzed the association between the ALL prognostic signature and ICB key targets ( Figure 17A ). The results showed that the ALL prognostic signature had a significant positive association with PDCD1 (r = 0.3; p = 6.1e−12; Figure 17B ) and CTLA4 (r = 0.33; p = 8.8e−15; Figure 17D ) and a significant negative association with HAVCR2 (r = −0.13; p = 0.0035; Figure 17C ) and CD274 (r = −0.13; p = 0.0025; Figure 17E ). Furthermore, 36 of the 47 (i.e., HHLA2, CD44, etc.) ICB key gene expression levels between the low- and high-risk groups were significantly dysregulated in the further correlation analysis ( Figure 17F ). These results suggested that the level of the ALL prognostic signature does affect the expression changes of ICB key genes, which could be a valuable factor.

In this study, we found only one prognostic AS-related gene, C4ORF19, whose expression level was significantly downregulated. According to the TCGA database, the expression level of C4ORF19 in normal adjacent tissues was higher than that in tumor tissues ( Figure 18A ). The IHC experiment showed that the expression level of C4ORF19 in ccRCC tissue was significantly lower than that in normal tissue, and the experimental results of qPCR also confirmed this point (p = 0.0115) ( Figure 19 ). It could be clearly seen that the expression levels of C4ORF19 in different tumor grades ( Figure 18C , almost p < 0.05), different pathological stages ( Figure 18D , almost p < 0.05), T state, M state, and gender ( Figures 18E–G , almost p < 0.05) had significant statistical significance. In order to further assess the prognostic value of C4ORF19 in ccRCC, K–M analyses were performed between patients with low and high expression of C4ORF19. A higher C4ORF19 expression level significantly correlated with a longer overall survival time, as illustrated in Figure 18B (p < 0.001). Moreover, in 28 of 47 immune check blockade-associated genes (i.e., PDCD1, CTLA4, etc.), there were significant dysregulations in the expression levels between the low C4ORF19 group and high C4ORF19 group in different subgroups ( Figure 18H ). Then, a possible role for C4ORF19 in ICB treatment of ccRCC was explored by analyzing the association between C4ORF19 and ICB key targets adjusted for tumor purity using TIMER. The TIMER results exhibited that C4ORF19 had a significant positive correlation with CD274 (r = 0.361; p = 1.21e−15) and HAVCR2 (r = 0.137; p = 3.15e−03) and a significant negative correlation with PDCD1 (r = −0.129; p = 5.46e−03) and CTLA4 (r = −0.095; p = 4.11e−02; Figure 18I ), suggesting that C4ORF19 may play a vital role in the ICB treatment of ccRCC.

Firstly, we classified ccRCC patients into high/low C4ORF19 groups for further study according to the median C4ORF19 expression level. The ESTIMATE results showed significantly higher stromal and immune scores in the low C4ORF19 group than in the high C4ORF19 group, suggesting more infiltration of stromal and immune cells and lower tumor purity in the low C4ORF19 group ( Figures 20A–D ). Moreover, the relationship between the gene copy number of the different mutation types and main immune cells is exhibited in Figure 20E . Afterward, a positive correlation was found between C4ORF19 expression level and B-cell infiltration, while a negative correlation was found between C4ORF19 expression level and CD8+ T-cell infiltration. There was no significant difference in the expression level of C4ORF19 when CD4+ T cells, macrophages, and neutrophils were infiltrated ( Figure 20F ). The consequences of ssGSEA presented that the infiltration fraction of aDCs, CD8+ T cells, DCs, macrophages, pDCs, Th1 cells, Th2 cells, NK cells, parainflammation, T helper cells, Tfh, TIL, APC co-stimulation, checkpoint, T-cell co-stimulation, CCR, cytolytic activity, inflammation-promoting, and IFN-response type-I were significantly increased when the C4ORF19 expression level was declining ( Figure 20G ). The consequences of the CIBERSORT analysis of the TCGA cohort presented that the proportions of plasma cells, Tregs, activated memory CD4 T cells, and M0 macrophages were significantly higher and the proportions of monocytes and resting dendritic cells were significantly lower in patients with low C4ORF19 expression ( Figure 20H ).

The upregulated and downregulated genes were the results of the correlation analysis with the corresponding gene expression levels in tumor samples ( Additional file 1: Table 4 ). Finally, to better explain the underlying mechanisms of AS regulation, we used 351 upregulated AS events (yellow diamond), 88 downregulated AS events (green triangle), and 31 SFs (blue hexagon; Figure 21 ) to establish the correlation network between SF expression level and PSI value of prognostic AS events. In the regulation network, the most important four nodes ( Additional file 1: Table 4 ) consisting of two upregulated AS events (METTL3|26596|RI and FADS3|16305|RI) and two SFs (DDX39B and LUC7L) were identified. As a result, these SFs had great potential to further mediate the occurrence and development of tumors in ccRCC as key regulatory factors involved in abnormal AS regulation.