Alternative splicing and alternative polyadenylation define tumor immune microenvironment and pharmacogenomic landscape in clear cell renal carcinoma.

Two major posttranscriptional mechanisms-alternative splicing (AS) and alternative polyadenylation (APA)-have attracted much attention in cancer research. Nevertheless, their roles in clear cell renal carcinoma (ccRCC) are still ill defined. Herein, this study was conducted to uncover the implications of AS and APA events in ccRCC progression. Through consensus molecular clustering analysis, two AS or APA RNA processing phenotypes were separately constructed with distinct prognosis, tumor-infiltrating immune cells, responses to immunotherapy, and chemotherapy. The AS or APA score was constructed to quantify AS or APA RNA processing patterns of individual ccRCCs with principal-component analysis. Both high AS and APA scores were characterized by undesirable survival outcomes, relatively high response to immunotherapy, and low sensitivity to targeted drugs, such as sorafenib and pazopanib. Moreover, several small molecular compounds were predicted for patients with a high AS or APA score. There was a positive correlation between AS and APA scores. Their interplay contributed to poor prognosis and reshaped the tumor immune microenvironment. Collectively, this study is the first to comprehensively analyze two major posttranscriptional events in ccRCC. Our findings uncovered the potential functions of AS and APA events and identified their therapeutic potential in immunotherapy and targeted therapy.

Introduction

Renal cell carcinoma (RCC) represents the most frequent malignancy that affects the adult kidney globally. It mainly includes three histological subtypes: clear cell RCC (ccRCC), papillary RCC (pRCC), and chromophobe RCC (chRCC). The most common subtype of RCC is ccRCC, which accounts for ∼75% of cases. The behavior of ccRCC is to "clear the cytoplasm" due to its capacity to accumulate glycogen and lipids. This malignancy is mostly asymptomatic in the early stage and diagnosed occasionally by imaging, with favorable clinical outcomes. Oppositely, in the late stage, the mortality of advanced patients is quite high owing to low sensitivity to radiotherapy and chemotherapy. Although advances have been made in tumor diagnosis and therapy, surgery is the first-line curative therapeutic strategy against ccRCC. Targeted drugs like axitinib, sunitinib, and sorafenib have been applied as the standard of care. Despite the early response in selected subjects, most of them ultimately develop resistance to chemotherapy. Thus, improving the understanding of the molecular mechanisms of ccRCC may facilitate the development of treatment options.

RNA processing, including alternative splicing (AS) and alternative polyadenylation (APA) of pre-mRNAs, is a highly specialized mechanism that allows organisms to enhance transcriptome and proteome diversity. Over 90% of transcripts undergo alternative RNA processing, which are estimated to generate approximately 100,000 different proteins from 20,000 human genes. AS causes skipping/inclusion of exons or retention of introns via differential selection of splice sites in pre-mRNAs. Eukaryotic 3′ end formation is a key step in mRNA maturation, involving requisite cleavage and polyadenylation events downstream of the polyadenylation signaling. APA leads to distinct 3′ ends through regulated processing of the 3′ end, including transcription termination, cleavage, and polyadenylation. In almost 70% of human genes, cleavage and polyadenylation occur in multiple locations through APA, which is a major post-transcriptional mechanism for gene regulation. Mounting evidence highlights the critical roles of AS and APA events in human cancers, such as driving oncogenic gene expression, chemotherapy resistance, and tumor microenvironment (TME). Furthermore, global AS and APA events possess predictive potential as prognostic biomarkers of human cancers. Nevertheless, much remains unknown regarding the regulatory mechanisms and the interactions of these two processes in ccRCC. Herein, this study comprehensively recognized the regulation mechanisms of AS and APA RNA processing patterns in survival outcomes, immune response, TME reshaping, and targeted therapeutic sensitivity in ccRCC. AS and APA scoring systems were developed to quantify the RNA processing patterns in individual tumors. Our findings could enhance the understanding of TME immune regulation and assist in developing more effective therapeutic strategies.

Results

Establishment of AS RNA processing patterns with distinct prognosis, TME, activation of signaling pathways, and chemosensitivity

