Dysregulated genes were identified in patients with distinct TP53 mutational statuses
As mentioned in Materials and methods, after rounds of filtrations, a total of 365 tumoral and 50 normal HCC samples were involved and the tumor samples were further divided into TP53 MT and TP53 WT group. Notably, patients with distinct TP53 mutational status exhibited significantly different OS and disease-free survival (DFS) (Suppl. Fig. 1A-B), connoting the detrimental effects of TP53 mutations in HCC. By later DEG identification procedures, a number of DEM, DEL and DEmi were identified by DESeq2, limma and edgeR tools in samples of two mutational statuses (Suppl. Fig. 2A-F, Suppl. Fig. 3A-F). When further inspecting the tool-wise gene intersection, 666 DEMs, 440 DELs and 34 DEmis were recurrently found up-regulated in TP53 MT group than normal control (Suppl. Fig. 2A, C, E) while 419 DEMs, 335 DELs and 10 DEmis showed higher expression in TP53 WT samples that the control set (Suppl. Fig. 2B, D, F). Besides, 276 DEMs, 49 DELs and 15 DEmis were down-regulated in TP53 MT samples (Suppl. Fig. 3A, C, E) while 167 DEMs, 34 DELs and 17 DEmis presented an expressional demotion in TP53 WT patients than the normal samples (Suppl. Fig. 3B, D, F).
With the aim of gaining insights into the accompanying changes granted by TP53 abnormality, DEMs obtained from the two groups were subjected to functional annotations by GO, KEGG and Reactome databases. For the TP53 MT subgroup, most of the up-regulated genes were enriched in biological processes or pathways related to cell migration, angiogenesis and apoptosis (Fig. 2A-C) while activities of metal iron response and transport-related pathways were suppressed (Fig. 2A-C). Similar pathway enrichments for the dysregulated genes were observed in TP53 WT group (Suppl. Fig. 4A-C). Interestingly, immune-related terms were found up-regulated in two subgroups and PD-1 checkpoint-related pathway was concurrently impaired. Concerning the impotence of anti-PD-1/PD-L1 blockers in major HCC patients43, investigation on the TME could possibly provide additional therapeutic benefits for HCC.
Construction of TP53 alternation-specific ceRNA network aided the elicit of novel prognostic markers in HCC
By combining the expressional correlation and miRNA targeting information, TP53 mutational status-specific ceRNA networks were constructed based on the ceRNA hypothesis. As shown in Figure 3A, multiple genes were found involved in the regulations provoked by TP53 mutation. Interestingly, no ceRNA network was identified using the same construction strategy and DEGs in TP53 WT group, which affirmed the causal effect of TP53 abnormality on tumorigenesis or tumor progression. More explicitly, for the TP53 MT ceRNA network, lncRNAs (AC017104.6, RP5-1092A11.5 and RP11-365O16.6) and mRNA (PTBP1) both exhibited expressional elevation in tumor samples while the only miRNA (hsa-miR-133b) showed opposing trend (Fig. 3B). When further examining tumor samples stratified by TP53 status, significant expressional increment of the aforementioned lncRNA and mRNA was observed in TP53 MT subgroup (Fig. 3C-F), again emphasized the importance of our identified regulatory network in biological process decoding.
To better understand the functional importance of the ceRNA network arose from TP53 abnormality, we performed cox regression analysis for the screening of prognosis-related biomarkers. As illustrated by the forest plot in Figure 4A, the univariate cox regression analysis confirmed the independent prognostic relevance of PTBP1, AC017104.6, RP5-1092A11.5 and RP11-365O16.6 with OS, all of which reached Hazard Ratio (HR)> 1 and p-value < 0.05. However, when inspecting the relationship between nodes’ expression and DSS, only PTBP1 showed significant relevance (Fig. 4B). Next, we conducted survival analysis on these independent prognostic factors to unveil their potential predictive power on patient survival and confirm their regulatory significance in HCC. As expected, the survival of TP53 MT patients presented a negative trend with the expression of PTBP1, both for OS (Fig. 4C) and DSS (Fig. 4D). In addition, reverse associations between expression and patient OS were analogously observed on lncRNAs including AC017104.6 (Suppl. Fig. 5A), RP5-1092A11.5 (Suppl. Fig. 5B) and RP11-365O16.6 (Suppl. Fig. 5C). Since the activation the above four ceRNA network nodes led to poor prognosis, we hypothesized their functional collaboration in the context of TP53 mutations. Indeed, patients in TP53 MT group could be separated by expressional patterns (Suppl. Fig. 5D) and the patients with higher expressions exhibited worse OS (Suppl. Fig. 5E) and DSS (Suppl. Fig. 5F). Through the ceRNA network construction in TP53 MT group, novel HCC prognosis-related genes were identified, which could deepen our understanding of TP53 gene’s mechanical role in cancer.
The ceRNA network established by TP53 mutations was associated with tumor microenvironment alternations in HCC
In the previous functional annotations on TP53-related DEGs, immune-related terms were uncovered. We next scrutinized the relationship between TP53-provoked regulatory alternations and TME. To be more precise, the prediction as well as quantification of the immune cell types were initially performed by the xCell method on TP53 MT samples and correlations between ceRNA network nodes and immune infiltration levels were further analyzed. As shown in Figure 5A, infiltration level of common lymphoid progenitor significantly correlated with the expression of PTBP1, AC017104.6, RP5-1092A11.5 and RP11-365O16.6. Likewise, infiltrations of T cell CD4+ T helper (Th) 2, T cell CD4+ memory and mast cell also formed a correlation module with ceRNA node genes (Fig. 5A). To be more precise, among the four cell types that reached statistical significance, PTBP1 (Fig. 5B-E) and RP5-1092A11.5 (Suppl. Fig. 6) showed the highest SCC with T cell CD4+ Th2 infiltrations while the most relevant immune cell types for AC017104.6 and RP11-365O16.6 were common lymphoid progenitor (Suppl. Fig. 7) and mast cell (Suppl. Fig. 8). Such associations deepened the understanding of the crosstalk between TME and TP53 alternation-provoked ceRNA network and possibly could benefit the immunotherapy development in HCC.
We next investigated whether the association between HCC TME perturbation and ceRNA network was TP53 mutant-specific. By comparing the predicted levels of immune cell infiltration between TP53 MT and WT group, 15 immune cell types reached statistical significance (adjusted p-value<0.05, Wilcoxon rank sum test adjusted by Benjamini & Hochberg method) (Suppl. Fig. 9A). Interestingly, infiltration levels of the immune cell types associated with ceRNA network nodes including T cell CD4+ Th2, mast cell, and T cell CD4+ memory were remarkably different between the groups stratified by TP53 mutational status (Suppl. Fig. 9A) and the elevated infiltration of these cell types conferred pernicious effect on TP53 MT patient survival, both for OS (Suppl. Fig. 9B, D, F) and DSS (Suppl. Fig. 9C, E, G). To conclude, our analyses confirmed the TP53 mutant-specific association between ceRNA network and TME variations, which could pose a pernicious effect on patient survival.
Additionally, the association between ceRNA network nodes and cancer-immunity cycle activity was examined by ssGSEA. The co-expressed network nodes exhibited a high correlation with the activity of “infiltration of immune cells (T cells) into tumors” (Suppl. Fig. 10A). The subsequent evaluation confirmed the statistical significance of the relationships between the enriched step and network nodes’ expression levels (Suppl. Fig. 10B-E). However, when further comparing the activity of “infiltration of immune cells (T cells) into tumors” between TP53 MT and WT group, no obvious difference was observed. Combining with the discoveries above, TP53 mutant-specific network could possibly trigger the infiltration of specific T cell subtypes, which offered a boost to the development of precision immuno-medicine.
TP53 mutant-specific ceRNA network genes were connected with chemotherapy resistance
We next scrutinized the relationship between ceRNA network nodes and possible chemotherapy resistance. Through the correlation analysis between chemotherapy resistance-related genes and ceRNA network nodes, we found that PTBP1, RP5-1092A11.5 and AC017104.6 exhibited high correlation with respective sets of chemo-related genes (Suppl. Fig. 11A) while RP11-365O16.6 exhibited a weak connection with chemotherapy resistance. More detailedly, PTBP1 co-expressed with multi-facet chemotherapy resistance-related genes (Fig. 6A), covering all four categories we defined. As for RP5-1092A11.5, it was closely related to chemotherapy resistance caused by transmembrane transporters, critical (specific) pathways, lipid metabolism and HCC target therapy (Fig. 6B) while AC017104.6 demonstrated associations with chemotherapy resistance factors related to TME, cytokine, cancer stem cell, heat shock proteins, specific pathways and epigenetics (Fig. 6C). Unsurprisingly, among other participants in ceRNA crosstalks, mRNA PTBP1 presented the strongest correlation with chemotherapy resistance factors. Interestingly, the expression of these resistant genes gradually elevated in normal, TP53 WT and TP53 MT groups (Suppl. Fig. 11B), again reflected the biological complexity prompted by the TP53 mutations and emphasized the difficulty as well as cruciality of treatment development in HCC.
The DNA methylation level of PTBP1 and DNA methylation regulator expressions changed in HCC patient groups
Transcription of genes could be affected in multifarious ways. Intending to uncover the cause of PTBP1 expressional elevation in TP53 MT patients, we collected the DNA methylation levels of PTBP1-related probes and checked their dynamics in different patient groups. Among all PTBP1-related probes, 7 of them demonstrated an observable group-wise difference (Suppl. Fig. 12A) while only probe cg02086742 possessed methylation level that anti-correlated with the expressional trend of PTBP1 (Suppl. Fig. 12A). As expected, methylation levels of cg02086742 were significantly lower in TP53 MT group (Suppl. Fig. 12B), implying the potential of applying it as a target of epigenetic therapeutics in HCC.
Moreover, regarding their pivotal regulatory effects on DNA methylation, the expression of DNA methylation writer (DNMT1, DNMT3A and DNMT3B), eraser (TET1, TET2 and TET3) and reader (MeCP2, MBD1 and MBD2) proteins were compared between groups. We only focused on proteins with significant group-wise difference as well as alternation trend conforming to their biological characteristics. As a result, two eraser proteins TET1 (Suppl. Fig. 12C) and TET3 (Suppl. Fig. 12D) met the selection criteria and were strongly correlated with the methylation of probe cg02086742 and expression of PTBP1 in three patient groups. These observations provided a possible novel strategy for targeted therapy in HCC.
Patients stratified by TP53 mutational status exhibited differences in drug sensitivity
With the aim of facilitating systematic therapy development, analyses on immunotherapy, chemotherapy and targeted therapy were conducted above. Next, we investigated the drug response diversity in HCC TP53 MT and WT groups. By utilizing the drug sensitivity data in GDSC, the half-maximal inhibitory concentration (IC50) value of patients in the two groups were predicted. Among all drugs tested, 7 drugs demonstrated a group-wise statistical difference (p-value<0.05) (Figure 7A). More specifically, drugs AZD8055, OSI.027 and RO.3306 showed higher effectiveness in TP53 MT patients (Fig. 7B-D), among which AZD8055 got the lowest IC50 values (around 1). Besides, the group-wise ID50 value difference exacerbated for the OSI.027 compound (Fig. 7C). As for drugs JQ1, NU7441, SB216763 and ZM447439, their IC50 values were significantly lower in TP53 WT group (Suppl. Fig. 13A-D). Compound SB216763 presented relatively higher IC50 values, while JQ1 got the lowest IC50 values in TP53 WT HCC, which warranted its therapeutic strength.
Multi-omics validation of PTBP1 in public databases
To verify the universality of our discoveries and eliminate possible biological biases, expressional data from ICGC and CTRP databases were collected. For ICGC database, the RNA expression level of PTBP1 in TP53 MT samples was significantly higher than the normal control group (Fig. 8A). When examining the protein expressional levels in CTRP database, though not reaching statistical significance, a certain reduction in TP53 WT samples was observed (Fig. 8B). These verifications shared high consistency with our discoveries in ceRNA network analyses and underpinned the connection between TP53 mutations and PTBP1 alternations.
In addition, the immunohistochemistry of PTBP1 in liver cancer was retrieved from the HPA database to provide an experimental point of view on protein expression. As illustrated in Figure 8C, the prevalence of PTBP1 expression in HCC was confirmed, again indicated its potential in HCC patient stratification and treatment.