Relationship of immune/stromal scores with ccRCC prognosis, identification and functional enrichment analysis of DEASs
Details of this study design are illustrated in Figure 1 as a flowchart. A total of 515 ccRCC patients with complete clinical data and transcriptome data from the TCGA database were included into the follow-up study. On the basis of ESTIMATE analysis, immune scores were distributed between -687.33 and 3339.36, stromal scores ranged from -1428.66 to 1976.04. The ccRCC patients were classified into high- and low- immune/stromal score groups on the boundary of median score. Kaplan-Meier curves displayed that the low immune-score group had longer OS than the high low immune-score group (Figure 2A). It was worth mentioning that the PFS of the low immune-score group seemed to be longer than PFS of the high immune-score group, the OS and PFS of the low stromal-score group seemed to be longer than those of the high stromal-score group, although there was no statistical significance (Figure 2B-2D). These results proved to some extent that immune/stromal score were related to the survival of ccRCC, which was worthy of further data mining.
Considering that a gene may have multiple types of AS events, UpSet plot was generated to visualize the intersecting sets of each AS type, which was illustrated in Supplementary Figure 1. We then conducted differential analysis of AS events expression in the high and low groups based on immune and stromal scores with the purpose of finding other explanations for survival difference. Specifically, the up-regulated and down-regulated AS events in high/low immune-score group and stromal-score group were showed in volcano plots (Figure 2E-2F). Considering the significant impact of intersection AS events on the prognosis of ccRCC, 400 AS events that were up-regulated in both immune and stromal groups and 461 AS events that were down-regulated in both immune and stromal groups were framed, as shown in Venn diagram, and selected as DEASs for subsequent study (Figure 2G-2H), which referred to Supplement 1 for details.
In order to assess the potential function and biological process of DEAS events, GO analysis and KEGG pathway were performed. We exhibited the top 30 results of GO analysis incorporated with the top 10 results in BP, the top 10 results in CC and the top 10 results in MF (Figure 2I). Interestingly, we found that in addition to the enrichment of immune-related response, such as " immune response−regulating cell surface receptor signaling pathway " and “Fc receptor signaling pathway”, there were many metabolism-related pathways, particularly GTPase-related, such as “GTPase binding”, “GTPase regulator activity” and “GTPase activator activity”. The KEGG pathways demonstrated that these DEASs were associated with immune-related pathways, like “Natural killer cell mediated cytotoxicity”, “T cell receptor signaling pathway”, “Chemokine signaling pathway”, “NOD−like receptor signaling pathway”, “Fc gamma R−mediated phagocytosis”, “TNF signaling pathway”, “Fc epsilon RI signaling pathway”, “B cell receptor signaling pathway”. Besides, some famous cancer-related pathways, such as “Ras signaling pathway”, “Rap1 signaling pathway”, “Apoptosis”, “ErbB signaling pathway” and “Ferroptosis”, were enriched (Figure 2J).
Establishment of two prognosis-related signatures
In order to identified survival-related DEASs, we next conducted the univariate Cox regression analysis to screen target AS events from the perspective of OS and PFS (Supplement 2, Supplement 3).In total, 455 OS-related AS events and 464 PFS-related AS events were put into LASSO analysis to obtain the the most prognostic DEASs, respectively (Figure 3A-D). Subsequently, two prognosis-related signatures related to final models were constructed by multivariate Cox regression, namely OS-related signature and PFS-related signature (Table1, Table 2). According to the results of multivariate Cox regression, the formula of OS-related signature was listed as follow: risk score= QKI*-4.617+ SAA2*-0.930+ EPB41L5*2.152+ METTL21A*-1.589+ EPC2*-3.423+ ARHGAP24*1.985+ CTNNB1*-2.365+ FOXN3*1.762+ FAHD2A*1.648+ SHC1*-1.683. Similarly, the formula of PFS-related signature was expanded as follow: risk score=QKI*-4.970+ SAA2*-0.611+ PDCD5*1.622+ C18orf32*-6.596+ SLK*-1.964+ RGS3*1.151+ METTL21A*-0.968+GLIPR1*-1.832+IDH1*-2.601+ARHGAP24*1.930+SF1*3.374. The risk scores of 507 ccRCC patients with completed OS information and 505 ccRCC patients with completed PFS information were calculated. Bounded by the medium risk score, all patients were split into high-risk group and low risk group. Kaplan-Meier curves reflected the fact that high risk group had worse OS than the low-risk group, which was the same as the result of PFS (Figure 3E, F). Additionally, ROC curves were then drawn with the aim of validating the predictive capacity of these models. The AUC value of OS-related signature for predicting 1-, 3-, and 5-year OS ranged 0.787 from 0.812, and the AUC value of PFS-related signature for predicting 1-, 3-, and 5-year OS were all greater than 0.775, which presented satisfactory performances (Figure 3G, H).
Inspired by the strong predictive power of the prediction models, we wanted to know whether the risk score had the ability to predict independently. The univariate and multivariate Cox regression analysis were performed to validate the independence of risk scores by combining relevant clinical information, such as Sex, Age, Grade, AJCC Stage, T, N, M, with risk scores (Table 3, Table 4). Surprisingly, risk score was an independent prognostic factor in both OS and PFS, independent of other clinical factors, which was consistent with the results of age.
AS-based clusters significantly corelated with prognosis and immune features
Driven by the outcomes of enrichment analysis, the unsupervised consensus analysis was conducted to assess the impacts of AS events on immune microenvironment of ccRCC. We then classified the TCGA ccRCC patients into three clusters in the help of the consensus matrix heatmap, specifically, cluster A (n = 146, 27.8%), cluster B (n = 182, 35.9%) and cluster C (n = 179, 35.3%) (Figure 4A). What’s more, Kaplan–Meier analysis of the relationship between clusters and prognosis revealed that different clusters kept different survival patterns. Obviously, cluster A was associated with the best prognosis, while cluster C tended to carry the poorest prognosis, either from the perspective of OS or PFS (Figure 4B-4C). The rare conclusion that the low immune-score group had a better prognosis further led us to wonder about the specific effect of AS events on the immune infiltration microenvironment and prognosis of renal cell carcinoma. Through GSVA analysis between clusters, we found that compared with cluster A, cluster C and cluster B were enriched more significantly in immune-related pathways, such as “Intestinal immune network for IgA production”, “Fc gamma R−mediated phagocytosis”, “T cell receptor signaling pathway”, “Chemokine signaling pathway”, “Cytokine-cytokine receptor interaction”, “Natural killer cell mediated cytotoxicity” (Figure 4D-4E). Additionally, cluster C were enriched more significantly in immune-related pathways than cluster B, such as “Intestinal immune network for IgA production” and “Natural killer cell mediated cytotoxicity”(Figure 4F).What’s more, the subsequent microscopic display of immune cell infiltration also showed that the immune infiltration level of cluster C was the highest, but that of cluster A was the lowest, whether it was immune activated cells or immunosuppressive cells, which was consistent with the previous conclusion (Figure 4G).
Identification of intersection SFs and the regulatory network of survival-related AS events and intersection SFs
A total of 52 SFs were investigated in the present study. To screen out SFs related to tumor immune microenvironment of ccRCC, differential analysis of SFs expression in the high and low groups based on immune and stromal scores were performed, whose analysis method was the same as the differential analysis of AS events. Theses heatmaps exhibited the SF expression between the high immune-score group and low immune-score group, the high stromal-score group and low stromal-score group (Figure 5A, B). Specifically, the up-regulated and down-regulated SF in high/low immune-score group and stromal-score group were showed in volcano plots (Figure 5C, D). Taking intersection SFs as the core of research direction, 3 SFs that were up-regulated in both immune and stromal groups and 1 SFs that was down-regulated in both immune and stromal groups were chose for subsequent study, as shown in Venn diagram (Figure 5E, F). Furthermore, we found the 4 intersection SFs, including YBX1, MBNL1, ELAVL4 and PCBP2, were positively correlated with immune cell infiltration, including the B cell, CD4+ T cell, CD8+ T cell, macrophage, neutrophil and dendritic cell (Supplementary Figure 2).
SFs are RNA-binding proteins that influence exon selection and splicing site selection by recognizing cis-regulatory elements in precursor mRNAs, suggesting that these survival-related AS events may be regulated by some key SFs. We collected 17 survival-related AS events according to the result of multivariate Cox regression analysis of OS and PFS. Subsequent correlation analysis between survival-related AS events and intersection SFs allowed us to finalize 16 SF-related AS events, as shown in Figure 5H. A total of 16 AS events including 9 AS events with low risk (favorable prognosis, green dots) and 7 AS events with high risk (inferior prognosis, red dots) were significantly correlated with expression levels of these 4 SFs (bule dots). The negative correlation between AS with good prognosis and SF was the majority (red line), and the positive correlation between AS with poor prognosis and SF was the majority (green line). Motivated by the regulatory network, correlation matrix was further constructed to figure out the relation between intersection AS events and intersection SFs (Figure 5I). Based on the results of the above regulatory network, correlation analysis and literature review, we selected MBNL1 from 4 hub SFs for subsequent validation in vitro.
Tumor suppressor phenotype of MBNL1 in 786O and A498 cells
Stable expression cell lines 786O and A498 were established, and the expression of MBNL1 were verified through western blot and RT-qPCR (Figure 6A). CCK8 assay showed that the proliferation capability of MBNL1 overexpressed cells was declined (Figure 6B). Meanwhile, the migration and invasion of MBNL1 overexpressed cells were found to be reduced in Transwell assay (Figure 6C). In the cell cycle and apoptosis experiments, it was found that the overexpressed MBNL1 cells were arrested in the G2/M phase, and the proportion of apoptosis was increased in the cells of transient transfection with MBNL1 expression (Figure 6D-6E).