To further understand the AS events involved in tumorigenesis, percent spliced in index (PSI) values for AS events on ccRCC samples were obtained from The Cancer Genome Atlas (TCGA) SpliceSeq database. We applied a non-negative matrix factorization (NMF) algorithm to cluster these ccRCC patients based on the top 1,000 AS events with the most variance. As a result, 520 ccRCC patients were classified into two AS RNA processing patterns, named AS cluster 1 (n = 264) and 2 (n = 256) (Figures S1A and S1B). The t-distributed stochastic neighbor embedding (t-SNE) results confirmed the accuracy of this classification (Figure S1C). The AS clustering patterns remained based on the top 5,000 AS events with the most variance (Figures S2A–S2C) or median absolute deviation (MAD) (Figures S3A–S3C). We conducted survival analysis on ccRCC patients to observe the prognostic difference between AS cluster 1 and 2. In Figure 1A, AS cluster 2 exhibited a significant survival advantage compared with AS cluster 1 (p = 1.203 × 10−5). The ESTIMATE method was presented to estimate the overall infiltration levels of immune and stromal cells in ccRCC samples. Notably, AS cluster 1 had elevated immune scores (p = 0.026; Figure 1B) and reduced stromal scores (p = 0.0052; Figure 1C). The differences in the TME cell infiltration were compared between clusters. As depicted in Figure 1D, antitumor lymphocyte cells, such as effector memory CD4+ T cells, central memory CD8+ T cells, natural killer T cells, and natural killer cells, exhibited higher infiltration levels in AS cluster 2, indicative of favorable prognosis. Moreover, there were distinct differences in the mRNA expression of immune checkpoints (Figure 1E) between clusters. Notably, ADORA2A, BTNL2, CD160, CD244, TNFRSF18, LAG3, and PDCD1 mRNAs were significantly up-regulated in AS cluster 1 whereas CD276 and NRP1 mRNAs were markedly up-regulated in AS cluster 2. In Figure 1F, the human leukocyte antigen (HLA) family HLA-3, HLA-DRB5, HLA-DRA, and HLA-DPA1 had significantly increased mRNA expression in AS cluster 2. These suggested that AS RNA was significantly associated with the degree and types of infiltrating immune cells in ccRCC. To further characterize the mechanisms underlying the two RNA processing patterns, we presented gene set variation analysis (GSVA). In Figure 1G, most of the pathways, such as carcinogenic pathways (RCC, pathway in cancer, and Wnt pathway) and metabolism pathways (such as glutathione metabolism, sphingolipid metabolism, and starch and sucrose metabolism), were significantly activated in AS cluster 2. We then investigated whether there were differences in the sensitivity to commonly used targeted drugs between clusters (Figures 1H–1K). The data showed that AS cluster 2 had the lower half-maximal inhibitory concentration (IC50) values of sorafenib (p = 4 × 10−22) and pazopanib (p = 7.22 × 10−17), indicating that patients in AS cluster 2 were more likely to benefit from sorafenib and pazopanib.

Figure 1

Construction of two AS RNA processing patterns characterized by distinct prognosis, TME, activation of signaling pathways and responses to targeted drugs

(A) Kaplan-Meier curves of OS for ccRCC patients in AS cluster 1 and 2 (log rank test). (B and C) The differences in immune score and stromal score between AS clusters by applying the ESTIMATE algorithm. (D) Heatmap of the infiltration levels in 28 immune cells between AS clusters via the ssGSEA algorithm. (E) Heatmap for visualizing the mRNA expression of immune checkpoints in AS cluster 1 and 2. (F) Heatmap showing the mRNA expression of HLA genes in AS cluster 1 and 2. (G) GSVA enrichment analysis for the differences in activation of KEGG pathways between AS cluster 1 and 2. (H–K) The differences in estimated IC50 values of sorafenib, sunitinib, axitinib, and pazopanib between AS clusters. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001; ∗∗∗∗p < 0.0001.

Quantification of AS score as a robust and independent prognostic indicator of ccRCC

Given the heterogeneity between AS RNA processing patterns, we identified 162 differential AS events between patterns, which were finally determined with the Boruta algorithm (Table S1). To further verify the AS post-transcriptional mechanism, we presented consensus molecular clustering analysis based on the expression profiling of the above AS genes. Consistent with the classification of AS RNA processing patterns, the ccRCC patients were clustered into two AS genomic phenotypes, named AS gene cluster 1 and 2 (Figures S4A–S4C). We observed that patients in AS gene cluster 2 had a more undesirable prognosis than those in AS gene cluster 1 (p = 1.129 × 10−5; Figure 2A). Considering the individual heterogeneity and complexity of AS post-transcriptional modification, on the basis of the above AS genes, this study developed an AS scoring system to quantify the AS RNA processing patterns of individual ccRCC patients, named the AS score. Figure 2B visualized the attribute variations of individual patients. We further evaluated the implication of the AS score in predicting patients' outcomes. With the median value, patients were stratified into high and low AS score groups. Survival analysis revealed that a high AS score was indicative of poorer overall survival (OS) (p = 1.279 × 10−5; Figure 2C), disease-free survival (DFS) (p = 4.622 × 10−1; Figure 2D), disease-specific survival (DSS) (p = 1.885 × 10−3; Figure 2E), and progression-free interval (PFI) (p = 1.31 × 10−1; Figure 2F) in comparison with a low AS score. However, no significant differences in AS were found among different clinicopathological characteristics (Figure 2G). Receiver operating characteristic (ROC) curves were conducted to better illustrate the clinical utility of AS score. In Figure 2H, we found that AS score had a prominent advantage in predicting long-term survival time. We evaluated whether the AS score could act as an independent prognostic biomarker for ccRCC. Multivariate Cox regression model analysis confirmed AS score, stage, and grade as robust and independent prognostic indicators for assessing a patient's OS time (Figure 2I). A nomogram that included AS score, stage, and grade was established for predicting 1-, 3-, and 5-year OS probabilities (Figure 2J). The excellent predictive efficacy of this nomogram was confirmed by ROC (Figure 2K), calibration (Figures S5A–S5C), and decisive curves (Figures S5D–S5F).

Figure 2

Development of AS score as a robust and independent prognostic indicator of ccRCC

