Comprehensive analysis of a ceRNA network associated with TP53 mutations in liver hepatocellular carcinoma

doi:10.21203/rs.3.rs-1461119/v1

PDF

Article

Comprehensive analysis of a ceRNA network associated with TP53 mutations in liver hepatocellular carcinoma

https://doi.org/10.21203/rs.3.rs-1461119/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted 17 Mar, 2022

You are reading this older preprint version

Read the latest preprint version →

Liver cancer is one of the most common causes of cancer-related death and is quantitatively dominated by the liver hepatocellular carcinoma (HCC) subtype. Although advances have been made in recent years, further exploration on new HCC carcinogenic mechanisms as well as diagnostic and therapeutic methods in the context of detrimental TP53 abnormalities is indispensable on the way to conquer HCC. In this study, competing endogenous RNA (ceRNA) network was constructed on tumor samples stratified by the TP53 mutational status. With demonstrated TP53 mutant-specificity, the network nodes including AC017104.6, RP5-1092A11.5, RP11-365O16.6 and PTBP1 were all prognosis-related. Further therapeutic correlation analysis unveiled the connection between their expressions and tumor microenvironment fluctuation as well as possible chemotherapy resistance. Additionally, we explored the drug sensitivity difference between TP53 wild type and mutant samples, complementing the therapeutic arsenal of HCC. The universality of the key node PTBP1 was verified by multi-source transcriptional and protein level expressions. Taken together, this is the first study focusing on the ceRNA network mediated by TP53 mutations and its provoked therapeutic associations in HCC. Our results deepened the understanding of the concomitant alternations brought by TP53 mutations and could benefit future therapeutic development in HCC.

Being the second most archetypal deadly cancer worldwide¹, liver cancer prevailed in underdeveloped areas and emerged as the fifth and ninth most common cancer in males and females respectively. When narrowing to Asian countries, the severity of the disease exacerbated due to the frequent chronic viral hepatitis infection or exposure to toxins. For example, the 5-year survival rate is as low as 12% in China². Hepatocellular carcinoma (HCC), one cancer type originated from liver cells, accounted for more than 80% of global primary liver cancer amount³ and standards of its diagnosis and treatment have been increasingly improved. For early-stage HCC, liver resection and transplantation were proven to be the surpassing treatment option. Disastrously, HCC early-stage surveillance and diagnosis still remains a challenge, resulting in the majority of HCC patients at initial diagnosis was at the intermediate stage. With the presence of vascular invasion and extrahepatic spread, curative treatment became impotent in late stage HCC⁴ and systemic therapy was proven an indispensable surrogate treatment option. Comparing to treatments like radiotherapy or chemotherapy, targeted therapy is now a more promising option for HCC⁵ regarding its improved efficacy. Multiple molecular therapies including tyrosine kinase inhibitors (TKIs), monoclonal antibodies and antisense oligonucleotides have been developed and researchers are paying indomitable efforts in continuous targeted drug development for live cancer.

Encoded by the human TP53 gene, P53 protein has long been known as the guardian of the genome, regarding its role in DNA repair and recombination⁶. Unfortunately, somatic mutations in TP53 are one of the most common alternations in human cancer, which could demolish the genome stability, enhance tumor growth, affect the metabolic function of mitochondria⁷ and cause drug resistance⁸. Unlike other tumor suppressor genes (TSGs) that were often debilitated by truncated mutations or deletions, missense mutations dominated the TP53 alternation landscape, producing full-length protein with a single amino acid substitution. In HCC, multiple oncogenes or TSGs including TERT, TP53, CTNNB1, etc⁹ were reported to harbor tumorgenetic mutations, among which the association between TP53 mutated phenotype and shorter overall survival (OS) and recurrence-free survival (RFS) was previously revealed¹⁰. Considering the ubiquitousness and severeness of TP53 mutations, researchers have been dedicated to developing drugs targeting TP53 aberrance. However, cancer therapies targeting mutant TP53 are mostly being tested in early phase clinical trails¹¹, hampering their extensive application in large and diverse populations. Owing to the inherent functional diversity of TP53 and its multi-faceted impact on tumor progression, the deriving of TP53 mutation concomitant biomarkers could undeniably facilitate the development of combination cancer therapy and benefit patients.

