Abnormal Pseudouridine-related Pathways
In this study, we downloaded the set of genes associated with pseudouridine processes from MSigDB and quantified the degree of activation of these processes using the ssGSEA algorithm. Abnormalities in five pathways, rRNA pseudouridine synthesis, pseudouridine synthase activity, pseudouridine synthesis, mRNA pseudouridine synthesis, and tRNA pseudouridine synthesis, were observed in 30, 28, 28, 28, and 24 of 31 tumors (MSEO and UVM without control tissues), respectively, compared to nontumor tissues. The ssGSEA scores of these five pathways showed a significant increase in 20 tumors, including glioma, lung, and CRC, suggesting aberrant activation of these pathways. However, these processes behaved inversely in LAML and PCPG (Fig. 1A). Interestingly, a low activation state of all five pseudouridine-related processes was associated with an OS benefit in SKCM. KIRC patients with consistently reduced activity of pathways for pseudouridine synthesis, pseudouridine synthase activity, and tRNA pseudouridine synthesis had a better prognosis. LIHC patients with high levels of pseudouridine synthase activity, tRNA pseudouridine synthesis, rRNA pseudouridine synthesis, and pseudouridine synthesis had lower OS. In addition, these pathways were also found to be associated with OS in BRCA, LGG, LUAD, LUSC, PAAD, and SARC (Fig. 1B and S1A-G). In conclusion, pseudouridine-associated processes are aberrant in tumors and affect OS in tumor patients.
Altered Expression Patterns of PUSs
To further explore the possible reasons for the abnormal pseudouridine process, we investigated the key enzymes regulating pseudouridine. First, in the TCGA dataset alone, we observed more than 11 tumors with varying increases or decreases in the expression of 8 PUSs compared to peritumoral tissues. In particular, PUS7 was significantly abnormal in 21 tumors. It is noteworthy that all eight genes were significantly increased in BRCA, HNSC, KICH, LUSC, and STAD (Fig. 2A, B, and S 2A-F). To further elucidate the changes in the expression patterns of the regulators, we explored the differential mapping of these PUSs in 31 cancer types (MSEO and UVM without control tissue) in conjunction with the GTEx database. In the heatmap, we can see that at least two regulators are abnormal in each tumor. The most notable aberrant expressions are DKC1 and TRUB1, which are reflected in differences in 29 of the 31 cancer types. Focusing on this, all PUSs remained elevated in CHOL, COAD, LUSC, and READ tumor tissues. In contrast, in LAML, the expression of seven of the eight genes was significantly lower in the tumor samples, with only TRUB1 showing higher expression. In conclusion, synthetases involved in the pseudouridine process were aberrantly expressed in almost all cancer types, contributing to the dysregulation of this RNA modification process (Fig. 2C).
Genomic Aberrations of PUSs
Alterations in genetic information cause irreversible tumor development. Gene expression landscapes are often associated with copy number variation (Sun et al., 2018). We found that PUS1, PUS7, and PUS10 maintained consistent amplification in most tumor types. In contrast, TRUB1 and RPUSD4 maintained CNV deletions in most types (Fig. S3A). Correlation analysis showed that CNV amplification was associated with high expression of most PUSs. In particular, the altered expression of PUS7 was attributed to CNV amplification, but the expression of DKC1 appeared to be less affected by CNV variation (Fig. 2D).
In addition, we explored mutations in pseudouridine-regulated proteases. Mutations were present in 394 of 10,540 samples from 33 tumors, with the highest number of patients having mutations in the PUS7 gene at 107 (Fig. 2E). In UCEC, PUS7 reached 5%, and TRUB2 reached a minimum of 1% (Fig. S3B). All genes were dominated by missense mutations.
During tumorigenesis, key oncogenes are regulated by DNA promoter methylation, and aberrant methylation leads to silencing of gene expression (Das and Singal, 2004). By comparing cancer tissues with normal tissues (tumor types with greater than 10 normal tissues in the TCGA database), we found that the majority of PUSs showed dysfunctional methylation levels in tumors. Among them, DKC1 showed a significant hypomethylation status in several cancer tissues, including LIHC, HNSC, UCEC, PRAD, LUSC, THCA, BLCA, BRCA, and COAD (Fig. S3C). There is a clear and close link between aberrant methylation and the expression of pseudouridine regulatory factors, whereby the majority of PUS expression is negatively regulated by the degree of methylation. Notably, DKC1 promoter methylation showed a stable association with DKC1 mRNA expression in 24 tumor types (Fig. 2F). In contrast to CNV variants, DKC1 expression was more inclined to be affected by aberrant methylation.
Correlation between PUSs Expression and Tumor Prognosis
To elucidate the impact of aberrantly expressed PUSs on tumor prognosis, we performed a K‒M survival analysis. The expression of PUSs in 21 out of 33 tumors correlated with OS, of which reduced expression of PUSs in ACC, CESC, ESCA, HNSC, KICH, KIRP, LIHC, MESO, PAAD, PRAD, SARC, UCEC, and UVM tumors was favorable for OS (Fig. 3A). The disease-specific survival (DSS) of an additional 20 tumors was affected by the expression of PUSs. The most affected tumor was PRAD, where five synthases, DKC1, PUS3, RPUSD4, PUS7, and TRUB2, showed mostly negative effects on DSS. In addition, DKC1 and PUS7 were negatively correlated with DSS time in nine of the tumors (Fig. 3B). In contrast, the disease-free interval (DFI) of the tumors appeared to be the least affected. However, PUSs also maintained significant values in 14 tumors. Among them, increased expression of more than two PUSs negatively affected the DFI in ACC, KIRP, LIHC, and PAAD but positively affected CHOL, KIRC, LUSC, and THCA (Fig. 3C). In ACC, KIRP, and PRAD, the expression of more than 4 genes negatively regulated the progression-free interval (PFI) in tumor patients. PUS7 affected PFI in more than 10 tumors (Fig. 3D). Notably, OS, DSS, and PFI were all susceptible to being modulated by PUSs in LGG but not DFI. Furthermore, high TRUB1 expression tended to result in better OS, DSS, and PFI in both types of gliomas, LGG and GBM. Interestingly, OS, DSS, DFI, and PFI of ACC, KIRP, and LIHC were all susceptible to being affected by PUSs-regulated expression.
Subsequently, univariate Cox regression analyses showed that DKC1 could serve as an independent risk element for PRAD, KICH, LGG, KIRP, ACC, MESO, LIHC, SARC, UCEC, LUAD, and HNSC (Fig. 3E). PUS7 acted as an independent rumor factor for ACC, CESC, KICH, THCA, MESO, KIRP, LGG, PAAD, SARC, and LIHC (Fig. 3F). PRUSD4 plays a risk predictor in ACC, BLCA, LIHC, SARC, and PAAD tumors. On the other hand, PRUSD4, DKC1, and PUS7 were also prognostic protective agents for tumor READ. PRUSD4 was a protective factor for KIRC along with TRUB1, TRUB2, and PUS3 (Fig. 3G). In contrast, PUS1, when used as an independent risk factor, only affected the development of ACC, KICH, PRAD, KIRC, MESO, LGG, LIHC, SARC, and UVM tumors (Fig. 3H). In contrast, TRUB1 was only protective in tumors (Fig. 3I). PUS10 could act as a risk factor for LGG and UCEC and a protectant for SKCM (Fig. 3J). Increased PUS3 expression suggested a worse prognosis for ACC and HNSC (Fig. 3K). The dual predictive ability of TRUB2 was demonstrated by acting as a positive factor for DLBC and KIRP, and as a risk factor affecting the prognosis of ACC, PRAD, and UVM (Fig. 3L). These results suggest that these eight proteases regulating pseudouridine are closely related to the prognosis of tumor patients.
PUSs and Signaling Pathways
Abnormally expressed genes often originate from or result in oncogenic events. To explore the effects of PUSs on tumor signaling pathways, we calculated 50 characteristic pathway enrichment scores for 33 tumors by the ssGSEA algorithm. Spearman correlation analysis showed that the expression of all seven PUSs in more than 20 tumors was robustly and positively correlated with the ssGSEA enrichment scores of MYC-TARGETS_V1, MYC_TARGETS_V2, E2F_TARGETS, and MTORC1. In all tumor types, the expression of DKC1 was strongly correlated with the activation of MYC-TARGETS_V1. The expression of PUS1 also maintained a stable positive correlation with MYC_TARGETS_V2. However, the expression pattern of PUS10 differed from that of other regulators in that it showed a stable negative correlation with MYC_TARGETS_V2, which was validated in 20 tumors (Fig. 4A-H). The enriched fractions of MYC_TARGETS_V1, MYC_TARGETS_V2, E2F_TARGETS, and MTORC1 were significantly increased in tumors compared to normal tissues (Fig. S4). We speculate that the aberrant activation of these pathways may be related to the aberrant regulation of PUSs.
PUSs modulate the TME
The potential link between RNA modification and immune infiltration has gradually attracted attention (Yang et al., 2023). The TIMER algorithm estimated the infiltration of six immune cells in tumor tissues using RNA-Seq expression profiles of TCGA. We calculated the correlation between PUSs and the infiltration of these immune cells. Overall, we observed that PUSs were mainly involved in immune infiltration in solid tumors. For example, in kidney cancer, most of the PUSs, including KICH, KIRC, and KIRP, were significantly positively correlated with six immune cells. In lung cancer (LUAD, LUSC), DKC1, PUS1, PUS3, PUS7, TRUB1, and TRUB2 were mostly negatively correlated with immune cell infiltration. Whereas PUS3, PUS10, TRUB1, and TRUB2 maintained a similar correlation with GBM and LGG in terms of immune cell infiltration. In addition, high expression of PUS7 and PUS10 in PAAD, PCPG, PRAD, READ, TGCT, and THCA was associated with increased B-cell and macrophage infiltration. Interestingly, in UVM, increased expression of DKC1, PUS1, PUS3, PUS7, PUS10, and TRUB2 showed a significant correlation with increased CD8 + cell infiltration and decreased MDC infiltration only (Fig. 5A).
Given the significant impact of PUSs on the TME in lung cancer, we focused on the immune landscapes of the two lung cancers. An algorithm based on ssGSEA, which has a more detailed classification of immune cell infiltration, computed results that were similar to the TIMER algorithm. Overall, the expression of PUSs showed a negative correlation with the total infiltration score in LUAD and LUSC. It was positively correlated with the immunosuppressive cell nTreg and negatively correlated with the infiltration score of NK and CD4 + T cells with tumor-killing effects. Increased expression of pseudouridine modulates an immunosuppressive TME (Fig. 5B and C).
With this in mind, we explored the effect of PUSs on the expression of ICPs. We found that in renal cancer, PUSs significantly affect CD8 + and CD4 + T-cell infiltration. However, this was accompanied by an increase in protease expression, with a steady rise in the expression of PD-1, PD-1, and CTLA4. This may limit the killing effect of these immune cells. This phenomenon was also observed in tumor UVM, PAAD, and PCPG (Fig. 5D).
Construction of the pseudouridine signature score
To better describe the pseudouridine modification within the tumor, we constructed the pseudouridine feature score. Facing the statistical methodological difficulties in modeling high-dimensional data, we applied the integration of multiple machine-learning modeling strategies to improve the accuracy and robustness of the feature model score. First, we performed Pearson correlation analysis of all genes in 33 tumors with eight PUSs to screen genes with potential interactions. Upset plots showed that the expression of 168 genes was strongly associated with PUSs in all types of tumors (Fig. 6A). Meanwhile, 21 genes maintained the same correlation with the corresponding PUSs in all tumors (Fig. 6B). Based on the expression of these 21 genes and 8 PUSs, we trained 112 prognostic models for each tumor using an aggregated machine-learning approach and applied them to other tumors. The prognostic power was assessed based on the average C-index of each model in the tumor. Ultimately, we found that the ACC-based genes screened with gloost + nets22 had the highest average C index (0.590). We used this as a template to construct a prognostic model and applied it to ACC, DLBC, GBM, KICH, MESO, THYM, TGCT, and PRAD with good prognostic ability (C-index > 0.6) (Fig. 6C and S5A).
The pseudouridine signature scores we constructed correlated with the mRNA pseudouridine synthesis pathway in all tumors (p < 0.05). The tumor with the strongest degree of correlation was LAML, with a correlation coefficient of 0.87. In addition, the pseudouridine signature score correlated with pseudouridine synthase activity and pseudouridine synthesis ssGSEA scores in 30 tumors, i.e., the signature model also better describes these two functional processes in tumors (Fig. 6D).
Predictive Chemotherapy Drugs
We divided the different tumor samples into high and low groups based on the pseudouridine signature score. We found that the signature score responded significantly to the therapeutic effect of etoposide. In 29 tumors, the IC50 of etoposide in the high group was much lower than that in the low group, implying that the high group was more dependent on etoposide treatment (Fig. 7A). In addition, two drugs, camptothecin and methotrexate, were more suitable for patients with BLCA, BRCA, LAML, KIRC, HNSC, GBM, LGG, LUAD, SKCM, SARC, THYM, UCSC, and UVM tumors with high pseudouridine status (Fig. S6 and S7). Low-scoring BLCA, BRCA, HNSC, KIRP, LGG, PAAD, SARC, SKCM, TGCT, and UCSC were more tolerant to cisplatin (Fig. S8). The high subgroups of CESC, KIRC, MESO, PRAD, SARC, and THCA had lower IC50 values for the chemotherapeutic agent bexarotene and were more amenable to bexarotene treatment. The low subgroups LAML, LUAD, STAD, UVM, and THYM were more suitable for bexarotene treatment (Fig. S9). The above results suggest that the pseudouridine signature score has the potential to guide the selection of clinical chemotherapeutic agents for different tumors.
Stratification analysis and mutation status
Analysis of tumor mutation information allows us to better understand patients with different characteristics. In addition, somatic mutation rates are often associated with tumor treatment. The mutation rate was 69.05% in the low pseudouridine-scoring group and 72.3% in the high-scoring group. This implies that the higher-scoring group with a higher somatic mutation rate possesses a high release of tumor neoantigens and more potent anticancer immune activity (Bekri et al., 2022; Rooney et al., 2015). In addition, the gene with the highest mutation rate in both sample groups was TP53 (33% of the low-scoring group and 44% of the high-scoring group) (Fig. 7B, C). In contrast, the signature scores were significantly higher in TP53-mutated ACC, BLCA, BRCA, COAD, HNSC, LIHC, LUAD, LUSC, PRAD, READ, SKCM, STAD, UCEC, and UCS samples (Fig. S5B).
TMB has emerged as a quantitative genomic biomarker of response to ICB (Goodman et al., 2017). We found that our constructed pseudouridine signature score was stably and positively correlated with TMB in ACC, BLCA, HNSC, LAML, LGG, LUAD, LUSC, OV, PRAD, READ, SKCM, STAD, THCA, and THYM and negatively correlated with KIRC and KIRP (Fig. 7D). This suggests a potential association between the characterization model and ICB treatment response. This conjecture has been confirmed in the exploration of MSI, which has emerged as a "predictive" biomarker of immunotherapy benefit (Schrock et al., 2019). We observed a stable positive correlation between pseudouridine scores and MSI in BRCA, BLCA, HNSC, LUAD, LUSC, MESO, READ, SARC, and THCA (Fig. 7E).
Predictive ICB Responsiveness
ICB has yielded surprising results in the treatment of a wide range of tumors, including melanoma and lung cancer (Sun et al., 2023). Considering that there may be some connection between the signature model and immunotherapy, we performed further exploration with the help of the TIDE algorithm. TIDE provides a new computational framework to quantify tumor immune dysfunction and rejection based on gene expression profiles of tumor samples and to further assess the likelihood of tumor immune escape and response to ICB. With this in mind, we performed a stratified predictive analysis of the pseudouridine signature model for ICB therapy. The results suggest that the pseudouridine signature model is a good predictor of the effect of immunotherapy on tumors. ACC, ESCA, GBM, KICH, KIRC, KIRP, LAML, LGG, LIHC, LUAD, MESO, PAAD, PCPG, PRAD, SARC, SKCM, STAD, THCA, and UCEC patients with lower pseudouridine scores were more sensitive to ICB therapy (Fig. 7F).