(A) OS analysis between AS gene cluster 1 and 2 across ccRCC patients (log rank test). (B) Alluvial diagram visualizing the changes of AS RNA processing patterns, AS genomic phenotypes, AS score, and survival status. (C–F) The differences in OS, DFS, DSS, and PFI between high and low AS score groups (log rank tests). (G) The distribution of AS score in groups of different clinicopathological characteristics, including age <65 and ≥65, female and male, grade 1–4, and stage I–IV. p values were determined with Wilcoxon or Kruskal-Wallis test. (H) Time-independent ROC curves of age, stage, grade, gender, and AS score across ccRCC patients. (I) Multivariate Cox regression analysis for assessing the independent prognostic indicators of OS and DFS. (J) Establishment of a nomogram that included stage, grade, and AS score in predicting 1-, 3-, and 5-year OS. (K) ROC curves under 1-, 3-, and 5-year OS based on AS score.

AS score in the role of survival outcomes and TME immune regulation

To uncover the role of AS score in the TME immune regulation, we quantified immune score and stromal score in ccRCC samples. Our results showed that a high AS score was characterized by increased immune score (p = 0.015; Figure 3A) and reduced stromal score (p = 0.012; Figure 3B). We observed that the mRNA expression of HLA genes was significantly down-regulated in high AS score samples (Figure 3C). We then investigated the specific differences of tumor-infiltrating immune cells between high and low AS score groups. Antitumor lymphocyte cells, such as activated CD4+ T cells, effector memory CD4+ T cells, central memory CD8+ T cells, effector memory CD4+ T cells, natural killer T cells, and natural killer cells, were significantly activated in the low AS group (Figure 3D), indicative of favorable survival outcomes. Increasing evidence has revealed the roles of AS events in change of tumor-infiltrating immune cell levels. For example, MALT1 controls signaling and activation of CD4+ T cells. Switch of HLA-G AS causes loss of HLA-G1 expression and sensitivity to natural killer lysis. In Figure 3E, immune checkpoints CD80, BTLA, CD40LG, CD244, TNFSF9, CD70, TMIGD2, PDCD1, BTNL2, LAG3, ADORA2A, CD160, and TNFRSF18 displayed significantly increased mRNA expression in high AS score samples, whereas others were significantly decreased in low AS score samples. For example, experimental evidence has confirmed that CD40 and PDCD1 expression is controlled by posttranscriptional and posttranslational regulation through AS. To uncover the potential mechanisms involved in the AS score, we carried out gene set enrichment analysis (GSEA) across ccRCC samples. As a result, adherens junction, focal adhesion, and RCC were significantly activated in the low AS score group (Figure 3F). We also examined the correlation between the known signatures and AS score. Our results showed that AS score was negatively correlated to epithelial-mesenchymal transition (EMT), Wnt, fibroblast growth factor receptor 3 (FGFR3), and angiogenesis (Figure 3G), indicating that a low AS score could be linked to stromal activation. Growing evidence shows that the activation of the above pathways is in relation to immunotherapy resistance. For example, WNT signaling activation mediates T cell exclusion and tumor cell escape from the immune system. EMT-mediated programmed death ligand-1 (PD-L1) accumulation on cancer stem cells facilitates immune evasion. Antiangiogenesis increases the infiltration of immune effector cells into tumors and converts the intrinsically immunosuppressive TME to an immunosupportive one. Thus, AS score might be negatively associated with immunotherapy resistance. Also, we found that AS score displayed the negative correlation to DNA replication, nucleotide excision repair, and mismatch repair.

Figure 3

Characterization of AS score in the role of survival outcomes and TME immune regulation

(A and B) Differences in immune score and stromal score between high and low AS score groups. (C) The mRNA expression of HLA genes in two groups. (D) The abundance of each tumor-infiltrating immune cell in two groups. (E) The mRNA expression of immune checkpoints in two groups. (F) GSEA for the activated signaling pathways in the low AS score group. (G) Correlation between the known signatures and AS score. Ns: not significant; ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001.

Characteristics of APA score in somatic mutation and prediction of candidate drugs against ccRCC