Recently, we have witnessed a surge of investigation into novel biomarker prediction for precision medicine, especially through bioinformatics approaches. Among multifarious biomarker identification methods based on mathematics, network topology and machine learning theories, the construction of competing endogenous RNA (ceRNA) networks has received particular attention due to its high mechanismic consistency with biological hypothesis and the affluence in knowledge discovery. Through the establishment of ceRNA network, the crosstalk between microRNA (miRNA), long non-coding RNA (lncRNA) and messenger RNA (mRNA) could be discovered, providing an unprecedented layer of miRNA regulatory network. Increasing evidence have shown that the ceRNA crosstalk participates in indispensable regulatory processes and its perturbation will scramble the balance of cell regulatory system, resulting in disease initiation and tumor occurrence. In HCC, ceRNA networks aided the prioritization of RNA biomarkers including miRNA¹² and lncRNA¹³ by comparing normal controls with tumor patients. As previously mentioned, the regulatory peculiarity of TP53 granted its capability in orchestrating complicated biological processes. By incorporating TP53 mutational information in ceRNA network construction, the consequent genetic fluctuations caused by TP53 could be uncovered and compared, which will provide unseen therapeutic targets and pave the way to systematic therapy advancement.

In this study, HCC expressional data was collected from The Cancer Genome Atlas (TCGA) database and stratified by the TP53 gene mutational status. By using these profiling data, a ceRNA network specifically existed in TP53 mutant samples was primarily established. After confirming the prognostic relevance of the network nodes, analyses on immune cell infiltration, chemotherapy resistance and drug response were further conducted. Generality of the key mRNA node was guaranteed by multi-source database validations. Clinically, these gene nodes could possibly benefit the development of systemic therapy as well as aid the treatment outcome prediction in HCC and improve the unsatisfactory prognosis of HCC patients.

TCGA data preparation