We conducted a comprehensive analysis to evaluate the somatic mutation difference between high and low AS score groups. A higher frequency of somatic mutation was observed in low AS score samples (Figures 4A and 4B). Moreover, a low AS score displayed markedly increased microsatellite instability (MSI) than did a high AS score (p = 5.9 × 10−5; Figure 4C). The mRNA expression-based stemness index (mRNAsi) was quantified for reflecting cancer stemness. In Figure 4D, there was a markedly elevated mRNAsi in low AS score samples (p = 2.89 × 10−42). We also investigated tumor mutation burden (TMB) difference between high and low AS score groups. A higher TMB was found in a low AS score (Figure 4E). The ccRCC patients were stratified into four groups, according to AS score, and TMB and survival analysis revealed that patients with a low AS score and low TMB exhibited a prominent survival advantage (p = 1.093 × 10−4; Figure 4F). Thus, AS score could be closely linked to somatic mutation. There is growing evidence showing that AS is tightly regulated and tightly interplays with genetic and epigenetic machinery. Although somatic mutations regulate AS patterns, AS also controls genomic stability, chromatin organization, and transcriptome. The response to common targeted drugs sorafenib, sunitinib, axitinib, and pazopanib was compared between groups (Figures 4G–4J). We observed that a low AS score was more likely to benefit from sorafenib (p = 8.33 × 10−22) and pazopanib (p = 3.02 × 10−17). Two main methods were applied for screening candidate compounds with higher drug sensitivity in high AS score patients utilizing Cancer Therapeutics Response Portal (CTRP)- and Profiling Relative Inhibition Simultaneously in Mixtures (PRISM)-derived drug response data. We first conducted differential drug response analysis between high and low AS score groups to identify compounds with lower estimated area under the curve (AUC) values in the high APA score group. Then, Spearman correlation analysis between AUC value and AS score was carried out to screen compounds with negative correlation coefficients. As a result, we identified eight CTRP-derived compounds (including ouabain [r = −0.42], SR-II-138A [r = −0.45], pevonedistat [r = −0.51], daporinad [r = −0.46], CR-1-31B [r = −0.49], narciclasine [r = −0.456], PHA-793887 [r = −0.47], and STF-31 [r = −0.45]; Figures 4K and 4L) and eight PRISM-derived compounds (including halcinonide [r = −0.38], VE-822 [r = −0.38], RKI-1447 [r = −0.67], ABT-702 [r = −0.47], BIBU-1361 [r = −0.42], flumethasone [r = −0.36], YM-155 [r = −0.49], and GW-583340 [r = −0.62]; Figures 4M and 4N). The above compounds displayed lower estimated AUC values in the high AS score group and were negatively associated with AS score. The therapeutic potential of these compounds in ccRCC was further investigated by Connectivity Map (CMap) analysis. Narciclasine, GW-583340, and BIBU-1361, in which gene expression increased in ccRCC tissues but decreased through treatment with certain compounds, had potential therapeutic effects in ccRCCs with a high AS score (Figure 4O).

Figure 4

Characteristics of APA score in somatic mutation and prediction of candidate drugs against ccRCC

(A and B) Landscape of somatic mutation in high and low AS score samples. The number on the right represents the mutation frequency in each gene. The right bar plot represents the proportion of each mutation type. (C–E) The differences in MSI, mRNAsi, and TMB in high and low AS score groups. (F) Survival analysis among four groups stratified by AS score and TMB score (log rank test). (G–J) Estimated IC50 values of sorafenib, sunitinib, axitinib, and pazopanib in high and low AS score samples. (K and L) Spearman's correlation analysis and differential drug response analysis of eight CTRP-derived compounds in high AS score samples. Lower AUC value indicates higher drug sensitivity. (M and N) Spearman's correlation analysis and differential drug response analysis of eight PRISM-derived compounds in high AS score samples. (O) Validation of the above compounds in the CMap database. Lower CMap score indicates higher therapeutic potential. ∗∗∗p < 0.001.

Construction of APA RNA processing patterns with distinct survival, activation of signaling pathways, TME, and chemosensitivity

Growing evidence suggests that an APA event has a strong efficacy in predicting clinical outcomes of human cancers. Herein, the distal polyA site usage index (PDUI) value was used to indicate the frequency of APA events. By applying the NMF algorithm, ccRCC patients were classified into two APA RNA processing patterns based on the top 1,000 APA events with the most variance (Figures S6A–S6C). We named the two patterns APA cluster 1 and 2. The APA clustering patterns remained based on the top 5,000 AS events with the most variance (Figures S7A–S7C) or MAD (Figures S8A–S8C). Survival analysis of the two patterns revealed a prominent survival advantage in APA cluster 1 (p = 3.643 × 10−4; Figure 5A). APA is a crucial event for gene regulation and is in relation to cancer progression. Thus, we determined the activation of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in two APA patterns through GSVA enrichment analysis. In Figure 5B, carcinogenic activation pathways (such as Notch signaling pathway, transforming growth factor [TGF]-β signaling pathway, mammalian target of rapamycin [mTOR] signaling pathway, and RCC and metabolism pathways [such as aminoacyl tRNA biosynthesis, lysine degradation, selenoamino acid metabolism, terpenoid backbone degradation, and citrate cycle tricarboxylic acid (TCA) cycle]) were distinctly activated in APA cluster 1. The ESTIMATE algorithm was employed for estimating the overall infiltration of stromal and immune cells in ccRCC. We found that APA cluster 2 was characterized by increased immune score (p = 0.034; Figure 5C) and decreased stromal score (p = 0.00042; Figure 5D). We then examined the expression of HLA genes and immune checkpoints in APA cluster 1 and 2. Most of HLA genes (Figure 5E) and immune checkpoints (Figure 5F) displayed higher expression in APA cluster 1 than 2. By single-sample GSEA (ssGSEA) algorithm, tumor-infiltrating immune cell components were determined. Activated CD4 and CD8 T cells, central memory CD8 T cells, and CD56bright natural killer cells displayed significantly higher infiltration levels in APA cluster 1 compared with 2 (Figure 5G). The above data indicated that APA was significantly associated with TME-infiltrating cell types and degree in ccRCC. Lower IC50 values of sorafenib (p = 0.009; Figure 5H), axitinib (p = 5.53 × 10−8; Figure 5I), and pazopanib (p = 5.65 × 10−6; Figure 5J) were found in APA cluster 1, indicating that patients in APA cluster 1 were more sensitive to the above targeted drugs. Meanwhile, APA cluster 2 displayed higher sensitivity to sunitinib (p = 1.67 × 10−5; Figure 5K).

Figure 5

Construction of APA RNA processing patterns with distinct survival, activation of signaling pathways, TME, and chemosensitivity

(A) Survival analysis of patients in APA cluster 1 and 2 (log rank test). (B) GSVA enrichment analysis for the differences in activation of KEGG pathways between APA clusters. (C and D) The differences in immune score and stromal score between APA clusters. (E and F) Heatmap for visualizing the mRNA expression of HLA and immune checkpoints in APA clusters. (G) Heatmap of the infiltration levels in 28 immune cells between APA clusters. (H–K) The differences in estimated IC50 values of sorafenib, sunitinib, axitinib, and pazopanib between APA clusters. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001; ∗∗∗∗p < 0.0001.

Development of APA scoring system and characteristics of clinical traits in APA-related phenotypes

To investigate the underlying biological behaviors of each APA molecular subtype, this study determined 161 APA-related genes through the linear models for microarray data (limma) package and Boruta algorithm (Table S2). Functional annotation analysis of the APA-related genes was performed via the clusterProfiler package. Gene ontology (GO) enrichment results revealed that the above genes displayed distinct enrichment in immunity-related biological processes and pathways, such as interleukin-27- and interleukin-35-mediated signaling pathways (Figures 6A–6D), confirming that APA events exerted a non-negligible role in regulating cancer immunity. To further observe this regulation mechanism, unsupervised clustering analysis was employed to cluster patients into two genomic subtypes according to the identified 161 APA-related genes (Figures S9A and S9B), which was consistent with APA RNA processing patterns. APA gene cluster 1 exhibited a prominent survival advantage compared with APA gene cluster 2 (p = 5.811 × 10−4; Figure 6E). To accurately predict APA patterns in individual patients, this study developed an APA scoring system for ccRCC, named the APA score. The alluvial diagram was depicted for visualizing the attribute change in each patient (Figure 6F). Figure 6G depicted the significant difference in APA score between gene clusters. APA gene cluster 1 had significantly higher APA score than APA gene cluster 2 (p < 2.2 × 10−16). The efficacy of the APA score in predicting ccRCC patients' prognosis was further evaluated. With the median value, patients were stratified into a high or low APA score group. Patients with a low APA score exhibited a markedly prominent advantage in OS (p = 8.326 × 10−4; Figure 6H), DFS (p = 7.847 × 10−4; Figure 6I), DSS (p = 3.061 × 10−4; Figure 6J), and PFI (p = 6.926 × 10−6; Figure 6K). To better illustrate the clinical characteristics of APA score, we examined the association between APA score and age, gender, grade, and stage of ccRCC patients (Figure 6L). No significant difference in APA score was found between age <65 and ≥65 patients. Male patients had a markedly higher APA score compared with female patients (p = 0.0096). With the increase of grade and stage, the APA score gradually increased. Furthermore, time-independent ROC curves confirmed the APA score showed a good permance in predicting ccRCC prognosis (Figure 6M). Multivariate Cox regression analysis confirmed APA score as an independent predictor for assessing patients' OS (hazard ratio [HR] 1.017 [1.003–1.032] and p = 1.78 × 10−2) and DFS (HR 1.018 [1.000–1.036] and p = 4.48 × 10−2), as shown in Figure 6N. To facilitate a personalized clinical application, we established a nomogram, including independent predictors (APA score, stage, and grade) for predicting 1-, 3-, and 5-year survival times (Figure S10A). The ROC, calibration, and decisive curves confirmed the excellent predictive efficacy of this nomogram in ccRCC patients' OS (Figures S10B–S10H).

Figure 6

Development of APA scoring system and characteristics of clinical traits in APA-related phenotypes

(A–D) Biological processes, cellular components, molecular functions, and KEGG pathways involved in the APA-related genes. (E) Survival analysis of patients in APA gene cluster 1 and 2 (log rank test). (F) The alluvial diagram for the attribute change of APA RNA processing pattern, APA gene cluster, APA score, and survival status. (G) The distribution of APA score in APA gene cluster 1 and 2. (H–K) OS, DFS, DSS, and PFI analysis for patients in high and low APA score groups. (L) The differences in APA score in different groups stratified by age, gender, grade, and stage. (M) Time-independent ROC curves of age, stage, gender, grade, and APA score. (N) Multivariate Cox regression analysis of age, stage, gender, grade, and APA score for OS and DFS.

Characteristics of APA score in TME immune regulation

High APA score had the characteristics of increased immune score (p = 0.0017; Figure 7A) and decreased stromal score (p = 0.00054; Figure 7B). As the APA score increased, the expression of HLA genes showed a decreasing trend (Figure 7C). By the ssGSEA method, we determined the immune cell types in ccRCC and compared the differences in immune cells between high and low APA scores. Our results showed that a high APA score was characterized by adaptive immune cell infiltration and immune activation, such as activated B cells, activated CD4 T cells, and activated CD8 T cells, whereas a low APA score was characterized by innate immune cell infiltration, such as eosinophils, immature dendritic cells, mast cells, natural killer cells, neutrophils, and plasmacytoid dendritic cells (Figure 7D). Furthermore, we found that the mRNAs of TNFSF9, PDCD1, LAG3, CD70, TMIGD2, and TNFRSF18 were up-regulated in the high APA score group whereas other immune checkpoints were up-regulated in the low APA score group (Figure 7E). As shown in the GSEA, oncogenic pathways like ERBB, JAK-STAT, MAPK, mTOR, pathway in cancer, TGF-β, and WNT pathways were distinctly activated in a low APA score (Figure 7F). Spearman correlation analysis also demonstrated the negative correlation between APA score and stroma activation-related pathways (WNT and EMT; Figure 7G). This indicated that a low APA score was distinctly correlated to stromal activation.