Firstly, the liver hepatocellular carcinoma (HCC) HTSeq-Count file from RNA-seq and isoform expression quantifications from miRNA-seq were retrieved from TCGA GDC Data Portal (https://portal.gdc.cancer.gov/) while the clinical information was downloaded from cBioPortal (https://www.cbioportal.org/) database. Patients met the following criteria were further discarded: i) harboring TP53 synonymous mutations; ii) lacking RNA and miRNA sequencing data. Additionally, we used miRBase¹⁴ MIMAT ID for miRNA quantifications and only mature miRNAs were kept and the summed read counts were calculated for further use.

Differentially-expressed gene identification procedure

The derivation of differentially-expressed mRNAs (DEMs), lncRNAs (DELs) and miRNAs (DEmis) were separately conducted between 365 primary solid tumor patients with two TP53 mutational status (111 TP53 mutant and 254 wild type) and 50 normal samples. More specifically, GENCODE v22 annotation files were primarily downloaded from GENCODE database (https://www.gencodegenes.org/) for gene type, length and gene symbol information retrieval. After eliminating lncRNAs, miRNAs and mRNAs with sample-wise average read count<1, a total of 8782 mRNAs, 3609 lncRNAs and 2166 miRNAs were retained for further analyses. For RNA-seq data, read counts of lncRNAs and mRNAs were normalized by TCGAbiolinks R package¹⁵ and DESeq2, limma and edgeR tools^16–18 were utilized for differentially-expressed gene (DEG) identification between TP53 mutant type (MT) and normal as well as TP53 wild type (WT) and normal samples. We kept the genes with p-value<0.05 and absolute value of log2(Fold Change)>1 as candidate DEGs. Regarding the intrinsic statistical differences between the algorithms, the DEG candidates from tools were intersected to get convergent lists for two TP53 mutational statuses. For miRNA-seq data, quantile normalization using TCGAanalyze_Filtering function was initially conducted and DESeq2, limma and edgeR tools similarly facilitated the DEmi derivation. Identical thresholds on p-value and log2(Fold Change) were applied on miRNAs and gene lists were intersected for further use. Conversions between MIMAT IDs and miRBase identifiers were completed by miRBaseVersions.db package (https://github.com/StefanHaunsberger/mirbaseversions.db) and statistics on the gene intersection between the three tools were performed by VennDiagram R package¹⁹.

Functional annotation on gene sets

Pathway enrichment analyses on the convergent DEM lists of two mutational statuses were conducted using Enrichr R package²⁰, which integrated comprehensive gene-set annotations. Resources including KEGG pathway²¹, Gene Ontology (GO) Biological Processes²² and Reactome database²³ were adopted in our enrichment analyses. Enriched GO terms with p-value<0.01 and enriched KEGG and Reactome pathways with p-value<0.05 were reported.

TP53 mutation-specific ceRNA network construction

In accordance with the definition of ceRNAs which regulates the expression of mRNAs by competing for shared miRNAs, the ceRNA network was constructed using following procedures: i) DEMs, DELs and DEmis from the two TP53 mutational statuses were collected and co-expression networks were firstly constructed based on the expressions of DEM/DEL, DEL/DEmi and DEM/DEmi pairs. Only pairs with Pearson Correlation Coefficient (PCC)>0.6 in DEM/DEL network and negatively-correlated pairs in DEL/DEmi or DEP/DEmi network were retained. ii) miRNA targets from miRNA interaction databases including miRDB²⁴ version 6.0, miRTarBase²⁵ release 8.0 and miRcode²⁶ version 11 were curated. Interactions with predicted score>80 were retained in miRDB database for confidence pursuance. Considering the incompatibility of gene identifiers in the multi-source interactions, R packages miRBaseVersions.db and org.Hs.eg.db (https://bioconductor.org/packages/release/data/annotation/html/org.Hs.eg.db.html) further aided the gene identifier conversion. Only DEM/DEmi interactions shared by miRDB and miRTarBase were used for further network construction. iii) For each TP53 mutational status-specific DEmi, the DEM targets in the shared interaction set were firstly retrieved and such interactions were intercepted with DEM/DEmi co-expression network, forming DEM/DEmi pairs in the ceRNA network. Analogously, the DEL targets of these DEmis were obtained from miRcode and intercepted with DEL/DEmi co-expression networks. When incorporating the mRNA-miRNA-lncRNA triplet in final ceRNA networks, only pairs existed in DEM/DEL co-expression networks were maintained. The Cytoscape²⁷ package was applied for ceRNA network visualization and the expressional variation of network nodes was illustrated by pheatmap R package (https://cran.r-project.org/web/packages/pheatmap/index.html).

Key prognosis-related node identification in ceRNA network

Survival analysis was conducted for the prioritization of key prognostic nodes in constructed ceRNA networks. To begin with, the normalized and log2-transformed counts of genes were integrated with survival information. Later univariate cox regression was conducted using survival R package²⁸ to examine the associations between gene expression levels and Overall Survival (OS) and Disease-Specific Survival (DSS). Significant nodes with p-value<0.05 were considered to be prognosis-related and were further subjected to Kaplan-Meier (K-M) curve generation by survminer package²⁹. A p-value<0.05 in log-rank test was considered with statistical difference.

Immune infiltration analyses for ceRNA network nodes

To detect potential associations between ceRNA network prognosis-related nodes and immune cell infiltration levels, the immune infiltration analyses were conducted on curated tumor samples. More specifically, gene-level FPKM values were downloaded from Xena database³⁰ and Ensembl gene IDs were converted to gene symbols using in-house scripts. Later R package Immunedeconv³¹ was used for the deconvolution of gene expressions and the derivation of the immune cell composition in tumor microenvironment. The xCell method³² implemented in Immunedeconv package was selected regarding its capability of portraying the enrichment landscape of the most immune and stroma cell types. Activities of cancer-immunity cycle³³ steps were addedly assessed by single sample Gene Set Enrichment Analysis (ssGSEA) from GSVA R package³⁴. Associations between nodes in ceRNA networks and immune infiltration analysis results as well as GSVA enrichment scores were further assessed by Spearman’s rank correlation.

The deduction of potential chemotherapy resistance for ceRNA network nodes

To elucidate the association between ceRNA network nodes and possible chemotherapy resistance in HCC with different TP53 mutational statuses, chemotherapy resistance-related genes were firstly collected from literatures^35,36. More detailedly, genes related to i) tumor microenvironment (TME), including TME remodeling, cytokine and cancer stem cell; ii) liver cancer specific events, including hepatitis B virus infection, HCC-related protein 1 (HCRP1) and HCC target therapy; iii) specific pathways and mechanisms, including transmembrane transporter, DNA damage and repair, telomerase, lipid metabolism, head shock protein, autophagy, apoptosis protein and chemotherapy critical pathway; iv) transcriptional and translational regulators, including topoisomerase, epigenetics and translational factor were incorporated in the analyses. Associations between the four gene categories and ceRNA network nodes were further quantified by Spearman’s rank correlation coefficients (SCC).