Figure 7

Characteristics of APA score in TME immune regulation

(A and B) Differences in immune score and stromal score between high and low APA score groups. (C–E) Differences in HLA mRNA expression, infiltration levels of 28 immune cells, and immune checkpoint mRNA expression between two groups. (F) GSEA for the activated KEGG pathways in low APA score group. (G) Correlation between APA score and the known signatures. Ns: not significant; ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001.

Characteristics of APA score in tumor somatic mutation

The distribution of somatic mutation in high and low APA score groups was visualized. VHL, PBTM1, and TTN were the most frequently mutated genes in the two groups (Figures 8A and 8B). A high APA score was significantly correlated to higher somatic copy-number alteration (SCNA) (p = 0.0227; Figure 8C). We also evaluated the MSI difference between groups. In Figure 8D, increased MSI was found in the high APA score group (p = 0.00055). The TMB quantification analysis showed that the high APA score group exhibited more extensive tumor somatic mutation compared with the low APA score group (p = 0.08; Figure 8E). Also, high somatic TMB indicated an unfavorable prognosis of ccRCC patients (p = 6.466 × 10−3; Figure 8F). With stratified patients by APA score and TMB, we found that patients with low TMB and low APA scores exhibited a distinct survival advantage (p = 4.504 × 10−4; Figure 8G).

Figure 8

Characteristics of APA score in tumor somatic mutation across ccRCC

(A and B) Landscape of somatic mutation in high and low APA score groups. (C–E) Differences in SCNA, MSI, and TMB between groups. (F) Survival differences between high and low TMB score groups (log rank test). (G) Survival differences among four groups stratified by APA score and TMB score (log rank test).

Estimation of candidate drugs for ccRCC patients based on APA score

The response to common targeted drugs sorafenib, sunitinib, pazopanib, and axitinib was estimated in the high and low APA score groups via the Genomics of Drug Sensitivity in Cancer (GDSC) database (Figures 9A–9D). We found that low APA score patients were more sensitive to sorafenib (p < 1 × 10−4), pazopanib (p = 0.00128), and axitinib (p = 0.000783). By CTRP and PRISM projects, candidate compounds were predicted for high APA score patients. This study identified two CTRP-derived compounds (including brivanib [r = −0.56] and ouabain [r = −0.76]); Figures 9E and 9F) and ten PRISM-derived compounds (including YM-201636 [r = −0.66], bosutinib [r = −0.62], GSK1904529A [r = −0.61], CNX-774 [r = −0.60], PD-0325901 [r = −0.45], mesna [r = −0.59], ciclesonide [r = −0.59], halobetasol-propionate [r = −0.57], Ro-4987655 [r = −0.49], and PD-0325901 [r = −0.45]; Figures 9G and 9H). These compounds had increased estimated AUC values in the high APA score group and were negatively correlated to APA score. The therapeutic potential of these compounds was verified through CMap analysis. PD-0325901 (CMap score <−65), in which gene expression increased in ccRCC tissues but reduced by treatment with certain compounds, possessed potential therapeutic utility in ccRCCs with a high APA score (Figure 9I).

Figure 9

Identification of candidate agents with high drug sensitivity in ccRCC patients based on APA score

(A–D) Drug sensitivity to sorafenib, sunitinib, pazopanib, and axitinib between high and low APA score patients. (E and F) Spearman's correlation analysis and differential drug response analysis of two CTRP-derived compounds in high APA score patients. Lower AUC value indicates higher drug sensitivity. (G and H) Spearman's correlation analysis and differential drug response analysis of ten PRISM-derived compounds in high APA score patients. (I) Validation of the therapeutic utility of the above compounds by CMap analysis. ∗∗∗p < 0.001.

Crosstalk between APA and AS modulates survival outcomes and tumor immune microenvironment

To explore the interaction between APA and AS, we calculated the correlation between APA score and AS score among ccRCC samples and found that there was a significantly positive correlation between APA and AS (R = 0.45 and p < 2.2 × 10−16; Figure 10A). Survival analysis showed that patients with high AS and high APA scores had the poorest survival outcomes whereas those with low AS and low APA scores had the most prominent survival advantage (p = 3.356 × 10−8; Figure 10B). An alluvial diagram depicted the changes of AS and APA scores in ccRCC samples (Figure 10C). Figure 10D visualized the landscape of AS and APA events across ccRCC patients. Antitumor lymphocyte cells exhibited the highest infiltration levels in patients with high AS and high APA scores and had the lowest infiltration levels in those with low AS and low APA scores (Figure 10E). Furthermore, most of HLA genes and immune checkpoints had the highest mRNA expression in the high AS and high APA score groups and had the lowest mRNA expression in the low AS and low APA score groups (Figure 10F). These indicated that the crosstalk between APA and AS may be important for modulating survival outcomes and TME immune regulation in ccRCC.

Figure 10

Crosstalk between APA and AS modulates survival outcomes and tumor immune microenvironment

(A) Correlation between AS score and APA score across ccRCC samples. (B) Survival analysis among four groups stratified by AS score and APA score (log rank test). (C) Alluvial diagram for the changes of AS score and APA score in ccRCC samples. (D) The landscape of AS and APA events across ccRCC patients. (E and F) Heatmap for the infiltration levels of tumor-infiltrating immune cells and the mRNA expression of HLA genes and immune checkpoints among four groups.

Discussion

Alternative RNA processing mechanisms, including AS and APA, are increasingly recognized as important regulators of proteome expansion and gene regulation., Growing evidence suggests that both AS and APA events play prominent roles in inflammatory response, innate immunity, and antitumor effects., For instance, immune-related genes with APA events in the TME may be predictive of risk stratification and clinical outcomes for patients with grade II/III gliomas. Because most studies center around individual tumor-infiltrating immune cells, the overall TME infiltration features induced by integrated implications of AS and APA events remain unrecognized. Nevertheless, the influence of AS and APA events in ccRCC has not been fully elucidated. Hence, identification of the roles of different AS or APA RNA processing patterns in the TME immune regulation may enhance the understanding of anticancer immune response, which may guide more effective immunotherapeutic strategies.

Pre-mRNA APA and AS play important roles during eukaryotic gene expression. Herein, on the basis of the top 1,000 AS or APA events with the most variance, we separately constructed two AS or APA RNA processing patterns. These patterns were characterized by distinct survival outcomes, immune response, immune cell infiltration, and activation of signaling pathways. This indicated that AS or APA events were of great significance in reshaping distinct TME landscape. Thus, comprehensive analysis of AS or APA patterns may strengthen the understanding of TME immune regulation. Furthermore, the mRNA transcriptome differences between distinct AS or APA RNA processing patterns were closely related to immune-related pathways. By the Boruta algorithm, we determined AS- or APA-related genes to develop the AS or APA scoring system. The AS or APA score was quantified to reflect the AS or APA patterns of an individual tumor. Both high AS and APA scores were indicative of undesirable survival outcomes. Our multivariate Cox regression analysis confirmed AS or APA score as robust and independent prognostic indicators of ccRCC patients. To facilitate clinical application, we developed an AS or APA score-based nomogram for predicting ccRCC prognosis. Moreover, a high AS or APA score was characterized by low mutation burden and high immune infiltration, indicative of higher sensitivity to immunotherapy. Hence, our AS and APA scores displayed a predictive advantage in precision immunotherapy in ccRCC. Our further analysis demonstrated that AS and APA scores could be also predictive of the efficacy of adjuvant chemotherapy. Due to poor prognosis for patients with a high AS or APA score, we screened several potential small molecular compounds against these patients. However, more experiments will be presented to verify the therapeutic effects of these compounds in ccRCC.

Although polyadenylation and splicing were initially considered different and independent events, it has recently been recognized that these mRNA processing events are all intertwined and interrelated. For instance, in some cases, APA is known to couple with AS to affect last intron removal. The remarkable correlation between AS and APA scores suggests that interactions between AS and APA may provide a shared activation mechanism for mRNA 3′ processing, splicing, and potentially other steps in RNA metabolism. Our study found that the interplay between AS and APA events played indispensable roles in prognosis and reshaping TME for ccRCC. These indicated that the crosstalk between APA and AS may be important for the generation of different RNA processing patterns between individual tumors.

Collectively, this study for the first time identified the roles of AS and APA events in ccRCC progression. Our findings uncovered the extensive regulatory mechanisms by which they affected the tumor immune microenvironment and their relationships with ccRCC survival outcomes. The AS and APA scores were separately constructed to quantify AS and APA RNA processing patterns in individual ccRCCs and identified their therapeutic utility in immunotherapy and targeted therapy. The interplay of the two posttranscriptional mechanisms displayed the prominent clinical implications of ccRCC. Our findings could help exploit personalized therapeutic strategies against ccRCC patients.

Materials and methods

Data collection and preprocessing