Methylation analysis on key mRNAs in ceRNA network

To explicate the expressional alternations of ceRNA network nodes, Illumina 450k methylation data was downloaded from Xena database. Experimental variations were further eliminated by discarding i) probes with not available (NA) value in more than 5% of samples; ii) samples with NA value in more than 20% of probes. Later BMIQ algorithm implemented in wateRmelon package³⁷ was used for the bias correction of probe-level methylation values and the normalized methylation intensities of ceRNA network nodes were scrutinized in samples with different TP53 mutational statuses. Finally, we conducted expressional comparisons on common 5-methylcytosine (5mC) methyltransferase (writer), demethylase (eraser) and methylation recognition (reader) proteins³⁸ in TP53 MT, TP53 WT and normal sample groups. Wilcoxon rank sum test was used in all group-wise comparisons.

Drug response difference investigation in patients with distinct TP53 mutational statuses

Data of drug response was download from the Genomics of Drug Sensitivity in Cancer database (GDSC, https://www.cancerrxgene.org/)³⁹, which contains drug sensitivity data measured by fluorescence-based cell viability assays on 805 cell lines subjected to 72 hour treatment of 198 drugs. The sensitivity data was used to construct a ridge regression model by oncoPredict R package⁴⁰, which further aided the imputation of drug response values on HCC expressional data. Group-specific drug candidates were selected by comparing the drug response differences between TP53 MT and TP53 WT groups. A p-value<0.05 in Wilcoxon rank sum test was considered significant.

Multi-omics validation on the expression levels of ceRNA network components

Multi-omics HCC datasets were collected from International Cancer Genome Consortium (ICGC, https://dcc.icgc.org/)⁴¹ and Cancer Therapeutics Response Portal (CTRP, https://portals.broadinstitute.org/ctrp)⁴² databases for ceRNA network component validation. After removing samples possessing TP53 synonymous mutations, RNA expression levels of 84 TP53 MT and 202 normal samples from ICGC database were used for verification. Similarly, the difference was measured by Wilcoxon rank sum test. As for the CTRP database, protein expression levels of ceRNA network key node in 9 TP53 MT and 5 TP53 WT samples were selected for comparison. In addition, immunohistochemistry data was retrieved from the Human Protein Atlas (HPA, https://www.proteinatlas.org/) database to validate the protein levels of key network nodes in normal and tumor liver tissues.

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 patients⁴³, 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.

Liver hepatocellular carcinoma (HCC), one of the most common and mortal cancer type, has attracted pullulating attention to mechanism exploration and therapy development. Owing to the perturbative effect and usualness of TP53 mutational abnormality in HCC, a comprehensive investigation on the altered microenvironment and its connection with systemic therapy could assist the development of novel medicinal strategies. In this study, HCC tumor samples from TCGA database were primarily stratified by TP53 mutational status. Based on the DEMs, DELs and DEmis identified from groups, TP53 mutant-specific ceRNA network was constructed. After confirming the prognostic value of key network node genes, therapeutic correlation analyses on node expression and immune cell infiltrations, cancer-immunity cycle activity, chemotherapy resistance were conducted. Group-wise drug sensitivity as well as the DNA methylome of node genes was also scrutinized. Finally, generality of the identified gene target was tested by public databases. To our knowledge, this is the first study focusing on the ceRNA network mediated by TP53 mutations and its therapeutic associations in HCC, which could undoubtedly complement the gap in precision medicine.

Basing on the premise that patients stratified by TP53 mutational status exhibited significantly different survival and the functional diversity of the gene, we utilized ceRNA network to unveil the genetic alternations as well as pathways provoked by TP53 abnormality. Most of the previous applications of ceRNA network focused on the discrepancy regulators between tumor and normal samples^44,45. However, the intrinsic heterogeneity of tumor often engenders homogeneous cancer subtypes possessing distinct molecular features and therapeutic responses, which could invalidate most of the biomarkers from previous studies. Emerging studies have incorporated tumor heterogeneity in ceRNA network analyses^46,47 but works considering key driver mutations were scarce. Interestingly, the DEGs identified in patients separated by TP53 mutational status were enriched in similar pathways, denoting that the effect conferred by TP53 anomaly should be investigated in a more microcosmic way, e.g. through the construction of regulatory networks. Noticing our ceRNA network only existed in TP53 mutated patients, the specificity, veracity and efficacy of the biomarkers in network was emphasized.

Multiple gene targets including hsa-miR-133b, AC017104.6, RP5-1092A11.5, RP11-365O16.6 and PTBP1 were prioritized in our ceRNA network construction procedure, some of which were reported to be involved in cancer-specific ceRNA networks. For example, hsa-miR-133b was identified as a prognostic marker in esophageal cancer⁴⁸ and cervical cancer⁴⁹. Our results italicized its role in TP53-mutant HCC homeostasis maintenance. As for the three lncRNA nodes, we uniquely reported their prognostic value as well as the associations with immune infiltration and chemotherapy resistance in HCC, again supplemented the therapeutic gene target set. Functioning as the only mRNA in our network, PTBP1 demonstrated the strongest correlation with immune cell infiltrations and chemoresistance. With the ability to control almost all steps of messenger RNA metabolism and processing⁵⁰, the elevation of PTBP1 caused an adverse effect in glioblastoma⁵¹ and colorectal cancer⁵². In HCC, relationships between PTBP1 up-regulation and tumor growth⁵⁰, tumor invasion as well as metastasis⁵³ were uncovered. Our work for the first time proved the detriment of PTBP1 on TP53-mutant patient survival. Considering the relationship between expressions of PTBP1 and TP53 in glioma⁵⁴, we speculate that the TP53 mutational deactivation in HCC triggered the higher expression of PTBP1, which quantity was maintained by our identified ceRNA network and acted as the accomplice of mutated TP53 in tumor development.

In the immune-related analyses, infiltration level of common lymphoid progenitor cells was found closely associated with expressions of ceRNA network nodes but the cell prevalence disparity between TP53 MT and TP53 WT group did not reach statistical significance. Since common lymphoid progenitor cells give rise to multiple types of lymphocytes, we suspect the key pathways’ activity as well as the infiltration of some differentiated lymphocyte in the lymphoid lineage could be different between TP53 MT and WT group. Indeed, the infiltrations of T cell CD4+ Th2, T cell CD4+ memory and mast cell were statistically different in the two groups. For the T cell CD4+ Th2 and CD4+ memory T cells, they demonstrated higher infiltration levels in TP53 MT patients. A recent study reported that Th2 cells initiate an anti-tumor response through the MHC II complex pathway⁵⁵ but the elimination of tumors is associated with CD8+ T cells. Regarding the observation on impaired PD-1 related pathway in HCC, adoptive CD8+ T cell therapy (ACT) could possibly benefit TP53 MT HCC patients. Besides, CD4+ T memory cells could instigate secondary immune responses by influencing CD8+ T cell response early events⁵⁶, again demonstrated the potential of applying ACT in specific HCC groups. Moreover, the two-sidedness of mast cells in tumor development was previously reported⁵⁷. An elevated infiltration of mast cells was found in our TP53 WT HCC samples while higher mast cell prevalence correlated with worse prognosis in TP53 MT group. We suspect TP53 mutations triggered such cellular functional switch and approaches targeting mast cells could be a new therapeutic option for HCC.

In our multi-faceted analyses on systematic therapy, chemotherapy resistance mechanisms were found closely associated with ceRNA network nodes, among which the key mRNA PTBP1 demonstrated the strongest correlation. Aside from prioritizing it as a novel inhibition therapy target, we also proposed possible directions of epigenetic therapeutics. Moreover, drugs including AZD8055, OSI.027 and RO.3306 exhibited higher efficacy in TP53-mutated HCC patients. All of them hamper the cell proliferation process. Considering the unrestricted cell cycle granted by TP53 mutational abnormality, these drugs could function as therapeutic surrogates in the way of conquering liver cancer harboring TP53 mutations.

However, our study still has some limitations. Firstly, only data from TCGA database was used in our TP53 mutant-specific ceRNA network construction procedure, which may cause possible racial bias on the results. With the advent of multi-omics sequencing data in multiracial liver cancer, these potential partialities could be eliminated. Secondly, relatively rigorous RNA inclusion and exclusion criteria were applied in our analyses, which could nurture our network’s specificity but limit the number of gene targets identified. Attempts pursuing a perfect balance between network breadth and depth should be made in the future studies. Thirdly, we confirmed the prognosis relativity of the key ceRNA network nodes in TCGA data. Further investigations incorporating multi-source samples with outcome information could better consolidate our discoveries. Moreover, our conclusions on immunotherapy and chemotherapy were mainly derived from correlation analyses. With the availability of future data measuring the transcriptomic dynamics conferred by therapies, the underlying association between our ceRNA network and precision medicine could be further elucidated in HCC. Last but not least, most of our conclusions were sorely data-based. Though multiple databases were involved in the validation of the network nodes, using experiments including quantitative polymerase chain reaction (qPCR), western blot as well as immunofluorescence (IF) to validate the expressions of key network members in cell lines or curated clinical samples could undoubtedly solidify our findings. Besides, techniques including luciferase assay and cross-linked immuno-precipitation could also aid the unveiling of biological mechanisms amongst the network nodes. By applying experiments in the subsequent studies, the role of our ceRNA network in TME modulation could be further illustrated.

In conclusion, the present study identified a ceRNA network mediated by TP53 mutations and focused on the therapeutic alternations associated with the network, providing unprecedented prognostic and diagnostic markers for HCC. Additionally, our comprehensive analyses complemented the understanding of systemic therapy in HCC, which could undoubtedly promote the development of precision medicine in HCC.

Data Availability Statement

The datasets analyzed during the current study are available in the TCGA GDC Data Portal (https://portal.gdc.cancer.gov/) and cBioPortal database (https://www.cbioportal.org/).

Acknowledgements

Not applicable.

Author Contributions

JY contributed to the study design. ZX, HG, RL, JJ and JY analyzed the data. ZX drafted the manuscript. JY and JP revised the manuscript. All authors read and approved the final version of the manuscript.

Funding

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

Ethical Approval

Not applicable.

Competing Interests

The authors declare that they have no competing interests.

1. Ferlay, J. et al. Cancer incidence and mortality worldwide: Sources, methods and major patterns in GLOBOCAN 2012. Int. J. Cancer 136, E359–E386 (2015).

2. Zheng, R. et al. Liver cancer incidence and mortality in China: temporal trends and projections to 2030. Chin. J. Cancer Res. 30, 571 (2018).

3. Chidambaranathan-Reghupaty, S., Fisher, P. B. & Sarkar, D. Hepatocellular carcinoma (HCC): Epidemiology, etiology and molecular classification. in Advances in Cancer Research vol. 149 1–61 (Elsevier, 2021).

4. Yang, J. D. et al. A global view of hepatocellular carcinoma: trends, risk, prevention and management. Nat. Rev. Gastroenterol. Hepatol. 16, 589–604 (2019).

5. Huang, A., Yang, X.-R., Chung, W.-Y., Dennison, A. R. & Zhou, J. Targeted therapy for hepatocellular carcinoma. Signal Transduct. Target. Ther. 5, 146 (2020).

6. Sengupta, S. & Harris, C. C. p53: traffic cop at the crossroads of DNA repair and recombination. Nat. Rev. Mol. Cell Biol. 6, 44–55 (2005).

7. Kamp, W. M., Wang, P. & Hwang, P. M. TP53 mutation, mitochondria and cancer. Curr. Opin. Genet. Dev. 38, 16–22 (2016).

8. Stiewe, T. & Haran, T. E. How mutations shape p53 interactions with the genome to promote tumorigenesis and drug resistance. Drug Resist. Updat. 38, 27–43 (2018).

9. Zucman-Rossi, J., Villanueva, A., Nault, J.-C. & Llovet, J. M. Genetic Landscape and Biomarkers of Hepatocellular Carcinoma. Gastroenterology 149, 1226-1239.e4 (2015).

10. Liu, J. et al. Alterations of TP53 are associated with a poor outcome for patients with hepatocellular carcinoma: Evidence from a systematic review and meta-analysis. Eur. J. Cancer 48, 2328–2338 (2012).

11. Hu, J. et al. Targeting mutant p53 for cancer therapy: direct and indirect strategies. J. Hematol. Oncol.J Hematol Oncol 14, 157 (2021).

12. Akula, S. M. et al. RAS/RAF/MEK/ERK, PI3K/PTEN/AKT/mTORC1 and TP53 pathways and regulatory miRs as therapeutic targets in hepatocellular carcinoma. Expert Opin. Ther. Targets 23, 915–929 (2019).

13. Chen, R. et al. Comprehensive analysis of an immune-related ceRNA network in identifying a novel lncRNA signature as a prognostic biomarker for hepatocellular carcinoma. Aging 13, 17607–17628 (2021).

14. Kozomara, A., Birgaoanu, M. & Griffiths-Jones, S. miRBase: from microRNA sequences to function. Nucleic Acids Res. 47, D155–D162 (2019).

15. Colaprico, A. et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 44, e71 (2016).

16. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).

17. Smyth, G. K. limma: Linear Models for Microarray Data. in Bioinformatics and Computational Biology Solutions Using R and Bioconductor (eds. Gentleman, R., Carey, V. J., Huber, W., Irizarry, R. A. & Dudoit, S.) 397–420 (Springer, 2005). doi:10.1007/0-387-29362-0_23.

18. Robinson, M. D., McCarthy, D. J. & Smyth, G. K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinforma. Oxf. Engl. 26, 139–140 (2010).

19. Chen, H. & Boutros, P. C. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics 12, 35 (2011).

20. Chen, E. Y. et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics 14, 128 (2013).

21. Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y. & Morishima, K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 45, D353–D361 (2017).

22. Ashburner, M. et al. Gene Ontology: tool for the unification of biology. Nat. Genet. 25, 25–29 (2000).

23. Jassal, B. et al. The reactome pathway knowledgebase. Nucleic Acids Res. 48, D498–D503 (2020).

24. Chen, Y. & Wang, X. miRDB: an online database for prediction of functional microRNA targets. Nucleic Acids Res. 48, D127–D131 (2020).

25. Huang, H.-Y. et al. miRTarBase 2020: updates to the experimentally validated microRNA-target interaction database. Nucleic Acids Res. 48, D148–D154 (2020).

26. Jeggari, A., Marks, D. S. & Larsson, E. miRcode: a map of putative microRNA target sites in the long non-coding transcriptome. Bioinforma. Oxf. Engl. 28, 2062–2063 (2012).

27. P, S. et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, (2003).

28. Therneau, T. M. & Lumley, T. Package ‘survival’. Surviv. Anal. Publ. CRAN 2, 119 (2014).

29. Kassambara, A., Kosinski, M., Biecek, P. & Fabian, S. Package ‘survminer’. Draw. Surviv. Curves Using ‘ggplot2’R Package Version 03 1 (2017).

30. Goldman, M. J. et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat. Biotechnol. 38, 675–678 (2020).

31. Sturm, G. et al. Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology. Bioinformatics 35, i436–i445 (2019).

32. Aran, D., Hu, Z. & Butte, A. J. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 18, 220 (2017).

33. Chen, D. S. & Mellman, I. Oncology meets immunology: the cancer-immunity cycle. Immunity 39, 1–10 (2013).

34. Hänzelmann, S., Castelo, R. & Guinney, J. GSVA: gene set variation analysis for microarray and RNA-Seq data. BMC Bioinformatics 14, 7 (2013).

35. Guo, J. et al. Mechanisms of resistance to chemotherapy and radiotherapy in hepatocellular carcinoma. Transl. Cancer Res. 7, (2018).

36. Lohitesh, K., Chowdhury, R. & Mukherjee, S. Resistance a major hindrance to chemotherapy in hepatocellular carcinoma: an insight. Cancer Cell Int. 18, 44 (2018).

37. Teschendorff, A. E. et al. A beta-mixture quantile normalization method for correcting probe design bias in Illumina Infinium 450 k DNA methylation data. Bioinforma. Oxf. Engl. 29, 189–196 (2013).

38. Jiang, S. & Guo, Y. Epigenetic Clock: DNA Methylation in Aging. Stem Cells Int. 2020, 1047896 (2020).

39. Iorio, F. et al. A Landscape of Pharmacogenomic Interactions in Cancer. Cell 166, 740–754 (2016).

40. Maeser, D., Gruener, R. F. & Huang, R. S. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief. Bioinform. 22, (2021).

41. Zhang, J. et al. The International Cancer Genome Consortium Data Portal. Nat. Biotechnol. 37, 367–369 (2019).

42. Basu, A. et al. An interactive resource to identify cancer genetic and lineage dependencies targeted by small molecules. Cell 154, 1151–1161 (2013).

43. Macek Jilkova, Z., Aspord, C. & Decaens, T. Predictive Factors for Response to PD-1/PD-L1 Checkpoint Inhibition in the Field of Hepatocellular Carcinoma: Current Status and Challenges. Cancers 11, 1554 (2019).

44. Jiang, J. et al. To construct a ceRNA regulatory network as prognostic biomarkers for bladder cancer. J. Cell. Mol. Med. 24, 5375–5386 (2020).

45. Tang, F. et al. Competitive endogenous RNA (ceRNA) regulation network of lncRNAs, miRNAs, and mRNAs in Wilms tumour. BMC Med. Genomics 12, 194 (2019).

46. Zhang, K. et al. A ceRNA network and a potential regulatory axis in gastric cancer with different degrees of immune cell infiltration. Cancer Sci. 111, 4041–4050 (2020).

47. Huang, R. et al. The construction and analysis of tumor-infiltrating immune cells and ceRNA networks in metastatic adrenal cortical carcinoma. Biosci. Rep. 40, BSR20200049 (2020).

48. Zhang, H.-C. & Tang, K.-F. Clinical value of integrated-signature miRNAs in esophageal cancer. Cancer Med. 6, 1893–1903 (2017).

49. Qin, S. et al. Identifying Molecular Markers of Cervical Cancer Based on Competing Endogenous RNA Network Analysis. Gynecol. Obstet. Invest. 84, 350–359 (2019).

50. Kang, H. et al. A miR‐194/PTBP1/CCND3 axis regulates tumor growth in human hepatocellular carcinoma. J. Pathol. 249, 395–408 (2019).

51. Kim, J.-H. et al. SON drives oncogenic RNA splicing in glioblastoma by regulating PTBP1/PTBP2 switching and RBFOX2 activity. Nat. Commun. 12, 5551 (2021).

52. Chen, J. et al. Circular RNA circRHOBTB3 represses metastasis by regulating the HuR-mediated mRNA stability of PTBP1 in colorectal cancer. Theranostics 11, 7507–7526 (2021).

53. Shen, L. et al. Skipping of exon 10 in Axl pre-mRNA regulated by PTBP1 mediates invasion and metastasis process of liver cancer cells. Theranostics 10, 5719–5735 (2020).

54. Sheng, J. et al. p53-targeted lncRNA ST7-AS1 acts as a tumour suppressor by interacting with PTBP1 to suppress the Wnt/β-catenin signalling pathway in glioma. Cancer Lett. 503, 54–68 (2021).

55. Schreiber, S. et al. Metabolic Interdependency of Th2 Cell-Mediated Type 2 Immunity and the Tumor Microenvironment. Front. Immunol. 12, 632581 (2021).

56. Krawczyk, C. M., Shen, H. & Pearce, E. J. Memory CD4 T Cells Enhance Primary CD8 T-Cell Responses. Infect. Immun. 75, 3556–3560 (2007).

57. Ribatti, D., Tamma, R. & Crivellato, E. The dual role of mast cells in tumor fate. Cancer Lett. 433, 252–258 (2018).

There is NO conflict of interest to disclose.

Supplementarymaterials.pdf
Supplementary Figures

PDF

Version 1

posted 17 Mar, 2022

You are reading this older preprint version

Read the latest preprint version →

Comprehensive analysis of a ceRNA network associated with TP53 mutations in liver hepatocellular carcinoma

Status:

Version 1

Abstract

Figures

Introduction

Materials And Methods

Results

Discussion

Declarations

References

Competing Interests

Supplementary Files

Status:

Version 1