The workflow of this study was depicted in Figure 11. RNA sequencing (RNA-seq), somatic mutation, SCNA, and clinical information of ccRCC samples were downloaded from TCGA (https://portal.gdc.cancer.gov/repository). APA events in transcriptome data were assessed according to the PDUI value. The PDUI values of all genes in ccRCC samples were retrieved from The Cancer 3′ Untranslated Region (UTR)Atlas (TC3A; http://tc3a.org/). The PDUI value (range: 0–1) indicated the frequency of APA events. The PDUI value was proportional to the distal polyadenylation site of the transcript. AS events indicated as PSI in ccRCC specimens were downloaded from the SpliceSeq database (https://bioinformatics.mdanderson.org/TCGASpliceSeq/). The PSI value (range: 0–100%) represented an intuitive ratio to calculate the splicing efficiency of a gene sequence into a transcript isoform. For generating a reliable set of AS events, a percentage of samples with a PSI value ≥75 was set as the cutoff value. The missing value of PDUI or PSI value was filled using the K-nearest neighbors (KNN) algorithm. Finally, data containing 8503 APA events in 529 ccRCC samples as well as 20,802 AS events in 516 ccRCC samples were retrieved in our study.

Figure 11

Overview of study design

Consensus molecular clustering analysis

The PDUI and PSI values were separately curated and APA and AS RNA processing patterns were clustered via the NMF algorithm based on the top 1,000 AS or APA genes with the most variance. The optimal k value was determined based on cophenetic, dispersion, and silhouette coefficients. The t-SNE was applied for evaluating the accuracy of classification. The clustering patterns were validated based on the top 5,000 AS or APA genes with the most variance or MAD.

Gene set variation analysis (GSVA)

GSVA enrichment analysis, a non-parametric and unsupervised method, was employed to estimate variations in pathway activity across ccRCC samples. The gene set of c2.cp.kegg.v6.2.symbols was retrieved from the Molecular Signatures Database (MSigDB). The screening criteria of significant pathways were |log2fold-change|>0.2 and adjusted p < 0.05.

Quantification of immune response predictors

The mRNA expression of HLA family genes and immune checkpoints was quantified in each ccRCC sample. Stromal and immune scores were calculated to infer the fractions of infiltrating stromal cells and immune cells in ccRCC tissues utilizing the Estimation of STromal and Immune cells in MAlignant Tumors using Expression data (ESTIMATE) algorithm. TMB and MSI status were calculated to predict the response to immunotherapy.

Estimation of tumor-infiltrating immune cells

The ssGSEA algorithm was utilized for quantifying the relative abundance of 28 tumor-infiltrating immune cells in ccRCC samples. The gene set for marking each immune cell type was retrieved from the study of Charoentong et al., The enrichment scores were calculated to represent the relative abundance of tumor-infiltrating immune cells in ccRCC specimens.

Prediction of chemotherapeutic response

The response to targeted drugs was estimated through the GDSC (https://www.cancerrxgene.org/) database. Four commonly applied drugs against ccRCC were selected, including sorafenib, sunitinib, axitinib, and pazopanib. By employing the pRRophetic package, the IC50 was assessed with ridge regression analysis. The prediction accuracy was estimated with 10-fold cross-validation.

Generation of AS or APA score

By the limma package, differential AS and APA events were determined between AS or APA RNA processing patterns. The differential AS and APA events with adjusted p < 0.05 were extracted for feature selection by the Boruta package. Gene ontology (GO) and KEGG enrichment analysis of genes within final determined AS and APA events was performed with the clusterProfiler package, with the cutoff value of adjusted p < 0.05. The PDUI or PSI values of the selected AS or APA events were curated to perform principal-component analysis (PCA). Principal components (PCs) 1 and 2 were extracted and used as the signature score. As previously reported, the formula was employed to define the AS or APA score: AS or APA score = ∑(PC1i + PC2i), where i represented the PDUI or PSI value of the selected AS or APA event.,

Gene set enrichment analysis (GSEA)

GSEA was presented for elucidating the molecular mechanisms involved in the AS or APA score. GSEA was run through javaGSEA 3.0 on the basis of the MSigDB project. The c5.bp.v6.2.symbols.gm gene set was retrieved to identify enriched KEGG pathways. Pathways with false discovery rate <0.05 were considered significantly enriched.

Correlation between AS or APA score and several biological processes

A set of gene sets that stored genes correlated to several biological processes was retrieved according to previous studies, including CD8 T effector, DNA damage repair, pan-fibroblast TGF-β response signature (pan-F-TBRS), antigen processing machinery, immune checkpoint, EMT markers (EMT1, EMT2, and EMT3), FGFR3-related genes, KEGG discovered histones, angiogenesis, Fanconi anemia, cell cycle, DNA replication, nucleotide excision repair, homologous recombination, mismatch repair, WNT target, and cell-cycle regulators.45, 46, 47 Pearson correlation analysis between AS or APA score and the above biological pathways was then presented.

Prediction of drug response

Drug sensitivity data of human cancer cell lines (CCLs) were obtained from the CTRP (https://portals.broadinstitute.org/ctrp) and PRISM (https://depmap.org/portal/prism/) datasets, which may provide the AUC values to evaluate the response to different compounds. Through the KNN method, the missing AUC values were imputed. The lower the AUC value, the higher the sensitivity to a specific compound. Because the CCLs in both datasets were obtained from the Cancer Cell Line Encyclopedia (CCLE) database (https://portals.broadinstitute.org/ccle/), the expression profiles in CCLE database were adopted for CTRP and PRISM analysis.

Assessment of mRNA expression-based stemness index (mRNAsi)

Cancer stemness of ccRCC was quantified as described by Malta et al. The mRNAsi was determined through a one-class logistic regression machine learning algorithm, ranging from 0 (no gene expression) to 1 (complete gene expression).

Statistical analysis

All analysis was generated by R 3.6.2. Kaplan-Meier curves and log rank tests were utilized to assess the differences in OS, DFS, DSS, and PFI between groups via the survminer package. Comparison between two groups was analyzed by student's t test or Wilcoxon test, while one-way analysis of variance or Kruskal-Wallis test was employed for comparison among three groups or more. The ROC curves were conducted to evaluate the prognosis prediction efficacy of AS score, APA score, age, stage, gender, and grade. Independent prognostic factors for OS and DFS were screened with multivariate Cox regression analysis, which were included for a nomogram model that could predict the 1-, 3-, and 5-year survival probability via the rms package. The predictive accuracy of the nomogram was assessed through ROC, calibration, and decisive curves; p < 0.05 was considered statistically significant.