Skip to main content
Download PDF

Potential prognosis index for m6A-related mRNA in cholangiocarcinoma

BMC Cancer volume 22, Article number: 620 (2022) Cite this article

Abstract

Background

Cholangiocarcinoma (CHOL) is a malignant tumor that originates in the extrahepatic bile duct and can extend from the hilar region to the lower end of the common bile duct. The prognosis of CHOL patients is particularly poor; therefore, in this study, we screened mRNAs correlated with N6-methyladenosine (m6A) to construct a risk model for prognosis in CHOL.

Methods

The TCGA-CHOL dataset was applied to obtain and analyze the coexpression of 1281 m6A-related mRNAs, from which 14 were selected for further analysis through univariate proportional hazards (cox) regression analysis. Aryl hydrocarbon receptor interacting protein (AIP), CCAAT/enhancer binding protein beta (CEBPB), syndecan1 (SDC1), vacuolar protein sorting 25 homolog (VPS25) and syntaxin binding protein 2 (STXBP2) were then screened out through the least absolute shrinkage and selection operator (LASSO) and multivariate Cox regression analysis to develop a precise m6A-related mRNA prognosis risk model (MRMRPM) with an area under curve (AUC) of 0.908 and 0.923 after 1 and 2 years, respectively. We divided the samples into high-risk and low-risk groups using the m6A-related mRNA prognosis risk model.

Results

Kaplan–Meier analysis indicated poor overall survival (OS) for the high-risk group. Two Gene Expression Omnibus (GEO) datasets (GSE89748 and GSE107943) were used to validate the risk model. The results of drug sensitivity and immune cell infiltration analysis showed that the risk model could serve as a prognosis index of potential immunotherapeutic characteristics and drug sensitivity. Furthermore, the proportion of resting dendritic cells and regulatory T cells was positively associated with an increased expression of four m6A-related mRNAs — AIP, CEBPB, SDC1, and VPS25 — in the high-risk CHOL group.

Conclusions

Our findings suggest that this model can be a prognostic indicator for CHOL patients.

Peer Review reports

Background

Cholangiocarcinoma (CHOL), which originates from epithelial cells in the bile duct, is anatomically divided into intrahepatic and extrahepatic cholangiocarcinoma [1]. The morbidity and mortality of CHOL continues to increase and accounts for 7–10% of all hepatobiliary malignancies [2]. The lack of specific early diagnostic markers and molecular targets for therapy means that most CHOL patients are diagnosed at an advanced stage and radical surgical resection can treat only 20% of patients [3]. These factors indicate that the 2-year OS is only 5.5% today [4]. The recent increase in research on cholangiocarcinoma has led to the discovery and elucidation of the expression characteristics and mechanisms of some key molecules involved in the development of cholangiocarcinoma, providing a new basis for targeted molecular therapies.

m6A is a widespread modification of the mRNA base that involves methylation of the sixth N position and generally affects adenosine (A). This modification maintains the stability of mRNA and is directly involved in the migration and proliferation of tumor cells [5,6,7]. The importance of this modification in gene regulation has generated considerable research interest in m6A.

Several studies have indicated that m6A genes regulate the formation and development of tumors. The high expression of methyltransferase-like proteins (METTL)3 in hepatocellular carcinoma (HCC) induces degradation of the tumor suppressor gene suppressor of cytokine signaling 2 (SOCS2) via m6A modification and is associated with poor prognosis [8]. METTL3 is also highly expressed in bladder cancer, and mediated target genes form a multilevel regulatory network that promotes tumor growth [9]. In recent years, several studies have suggested that the abnormal regulation of the m6A genes is related to CHOL [10, 11]. As the most common mRNA modification, transcriptome studies had linked m6A to cancer development, affecting the self-renewal and differentiation, proliferation, apoptosis, invasion and metastasis, resistance, and immune suppression processes of tumor stem cells [12,13,14,15]. Therefore, m6A-related mRNA is a potential molecular target for cancer treatment and drug development. However, the correlation between m6A-related mRNA, immune infiltration, and clinical prognosis in CHOL patients remains unclear.

In our study, the mRNA expression profiles, including expression data for 21 m6A genes, related to 36 CHOL samples were obtained from The Cancer Genome Atlas (TCGA) database. mRNAs related to the m6A genes were first screened based on Pearson’s correlation analysis. The m6A-related mRNA expression and prognostic information of the 36 TCGA-CHOL cohorts were then applied to establish an individual prognostic risk model for CHOL and the prognostic risk model was validated using two GEO datasets. We then analyzed the tumor mutational burden (TMB) of each TCGA-CHOL cohort. The sensitivity of 83 candidate drugs was explored using external cell lines from the Cancer Cell Line Encyclopedia (CCLE) to obtain drug targets for a risk model based on the relationship between the five obtained mRNAs and the m6A genes. We found that the relationship between the immunotherapy response and m6A-related mRNAs could promote infiltration by immune cells in CHOL. Thereafter, a model was constructed using the five m6A-related mRNAs to predict the OS of CHOL patients through nomogram analysis.

Methods

Database of CHOL patients

All data used to construct the m6A-related mRNA risk model were obtained from the TCGA database (https://portal.gdc.cancer.gov/). Testing sets were downloaded from the GEO database on the GEO website (http://www.ncbi.nlm.nih.gov/geo).

Identification of m6A genes and m6A-related mRNAs

The expression database of 21 m6A genes and all protein-coding mRNAs were extracted from the TCGA database. m6A-related mRNAs were selected using the Pearson correlation coefficient method to calculate correlations between m6A and other genes (R > 0.3 and p < 0.001). A total of 1281 m6A-related mRNAs were screened.

Analysis of function and protein-protein interaction

The functions of m6A-related mRNAs with different expression were determined by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyzes with p-values < 0.05 [16,17,18]. The online database STRING (https://string-db.org/) was used to construct the PPI network.

Identification and validation of the MRMRPM model

The m6A-related mRNA model was developed utilizing 36 TCGA-CHOL samples for training and validated using the GSE89748 and GSE107943 datasets that were downloaded from GEO. We screened m6A-related mRNA from 1281 m6A-related mRNAs (p < 0.05) via univariate Cox regression [19] and used LASSO with 10-fold cross validation to find seven m6A-related mRNAs [20]. The five-m6A-related mRNA risk model was then established through multifactor cox regression. The formula MRMRPM = Ʃ (βi × EXPi) was used to calculate the risk score based on multivariate Cox analysis.

Analysis of candidate target genes for transcription factors

The candidata target genes of each transcription factor were extracted according to cotarget genes from the ChEA ChIP-X target gene database (https://maayanlab.cloud/Harmonizome/dataset/CHEA+Transcription+Factor+Targets) [21].

Principal component analysis

PCA was used to perform model recognition, dimensionality reduction, and visualization of the groups associated with the whole gene expression, m6A gene expression, and the expression of five m6A-related mRNA to further verify the discriminative ability of the model [22].

Kaplan–Meier survival analysis

Kaplan–Meier analysis was performed using the “survival” package. We did not use any value to exclude the survival time, because the sample was too small.

Drug sensitivity prediction targeting for the m6A-related mRNA model in clinical treatment

The CCLE gene expression data were used to predict the drug sensitivity for MRMRPM and the sensitivity of patients to 83 drugs embedded in the pRRophetic package was predicted by the gene expression in the CHOL samples. The difference in the prognostic groups predicted by the model in terms of their sensitivity to drugs was then analyzed, combined with the two risk groups.

Investigation of the m6A-related mRNA model regarding immunotherapeutic treatment

First, we counted the TMB in each TCGA-CHOL sample to produce the total amount of coding variants/the length of exons (38 million), and we then divided the 36 samples into high- and low-TMB groups. Second, we evaluated the abundance of immune cells in the 36 TCGA-CHOL samples using CIBERSORT (https://cibersort.stanford.edu) with p < 0.05. We then analyzed the correlation between m6A-related mRNA expression and the infiltration by immune cells.

Establishing and validating the predictive nomogram

Comparing other clinical characteristics (clinical stage, age, sex, and risk value), we validated MRMRPM by multivariate and univariate Cox regression analyses of the TCGA-CHOL patients.

Furthermore, a nomogram for MRMRPM was established using 1- and 2-year OS, combined with other clinical indicators (clinical stage, age, sex, and risk value).

Statistical analysis

The R software clusterProfiler package and the survival ROC R software package were used to perform gene functional enrichment analyses and the ROC curve, respectively. The Student’s t-test was used to assess significant differences.

Results

Profiles of m6A-related mRNAs in CHOL

The research process is presented in Fig. 1. The profiles of 21 m6A gene expressions and the mRNA expression in 36 samples were downloaded from the TCGA-CHOL dataset (Table 1). Data from 16 (44%) males and 20 (56%) females were included in the study. There are no differences between male and female in the sample, and indicated that there is no specific bias in gender and tumor stage for our study (Table 1). m6A-related mRNAs associated with at least one of the 21 m6A genes were defined using the Pearson correlation coefficient (R > 0.35 and p < 0.001). The coexpression network obtained for m6A-mRNA is presented in Fig. 2A, and the 1281 identified m6A-related mRNAs are presented in Supplementary Table S1. GO analysis indicated that in terms of biological processes (Fig. 2B), RNA catabolic process, protein targeting, mRNA catabolic processes, the establishment of proteins on membranes, small molecule catabolic processes, ATP metabolic processes, and cellular amino acid metabolic processes were enriched in samples, whereas m6A-related mRNAs were found mainly in the mitochondrial protein complex, the inner mitochondrial inner membrane, inner mitochondrial membrane protein complex, mitochondrial matrix, ribosomal subunit, cytosolic ribosome, large ribosomal subunits, and mitochondrial ribosomes within cells (Fig. 2C). In terms of molecular function, m6A-related mRNAs were mainly involved in producing the structural constituents of ribosomes (Fig. 2D). Following GO analysis, we analyzed data from the Kyoto Encyclopedia of Genes and Genomes (KEGG) to show that m6A-related mRNA was also implicated in COVID-19 and amyotrophic lateral sclerosis. Therefore, the datasets obtained indicate that m6A-related mRNAs were mainly related to biological functions and metabolic pathways of ribosomes and mitochondria (Fig. 2E). We then analyzed PPI networks using the STRING database to produce the interaction network in Fig. 2F. Isolated points, namely m6A-related mRNAs without interactions, were removed from the network. The maximum connected subgraph in the network was then extracted using the igraph [23] package to obtain 57 hub genes.

Fig. 1
figure 1

Workflow chart for this study

Full size image
Table 1 Summary of patient characteristics
Full size table
Fig. 2
figure 2

Identification of m6A-related mRNAs in Cholangiocarcinoma patients and enrichment pathway analysis. A Sankey relational diagram of 21 m6A genes and m6A-related mRNAs from 36 samples in TCGA, and |Pearson R| > 0.35 and p < 0.001 analysis with the m6A-related mRNA significance filtering threshold. B-D Function analysis of GO enrichment in differentially expressed m6A-related mRNAs. E Immune-related pathways involved in the differential expression of m6A-related mRNAs by KEGG pathway analysis. F PPI network of 57 proteins. TCGA: The Cancer Genome Atlas; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; PPI: protein-protein interaction

Full size image

Construction prognostic model according to m6A-related mRNAs in CHOL

Univariate Cox regression analysis of m6A-related mRNAs was performed with the R package “survival.” Them6A-related mRNAs were considered associated with a prognosis and therefore applicable for subsequent analysis when p < 0.05. Fourteen m6A-related mRNA genes from the 1281 m6A-related mRNAs associated with CHOL were significantly related to OS (Table 2). LASSO was used to select a variable optimization model and LASSO Cox analysis was conducted using the significant m6A-related mRNAs to construct a risk prognosis model [20]. The 10-fold cross-validation method was used to eliminate collinearity in gene optimization and simplified models. The variable λ in LASSO was utilized to find the best model. As λ increased, the regression coefficient β for each variable decreased, and a λ of 0 indicates that a variable contributes little to the model and could be eliminated. The LASSO screening indicated that seven m6A-related mRNAs were suitable for constructing the risk prognosis model and were used in subsequent studies (Fig. 3A, B). LASSO analysis was followed by multivariate Cox analysis and five characteristic m6A-related mRNA prognosis biomarkers for CHOL patients were obtained (Fig. 3C, Table 3). The relative hazard ratio (HR) for the five m6A-related mRNAs are presented in a forest plot (Fig. 3C). The MRMRPM was constructed based on the expression of the five m6A-related mRNAs, from which the risk to each patient was obtained. The expression of the five m6A-related mRNAs was statistically related to clinical characteristics; neither SDC1 nor the CEBPB showed any such correlation. Subsequently, we found that AIP and SDC1 were target genes of the transcription factor CEBPB from the ChEA ChIP-X target gene database (Supplementary Table S2). We also analysed the function of the five genes (Fig. 3D, E), The result of GO enrichment showed that all of the five mRNAs were involved in protein binding (Fig. 3D), and methylation analysis indecated that normoal tissue has higher methylation level of AIP, the methylation level of CEPBP and SDC1 was increased in the normoal tissue, and there was no difference of methylation level of VPS25 and STXEBP2 (Fig. 3E).

Table 2 Univariate Cox analysis reveals 14 m6A-related mRNA results significantly associated with prognosis
Full size table
Fig. 3
figure 3

Risk model for patients with Cholangiocarcinoma based on m6A-related mRNAs. A LASSO regression coefficient and lambda correspondence diagram, with LASSO coefficient profiles of the characteristics against the log (λ). B LASSO regression, partial-likelihood deviance curve with log (λ). The figure above shows lambdas. Min (collimated dashed line on the left) and lambdas.1se (vertical dashed line on the right). The value above represents the number of features included in the model with the associated λ value. C) The patient risk ratio of the five m6A-related mRNAs is shown in the forest map. D) GO enrichment analysis of the five m6A-related mRNAs. E methylation analysis of the five m6A-related mRNAs. LASSO: least absolute shrinkage and selection operator; GO: Gene Ontology

Full size image
Table 3 Five m6A-related mRNAs used to construct the risk prognosis model
Full size table

CHOL patients were then divided into high-risk (n = 18) and low-risk groups (n = 18) using the median risk as a threshold. The distribution of the low-risk and high-risk groups is presented in Fig. 4A, and the interrelation between the risk score and the survival period of CHOL patients in the two risk groups is reflected along with the survival status in Fig. 4B. The relative expression levels of the five m6A-related mRNAs for each patient in each risk group are shown in Fig. 4C. The expression of four of the five m6A-related mRNAs (excluding STXBP2) was increased in the high-risk CHOL patients (Fig. 4C). The patterns observed in the m6A-related mRNA expression were in accordance with the HR values (Table 3). Kaplan–Meier analysis indicated that high-risk CHOL patients had worse OS than low-risk CHOL patients (Fig. 4D). Supplementary Fig. S1 shows that high expression of the four m6A mRNAs (except STXBP2) (Supplementary Fig. S1E) was associated with worse OS in CHOL patients. The curves for receiver operating characteristics (ROC) that predict 1 and 2-year OS indicate that the model has an accurate risk prediction, with AUCs of 0.908 and 0.923, respectively (Fig. 4E). The 1- and 2-year OS predictions were calibrated using nomogram analysis, and all calibrated curves were well fit (Fig. 4F, G).

Fig. 4
figure 4

Estimation, evaluation, and calibration in training set for the five m6A-related mRNAs risk prognostic model A, B Distribution of the risk score for the five m6A-related mRNAs in the risk model. B Dot plot of the different survival time patterns for each risk group. C The different expression of m6 A-related mRNAs in each patient are shown in the heatmap. D Kaplan–Meier analysis of CHOL patients indicated that patients in the low-risk group had better OS. E The 1-year (AUC = 0.928) and 2-year (AUC = 0.903) ROC curves show the superior accuracy prediction of five independent prognoses in the m6A-related mRNA model. F, G) Calibration curves for survival after 1 y (F) and 2 y (G using MRMRPM. CHOL: Cholangiocarcinoma; OS: overall survival; AUC: area under curve; ROC: receiver operating characteristics; MRMRPM; five m6A-related mRNAs risk prognostic model

Full size image

Testing the prognostic capability of the constructed model

Two GEO datasets; GSE89748 with 71 patients and GSE107943 with 22 patients, were used to validate the prognostic model in this study. The risk value for each CHOL patient in each dataset was calculated according to the uniform formula. CHOL patients in the GSE89748 (Fig. 5A-C) and GSE107943 datasets (Fig. 6A-C) were divided into high-risk and low-risk groups for the following analysis based on the expression of the five m6A-related mRNAs (Fig. 5A and 6A). The same expression pattern for the five m6A-related mRNAs was observed in the TCGA (Fig. 4C), GSE89748 (Fig. 5C), and GSE107943 (Fig. 6C) groups. Only STXBP2 was upregulated in the low-risk group. Survival curves based on the two datasets also showed that the OS for CHOL patients was similar to that observed in the TCGA dataset and that the higher risk CHOL patients were associated with a worse OS than the lower risk CHOL patients (p = 0.022, Fig. 5D and p = 0.047, Fig. 6D). To validate the ability of the five m6A-related mRNAs risk prognostic model (MRMRPM), the conformance index and the area under the ROC curve (AUC) were also evaluated to achieve a risk score for the two datasets. The ROC curve for 1-, 3-, and 5-year OS prediction in the GSE89748 dataset suggest that the model possesses predictive accuracy with AUC values of 0.601, 0.688 and 0.681, respectively (Fig. 5E). Furthermore, the results of the analysis of the OS ROC curve at 1, 3, and 5 years indicated that the MRMRPM could accurately predict the prognosis of CHOL patients (AUC = 0.898, 0.725 and 0.722, respectively, Fig. 6E). As predicted, the OS and the ROC curves differed between the two risk groups, suggesting that the level of risk predicted by the MRMRPM accurately described the prognosis of CHOL.

Fig. 5
figure 5

Validation of five m6A-related mRNAs risk prognostic model using the GEO testing dataset. A, B Distribution of the risk score of the m6A-related mRNA risk model in GSE89748. B Dot plot of different survival time patterns for each risk group in GSE89748. C The different m6 A-related mRNA expression in each patient in GSE89748 is shown in the heatmap. D Kaplan–Meier analysis of CHOL patients shows that low-risk groups in GSE89748 have better OS. E The 1-year (AUC = 0.601), 3-year (AUC = 0.688) and 5-year (AUC = 0.686) ROC curves show the accuracy in the prediction of the prognostic m6A-related mRNA model in GSE89748. GEO: Gene Expression Omnibus; GSE:gene set enrichment; CHOL: Cholangiocarcinoma; OS: overall survival; AUC: area under curve; ROC: receiver operating characteristics

Full size image
Fig. 6
figure 6

Validation of MRMRPM using the GEO testing dataset (GSE107943, n = 22). A, B Distribution of the risk score of five m6A-related mRNA risk models in GSE107943. B Dot plot of different survival time patterns for each risk group in GSE107943. C The different m6 A-related mRNA expression in each patient in GSE107943 is shown in the heatmap. D Kaplan–Meier analysis of CHOL patients shown that low-risk groups had better OS in GSE89748. E The 1-year (AUC = 0.898), 3-year (AUC = 0.725) and 5-year (AUC = 0.722) ROC curves show the accuracy in the prediction of the m6A-related mRNA model in GSE107943. MRMRPM; five m6A-related mRNAs risk prognostic model; GEO: Gene Expression Omnibus; GSE:gene set enrichment; CHOL: Cholangiocarcinoma; OS: overall survival; AUC: area under curve; ROC: receiver operating characteristics

Full size image

Verification of the five m6A-related mRNA risk model by principal component analysis (PCA)

The entire gene expression profile, 21 m6A genes, the expression data for the five m6A-related mRNAs, the risk prognosis model generated using the m6A-related mRNA expression data, and the discrimination ability of the model were further verified by PCA (Fig. 7A–C). PCA performed on the expression levels of the characteristic mRNA (Fig. 7C) showed that the two risk groups differed, whereas the high-risk and low-risk individuals were randomly distributed over the two GEO datasets. These results verify the effective differentiating capacity of the characteristic five m6A-related mRNAs in the prognostic risk model.

Fig. 7
figure 7

PCA between high-risk and low-risk groups. A 21 m6A genes. B entire gene expression profiles. C) 5 m6A-related mRNAs. PCA:principal component analysis

Full size image

Testing drug susceptibility in the two risk groups

After the model was constructed and the differentiation ability of the model was verified by PCA, drug sensitivity was predicted to explore the difference between the prognosis in each group of patients as predicted by both the prognostic model and the drug sensitivity of the individuals involved.

Gene expression data detailing the drug sensitivity of external cell lines (CCLE) [24] were used to predict drug susceptibility by estimating the sensitivity of patients to 83 drugs embedded in pRRophetic by genetic expression in CHOL samples (Supplementary Table S3). The difference between the groups of patients predicted by the prognostic model and their sensitivity to drugs was then analyzed with the high-low-risk groups. The low-risk TCGA-CHOL cohorts showed more sensitivity to the 83 drug compounds. Supplementary Table S3 and Fig. 8A-H show the results for c-Jun N-terminal kinase, Bleomycin, Doxorubicin, EHT.1864, Elesclomol, FH535, Gefitinib, Imatinib. The results suggest that 9 L (LNK. 9 L) might be used for patients with CHOL (p < 0.05, Supplementary Table S3 and Fig. 8A).

Fig. 8
figure 8

Drug sensitivity analysis for the m6A-related mRNA model. A-H Diagram of different sensitivities in the high- and low-risk groups for JNK.9 L (A), Bleomycin (B), Doxorubicin (C), EHT.1864 (D), Elesclomol (E), FH535 (F), Gefitinib (G) and Imatinib (H)

Full size image

TMB analysis of the two risk groups

The somatic mutations of 36 CHOL patients were downloaded from the TCGA database and the mutations in the two risk groups analyzed (Fig. 9). Missense mutations were most frequently observed in both risk groups (Fig. 9A, B). The genes BAP1 (24%), PBRM1 (24%), CHD7 (18%), IDH (18%), MUC16 (18%), ADAM30 (12%), ANGEL (12%), APPL1 (12%), ARID1A (12%), and ARMC9 (12%) (Fig. 9A) had the highest number of mutations in the low-risk group, whereas ARID1A (22%), PBRM1 (22%), AHNAK (17%), ALB (17%), BAP1 (17%), BCOR (17%), BPHA2 (17%), KMT2C (17%), SRCAP (17%), and SRCAP (11%) had the highest number of mutations in the high-risk group (Fig. 9B). Meanwhile mutation in one of ARID1A and PBRM1 were found in almost half of 32 intrahepatic cholangiocarcinomas through exomic sequencing [25]. ARID1A and PBRM1 are involved in the formation of the BAF complex, and DNA damage in cells lacking PBRM1 leads to cellular dynamic chromatin instability [26, 27]. Whether in high- and low- risk, the frequent alterations in ARID1A and PBRM1 in cholangiocarcinomas highlight the key role of chromatin remodeling which in cholangiocarcinomas.

Fig. 9
figure 9

Statistics describing mutation and TMB in TCGA-CHOL samples. A Mutation information for the low-risk group is presented in the waterfall plot displays. B Waterfall plot displays mutation information for the high-risk group. C TMB differences observed in high- and low-risk patients. D Kaplan–Meier survival of TMB in CHOL patients. TMB: tumor mutational burden; TCGA: The Cancer Genome Atlas; CHOL: Cholangiocarcinoma

Full size image

The TMB per million bases was then calculated for each of the 36 TCGA-CHOL patients; however, no significant differences were observed between the two groups in terms of this factor (Fig. 9C, Table 4) and no correlation was observed between the m6A-based classification index and the TMB. The survival curves suggest that the TMB of CHOL patients is not associated with OS (p = 0.51) (Fig. 9D).

Table 4 Results of the t-test and Wilcoxon analysis of the tumor mutation burden between the high-and low-risk groups in TCGA-CHOL cohort
Full size table

To validate the sensitivity of MRMRPM to the immunotherapy response in CHOL, we next analyzed the relationship between the tumor immune microenvironment (IME) and the m6A-related mRNA model.

Sensitivity of the m6A-related mRNA risk model to the immunotherapy response

To further confirm the correlation between m6A-related mRNA expression and tumor IME, CIBERSORT software [28] was used to analyze the infiltration profiles of 22 types of tumor-infiltrating immune cells (TICs) in the high-risk and low-risk groups. The fractions of the 22 immune cells were obtained by filtering with the “CIBERSORT” package, with P > 0.05 for each of the 36 TCGA-CHOL patients (Fig. 10A).

Fig. 10
figure 10

Relationship between abundance of TICs and five m6A-related mRNA expression in the two risk groups. A Violin plot shows the proportion of TICs in the high- and low-risk groups for CHOL. B Scatter plot shows the correlation between the 22 types of TICs and AIP expression. C Scatter plot shows the correlation between 22 types of TICs and SDC1 expression. TICs: tumor-infiltrating immune cells; CHOL: Cholangiocarcinoma; AIP: aryl hydrocarbon receptor interacting protein; SDC1: syndecan1

Full size image

Moreover, our bioinformation analysis verified that the low-risk group had higher T cell regulatory (Tregs) infiltration levels (p = 0.023) than the high-risk group (Fig. 10A, Supplementary Table S4) and 22 kinds of TICs were associated with the expression of the five m6A-related mRNAs. Of the TICs, four types were positively associated with AIP expression, which is associated with memory B-cells, naïve B-cells, plasma cells, and Tregs. No TICs were negatively associated with AIP expression (Fig. 10B), whereas Tregs was negatively correlated with SDC1 expression (Fig. 10C). Three kinds of TICs; monocytes, dendritic cells (DCs), and eosinophils, were positively associated with CEBPB expression. B-cell memory and naïve B-cells positively correlated with AIP expression showed a negative association with VPS25 expression. Immune cell infiltration was not found to be related to the expression of the STXBP2 gene. Only one incidence of correlation (positive or negative) was observed between the expression of 5 m6A-related mRNAs and immune cell infiltration in this study.

Estimation of the m6A-related mRNA risk prognosis model for CHOL

The capability of independent prognostic signatures of the MRMRPM for CHOL was then estimated using univariate and multivariate Cox analyses. In this study, the higher risk score, HR, was based on univariate Cox analysis (HR = 2.72, 95% confidence interval (CI) = 1.77–4.18, p < 0.001) (Table 5) and multivariate Cox analysis (HR = 2.77, 95% CI = 1.81–4.25, p < 0.001) (Table 5, Fig. 11A), and the results obtained using the MRMRPM were not associated with tumor stage, gender, or age. The ROC curves for 1- and 2-year OS revealed that the MRMRPM could accurately predict the OS for the two CHOL groups with different risks (AUC = 0.913 and 0.903) (Fig. 11B). The AUC for risk was higher than other clinicopathological characteristics, and the risk was the only HR characteristic. The MRMRPM is therefore more reliable for use with CHOL (Fig. 11C). Compared to the other clinical signatures, the risk score obtained using the MRMRPM demonstrated significant predictive ability with nomogram analysis (Fig. 11D) and the calibrated curves suggest that the 1- and 2-year OS correction curves were consistent (Fig. 11E, F).

Table 5 Hazard ratio analyses of the clinical characteristics and risk score with the OS
Full size table
Fig. 11
figure 11

Estimate of MRMRPM and clinical characteristics in CHOL. A Hazard ratio of each independent risk signature is shown in the forest plot. B, C ROC curves of clinical characteristics and risk score to estimate prognostic sensitivity in CHOL. D Probability of 1- and 2-year OS by nomogram analysis. E, F) Calibration curves for 1-year (E) and 2-year (F) OS by MRMRPM. MRMRPM; five m6A-related mRNAs risk prognostic model; CHOL: Cholangiocarcinoma; ROC: receiver operating characteristics; OS: overall survival

Full size image

Discussion

Most patients with cholangiocarcinoma are diagnosed at an advanced stage. Over the past 25 years, the incidence and mortality of cholangiocarcinoma has increased, the prognosis of patients with cholangiocarcinoma has not improved significantly, and treatment is limited. The 5-year OS of cholangiocarcinoma is only 5–10%, and the median survival time of an individual with advanced cholangiocarcinoma is less than 12 months [2]. Therefore, it is vitally important that prognostic factors for CHOL are found. The Joint American Council on Cancer (AJCC) staging manual has become the benchmark for classifying cancers, predicting the prognosis, and determining the best treatment for cancer patients. However, the acceleration of cancer research has indicated that the AJCC-TNM stage system is not sufficient to estimate the prognosis or respond to the heterogeneity of these tumors. Even for patients with tumors at the same stage, the treatment response is highly heterogeneous, and other factors such as age, presentation status, and tumor site affect patient survival, which means that the information provided by the TNM staging system is limited in terms of prognosis [29]. Therefore, there is an urgent need to establish reliable prognostic biomarkers so prognoses can be improved. Recent studies investigating the use of screening for clinical prognostic mRNA indicators and other biomarkers have predicted the OS and the response of CHOL [30,31,32] and other cancers [33,34,35,36] to immunotherapy.

Higher TMB patients had better OS [35] when treated with ICIs. TMB is accompanied by an increase of tumor immune antigen, which is positively associated with an increasing CD8+ T cell count, and CD8+ T infiltration is the basis of immunotherapy. Previous studies have confirmed that an increase in the number of CD8+ T cells is not accompanied by an increase in the number of tumor antigens [37] for a variety of cancers, including breast and prostate. Analysis of data obtained from TCGA, the largest cancer gene database, indicates that not all cancers are accompanied by high CD8+ T-cell counts [38]. Thus, TMB-H tumor CD8+ T cells do not necessarily increase in response to the presence of tumors. McGrail et al. divided tumors into two categories based on the correlation with neoantigens; I: CD8 is positively correlated with neoantigens, and II: CD8 is not positively correlated with neoantigens. The results showed a significantly stronger immune checkpoint inhibitor response for high tumor mutation load (TMB-H) compared to low tumor mutation load (TMB-L) in category I cancers, whereas such effects were difficult to observe in category II cancers. In renal clear cell carcinoma and metastatic lung squamous cell carcinoma, TMB-L was found to lead to a prognosis better than TMB-H treated with immunotherapy [38]. After combining the data describing category II cancer, the objective response rate was only 15.3%, whereas the odds ratio of the objective response rate for TMB-L cancer was 0.46. TMB was significantly associated with the objective response rate in category I tumors, but not in category II tumors. The same result was observed in terms of OS, with a significant association detected for category I cancers and none for category II cancers. CHOL is a category II cancer [37], indicating that TMB is not prominently correlated with the objective response rate to immunotherapy in patients with this cancer. Marabelle et al. also found that TMB did not predict the efficacy of immunotherapy for CHOL patients [39]. Our study also revealed that TMB is not associated with OS (Fig. 9D), and no differences were observed in the OS of the high and low- TMB groups. These results indicate that TMB is not the only prognostic indicator of immunotherapy for CHOL patients and that other indicators should be considered in full. Therefore, we attempted to detect m6A-related mRNAs that contribute to survival and the classification of TNM stages in CHOL patients using the TCGA database.

The abnormal methylation of m6A RNA is closely related to the occurrence and development of tumors [40,41,42]. Although the molecular mechanism and role of m6A mediation in different tumors have not been fully elucidated, improvements in high-throughput sequencing and bioinformatics have led to the development of multiple methods for detecting and analyzing m6A methylation. The study of m6A and the associated participants that induce the process of reversible regulation (m6A-modifying enzymes and m6a-binding proteins) and the mechanisms involved in the occurrence and development of tumors have become a target for research [43, 44]. Data from the TCGA database showed that the mRNA expression of METTL3 is significantly elevated in lung adenocarcinomas (LUAD) [45]. The mRNA metabolism includes mRNA processing, nuclear export, translation, degradation, and other processes and the mechanisms by which m6A is involved in mRNA metabolism are still being revealed. The mRNA processing promotes the maturation of precursor mRNA via three steps and studies have found that m6A is more abundant in the precursor than mature mRNA [46]. Many m6A writers, erasers, and readers are primarily located in cellular nuclear spots [47,48,49] and the selective splicing factors that play an important role in the splicing of precursor mRNA are also enriched in this subnuclear structure. Ribosomal development occurs in the cellular nucleoli, with cellular nucleolar proteins playing an important role in the transcription and processing of ribosomal RNA (rRNA). Dysfunctional ribosomal function usually leads to changes in the cell cycle and the occurrence of tumors. METTL3 and FTO can influence the function of selective splicing factor 2 (SRSF2) by influencing the level of m6A, leading to regulation of the selective splicing process of mRNA [50]. In pancreatic adenocarcinoma research, a 16 m6A-related mRNA signature score system has been established for predicting the prognosis of pancreatic cancer patients based on the TCGA database [51]. In our study, we screened five m6A-related mRNAs from 1281 m6A-related mRNAs that utilized TCGA-CHOL cohorts for developing a MRMRPM. MRMRPM showed an accurate predictive ability with AUCs of 0.908 and 0.923 for the prediction of 1 and 2-year OS, respectively (Fig. 4E). Furthermore, we used two GEO datasets and visualization by PCA to validate the risk model and found that the MRMRPM model was accurate.

Among the five m6A-related mRNAs, CEBPB is involved in several of the most active biological processes associated with proliferation, such as biogenesis, the translation of ribosomes, and RNA processing [52]. The overexpression of CEBPB in C-REL-deficient myeloid-derived suppressor cell (MDSC) facilitates the expression of downstream CEBPB genes, enhancing mitochondrial respiration, weakening glycolysis, and inhibiting T cells and the immunotherapeutic response [53].

In this study, we found AIP and SDC1 were downstream target genes of CEBP (Supplementary Table S2). SDC1 is key in maintaining the survival and development of cancer cells [54], and SDC1 expression has been related to tumor differentiation, depth of invasion, lymph node metastasis, distant metastasis, TNM stage, and prognosis [55, 56]. The AIP protein could induced Epstein-Barr virus [57] and the hepatitis B virus [58] by involveing in cell transformation. We found high-risk group had higher expression of CEBPB, AIP and SDC1. the overexpression of CEBPB might facilitated the overexpression of SDC1 and AIP (Fig. 4C), resulted in worse OS in high-risk group in CHOL. Meanwhile, the expression of CEBPB and STXEP2 was related TICs function [58]. Fto knockdown suppress the expression of transcription factors c-Jun, JunB, and CEBPB, restores the function of CD8+ T cells by impairing the glycolytic activity of tumor cells, thereby inhibiting tumor growth [59]. NK cells have severely reduced/absent degranulation and cytotoxicity in STXBP2-deficient patients [60]. In our study, the high-risk group had higher expression of CEBPB and lower STXBP2 (Fig. 4C). The opposite pattern of expression might inhibited the activity of immune cells, thereby promoted tumor growth, resulted in worse OS in high-risk group in CHOL. More experiments are needed to confirm these conjecture.

Essential endocytoid sorting complex-II (ESCRT-II) is a preformed complex consisting of Vps36, Vps25 and Vps22 molecules and is recruited to endosomes [61, 62]. As an RNA-binding complex [63], ESCRT-II subcomplex can sort circRHOBTB3 into exosomes and secrete it out of the cell, resulting in tumor exosome escape mechanism [64]. In this study, the overexpression of VPS25 in the high-risk group (Fig. 4C) promoted the formation of ESCRT-II and delayed or prevented the necrotic apoptosis of tumor cells, thus contributed to worse OS in high-risk groups (Fig. 4D). More experiments are needed to confirm this conjecture.

In the present study, CHOL patients with high expression of the four m6A-related mRNAs (except STXBP2) were associated with worse OS (Supplementary Fig. S1A-E), and CHOL patients with a higher risk score were also associated with worse OS (Fig. 4D). These results indicate that both the presence of m6A-related mRNA and the MRMRPM can be used as prognostic indices for CHOL patients.

We also compared the infiltration levels of different immune cells to further confirm the relationship between the expression of m6A genes, m6A-related mRNAs, and tumor IME. Early studies confirmed that higher expression of the m6A genes (FTO, MELLT14, METTTL3, YTHDC1, YTHDC2, and ZC3H13) leads to higher immune infiltration, whereas patients with low immune infiltration showed higher expression of ALKBH5, HNRNPC, WTAP, YTHDF1, and YTHDF2 [59]. We found that high m6A-related mRNA expression promoted the infiltration of memory B-cells, naïve B-cells, plasma cells, Tregs, monocytes, resting dendritic cells, and eosinophils (Fig. 10, Supplementary Fig. S2). Reducing the level of infiltration of Tregs cells that inhibit the proliferation of CD4+ T cells [65, 66] can effectively improve the OS of cancer patients [67], and a higher level of infiltration of resting DCs is correlated with low OS [68]. In the present study, Treg infiltration levels of Tregs and resting dendritic cells were positively associated with high expression of the four m6A-related mRNAs (AIP, CEBPB, SDC1 and VPS25) (Fig. 10, Supplementary Fig. S2), which also showed high levels of expression in the high-risk CHOL group. These studies indicate that infiltration with Tregs and resting dendritic cells, which is promoted by the four m6A-related mRNAs, leads to the adverse outcomes observed in the high m6A-related mRNA risk group. This conclusion requires further verification in large-scale experiments.

Conclusions

In summary, our results provide evidence for the predictive prognosis of CHOL and demonstrated that TMB is not the only prognostic model that can be used to predict the prognosis for patients with CHOL. The m6A-related mRNA risk predictive model shows promise for screening patients with CHOL with better immunotherapy response.

Availability of data and materials

The datasets GSE89748 and GSE107943 for this study can be found in the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The datasets TCGA of CHOL for this study can be found in the TCGA project (https://portal.gdc.cancer.gov/).

Abbreviations

CHOL:

Cholangiocarcinoma

OS:

Overall survival

AUC:

Area under curve

ROC:

Receiver operating characteristics

MRMRPM:

Five m6A-related mRNAs risk prognostic model

CCLE:

Cancer Cell Line Encyclopedia

HR:

Hazard ratio

IME:

Immune microenvironment

AIP :

Aryl hydrocarbon receptor interacting protein

CEBPB :

CCAAT/enhancer binding protein beta

SDC1 :

Syndecan1

VPS25 :

Vacuolar protein sorting 25 homolog

STXBP2 :

Syntaxin Binding Protein 2

TCGA:

The Cancer Genome Atlas

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

PPI:

Protein-protein interaction

LASSO:

Least absolute shrinkage and selection operator

GSE:

Gene set enrichment

GEO:

Gene Expression Omnibus

PCA:

Principal component analysis

TMB:

Tumor mutational burden

DC:

Endritic cell

Tregs:

T cells regulatory

ESCRT:

Essential endocytoid sorting complex

TICs:

Tumor-infiltrating immune cells

NK:

Natural killer

METTL :

Methyltransferase-like proteins

SOCS2 :

Suppressor of cytokine signaling 2

References

  1. Turgeon MK, Maithel SK. Cholangiocarcinoma: a site-specific update on the current state of surgical management and multi-modality therapy. Chin Clin Oncol. 2020;9:4. https://doi.org/10.21037/cco.2019.08.09. 

    Article  PubMed  Google Scholar 

  2. Razumilava N, Gores GJ. Cholangiocarcinoma Lancet. 2014;383:2168–79. https://doi.org/10.1016/S0140-6736(13)61903-0.

    Article  PubMed  Google Scholar 

  3. Singhal D, van Gulik TM, Gouma DJ. Palliative management of hilar cholangiocarcinoma. Surg Oncol. 2005;14:59–74. https://doi.org/10.1016/j.suronc.2005.05.004.

    Article  CAS  PubMed  Google Scholar 

  4. Mihalache F, Tantau M, Diaconu B, Acalovschi M. Survival and quality of life of cholangiocarcinoma patients: a prospective study over a 4 year period. J Gastrointestin Liver Dis. 2010;19:285–90.

    PubMed  Google Scholar 

  5. Sun T, Wu R, Ming L. The role of m6A RNA methylation in cancer. Biomed Pharmacother. 2019;112:108613. https://doi.org/10.1016/j.biopha.2019.108613.

    Article  CAS  PubMed  Google Scholar 

  6. Chen XY, Zhang J, Zhu JS. The role of m6A RNA methylation in human cancer. Mol Cancer. 2019;18:103. https://doi.org/10.1186/s12943-019-1033-z.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Liu ZX, Li LM, Sun HL, Liu SM. Link between m6A modification and cancers. Front Bioeng Biotechnol. 2018;6:89. https://doi.org/10.3389/fbioe.2018.00089.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Chen M, Wei L, Law CT, Tsang FH, Shen J, Cheng CL, et al. RNA N6-methyladenosine methyltransferase-like 3 promotes liver cancer progression through YTHDF2-dependent posttranscriptional silencing of SOCS2. Hepatology. 2018;67:2254–70. https://doi.org/10.1002/hep.29683.

    Article  CAS  PubMed  Google Scholar 

  9. Cheng M, Sheng L, Gao Q, Xiong Q, Zhang H, Wu M, et al. The m6A methyltransferase METTL3 promotes bladder cancer progression via AFF4/NF-κB/MYC signaling network. Oncogene. 2019;38:3667–80. https://doi.org/10.1038/s41388-019-0683-z.

    Article  CAS  PubMed  Google Scholar 

  10. Jo HJ, Shim HE, Han ME, Kim HJ, Kim KS, Baek S, et al. WTAP regulates migration and invasion of cholangiocarcinoma cells. J Gastroenterol. 2013;48:1271–82. https://doi.org/10.1007/s00535-013-0748-7.

    Article  CAS  PubMed  Google Scholar 

  11. Rong ZX, Li Z, He JJ, Liu LY, Ren XX, Gao J, et al. Downregulation of fat mass and obesity associated (FTO) promotes the progression of intrahepatic cholangiocarcinoma. Front Oncol. 2019;9:369. https://doi.org/10.3389/fonc.2019.00369.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Cui Q, Shi H, Ye P, Li L, Qu Q, Sun G, et al. m6A RNA methylation regulates the self-renewal and tumorigenesis of glioblastoma stem cells. Cell Rep. 2017;18:2622–34. https://doi.org/10.1016/j.celrep.2017.02.059.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Zhou S, Bai ZL, Xia D, Zhao ZJ, Zhao R, Wang YY, et al. FTO regulates the chemo-radiotherapy resistance of cervical squamous cell carcinoma (CSCC) by targeting β-catenin through mRNA demethylation. Mol Carcinog. 2018;57:590–7. https://doi.org/10.1002/mc.22782.

    Article  CAS  PubMed  Google Scholar 

  14. Niu Y, Lin Z, Wan A, Chen H, Liang H, Sun L, et al. RNA N6-methyladenosine demethylase FTO promotes breast tumor progression through inhibiting BNIP3. Mol Cancer. 2019;18:46. https://doi.org/10.1186/s12943-019-1004-4.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Li Z, Weng H, Su R, Weng X, Zuo Z, Li C, et al. FTO plays an oncogenic role in acute myeloid leukemia as a N6-Methyladenosine RNA demethylase. Cancer Cell. 2017;31:127–41. https://doi.org/10.1016/j.ccell.2016.11.017.

    Article  CAS  PubMed  Google Scholar 

  16. Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30. https://doi.org/10.1093/nar/28.1.27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51. https://doi.org/10.1002/pro.3715.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2021;49(D1):D545–51. https://doi.org/10.1093/nar/gkaa970.

    Article  CAS  PubMed  Google Scholar 

  19. Xu F, Zhan X, Zheng X, Xu H, Li Y, Huang X, et al. A signature of immune-related gene pairs predicts oncologic outcomes and response to immunotherapy in lung adenocarcinoma. Genomics. 2020;112:4675–83. https://doi.org/10.1016/j.ygeno.2020.08.014.

    Article  CAS  PubMed  Google Scholar 

  20. Xu F, Lin H, He P, He L, Chen J, Lin L, et al. A TP53-associated gene signature for prediction of prognosis and therapeutic responses in lung squamous cell carcinoma. Oncoimmunology. 2020;9:1731943. https://doi.org/10.1080/2162402X.2020.1731943.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Lachmann A, Xu H, Krishnan J, Berger SI, Mazloom AR, Ma'ayan A. ChEA: transcription factor regulation inferred from integrating genome-wide ChIP-X experiments. Bioinformatics. 2010;26(19):2438–44. https://doi.org/10.1093/bioinformatics/btq466.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Li X, Li Y, Yu X, Jin F. Identification and validation of stemness-related lncRNA prognostic signature for breast cancer. J Transl Med. 2020;18:331. https://doi.org/10.1186/s12967-020-02497-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Mora A, Donaldson IM. iRefR: an R package to manipulate the iRefIndex consolidated protein interaction database. BMC Bioinformatics. 2011;12:455. https://doi.org/10.1186/1471-2105-12-455.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One. 2014;9:e107468. https://doi.org/10.1371/journal.pone.0107468.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Jiao Y, Pawlik TM, Anders RA, Selaru FM, Streppel MM, Lucas DJ, et al. Exome sequencing identifies frequent inactivating mutations in BAP1, ARID1A and PBRM1 in intrahepatic cholangiocarcinomas. Nat Genet. 2013;45(12):1470–3. https://doi.org/10.1038/ng.2813.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Masliah-Planchon J, Bièche I, Guinebretière JM, Bourdeaut F, Delattre O. SWI/SNF chromatin remodeling and human malignancies. Annu Rev Pathol. 2015;10:145–71 annurev-pathol-012414-040445.

    Article  CAS  PubMed  Google Scholar 

  27. Miao D, Margolis CA, Gao W, Voss MH, Li W, Martini DJ, et al. Genomic correlates of response to immune checkpoint therapies in clear cell renal cell carcinoma. Science. 2018;359(6377):801–6. https://doi.org/10.1126/science.aan5951.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12:453–7. https://doi.org/10.1038/nmeth.3337.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Amin MB, Greene FL, Edge SB, Compton CC, Gershenwald JE, Brookland RK, et al. The Eighth Edition AJCC Cancer Staging Manual: continuing to build a bridge from a population-based to a more “personalized” approach to cancer staging. CA Cancer J Clin. 8th ed. 8th ed. 8th ed. 2017;67:93–9. https://doi.org/10.3322/caac.21388.

    Article  PubMed  Google Scholar 

  30. He W, Zhang A, Qi L, Na C, Jiang R, Fan Z, et al. FOXO1, a potential therapeutic target, regulates autophagic flux, oxidative stress, mitochondrial dysfunction, and apoptosis in human cholangiocarcinoma QBC939 cells. Cell Physiol Biochem. 2018;45:1506–14. https://doi.org/10.1159/000487576.

    Article  CAS  PubMed  Google Scholar 

  31. Jones RP, Bird NT, Smith RA, Palmer DH, Fenwick SW, Poston GJ, et al. Prognostic molecular markers in resected extrahepatic biliary tract cancers; a systematic review and meta-analysis of immunohistochemically detected biomarkers. Biomark Med. 2015;9:763–75. https://doi.org/10.2217/BMM.15.48.

    Article  CAS  PubMed  Google Scholar 

  32. Chen TC, Jan YY, Yeh TS. K-ras mutation is strongly associated with perineural invasion and represents an independent prognostic factor of intrahepatic cholangiocarcinoma after hepatectomy. Ann Surg Oncol. 2012;19(Suppl 3):S675–81. https://doi.org/10.1245/s10434-012-2224-7.

    Article  PubMed  Google Scholar 

  33. Chen J, Chen Z, Huang Z, Yu H, Li Y, Huang W. Formiminotransferase cyclodeaminase suppresses hepatocellular carcinoma by modulating cell apoptosis, DNA damage, and phosphatidylinositol 3-kinases (PI3K)/Akt signaling pathway. Med Sci Monit. 2019;25:4474–84. https://doi.org/10.12659/MSM.916202.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Park SE, Park K, Lee E, Kim JY, Ahn JS, Im YH, et al. Clinical implication of tumor mutational burden in patients with HER2-positive refractory metastatic breast cancer. Oncoimmunology. 2018;7:e1466768. https://doi.org/10.1080/2162402X.2018.1466768.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Samstein RM, Lee CH, Shoushtari AN, Hellmann MD, Shen R, Janjigian YY, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet. 2019;51:202–6. https://doi.org/10.1038/s41588-018-0312-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Yarchoan M, Hopkins A, Jaffee EM. Tumor mutational burden and response rate to PD-1 inhibition. N Engl J Med. 2017;377:2500–1. https://doi.org/10.1056/NEJMc1713444.

    Article  PubMed  PubMed Central  Google Scholar 

  37. McGrail DJ, Pilié PG, Rashid NU, Voorwerk L, Slagter M, Kok M, et al. High tumor mutation burden fails to predict immune checkpoint blockade response across all cancer types. Ann Oncol. 2021;32:661–72. https://doi.org/10.1016/j.annonc.2021.02.006.

    Article  CAS  PubMed  Google Scholar 

  38. Vokes NI, Liu D, Ricciuti B, Jimenez-Aguilar E, Rizvi H, Dietlein F, et al. Harmonization of tumor mutational burden quantification and association with response to immune checkpoint blockade in non-small-cell lung cancer. JCO precis. Oncol. 2019;3. https://doi.org/10.1200/PO.19.00171.

  39. Marabelle A, Fakih M, Lopez J, Shah M, Shapira-Frommer R, Nakagawa K, et al. Association of tumour mutational burden with outcomes in patients with advanced solid tumours treated with pembrolizumab: prospective biomarker analysis of the multicohort, open-label, phase 2 KEYNOTE-158 study. Lancet Oncol. 2020;21:1353–65. https://doi.org/10.1016/S1470-2045(20)30445-9.

    Article  CAS  PubMed  Google Scholar 

  40. Fu Y, Dominissini D, Rechavi G, He C. Gene expression regulation mediated through reversible m6A RNA methylation. Nat Rev Genet. 2014;15:293–306. https://doi.org/10.1038/nrg3724.

    Article  CAS  PubMed  Google Scholar 

  41. Visvanathan A, Patil V, Arora A, Hegde AS, Arivazhagan A, Santosh V, et al. Essential role of METTL3-mediated m6A modification in glioma stem-like cells maintenance and radioresistance. Oncogene. 2018;37:522–33. https://doi.org/10.1038/onc.2017.351.

    Article  CAS  PubMed  Google Scholar 

  42. Liu N, Parisien M, Dai Q, Zheng G, He C, Pan T. Probing N6-methyladenosine RNA modification status at single nucleotide resolution in mRNA and long noncoding RNA. RNA. 2013;19:1848–56. https://doi.org/10.1261/rna.041178.113.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Panneerdoss S, Eedunuri VK, Yadav P, Timilsina S, Rajamanickam S, Viswanadhapalli S, et al. Cross-talk among writers, readers, and erasers of m6A regulates cancer growth and progression. Sci Adv. 2018;4:eaar8263. https://doi.org/10.1126/sciadv.aar8263.

  44. Liu X, Wang P, Teng X, Zhang Z, Song S. Comprehensive analysis of expression regulation for RNA m6A regulators with clinical significance in human cancers. Front Oncol. 2021;11:624395. https://doi.org/10.3389/fonc.2021.624395.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Lin S, Choe J, Du P, Triboulet R, Gregory RI, The M, et al. The m (6) a methyltransferase METTL3 promotes translation in human cancer cells. Mol Cell. 2016;62:335–45. https://doi.org/10.1016/j.molcel.2016.03.021.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Salditt-Georgieff M, Jelinek W, Darnell JE, Furuichi Y, Morgan M, Shatkin A. Methyl labeling of HeLa cell hnRNA: a comparison with mRNA. Cell. 1976;7:227–37. https://doi.org/10.1016/0092-8674(76)90022-2.

    Article  CAS  PubMed  Google Scholar 

  47. Ping XL, Sun BF, Wang L, Xiao W, Yang X, Wang WJ, et al. Mammalian WTAP is a regulatory subunit of the RNA N6-methyladenosine methyltransferase. Cell Res. 2014;24:177–89. https://doi.org/10.1038/cr.2014.3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Xiao W, Adhikari S, Dahal U, Chen YS, Hao YJ, Sun BF, et al. Nuclear m (6) a reader YTHDC1 regulates mRNA splicing. Mol Cell. 2016;61:507–19. https://doi.org/10.1016/j.molcel.2016.01.012.

    Article  CAS  PubMed  Google Scholar 

  49. Jia G, Fu Y, Zhao X, Dai Q, Zheng G, Yang Y, et al. N6-methyladenosine in nuclear RNA is a major substrate of the obesity-associated FTO. Nat Chem Biol. 2011;7:885–7. https://doi.org/10.1038/nchembio.687.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Zheng G, Dahl JA, Niu Y, Fedorcsak P, Huang CM, Li CJ, et al. ALKBH5 is a mammalian RNA demethylase that impacts RNA metabolism and mouse fertility. Mol Cell. 2013;49:18–29. https://doi.org/10.1016/j.molcel.2012.10.015.

    Article  CAS  PubMed  Google Scholar 

  51. Meng Z, Yuan Q, Zhao J, Wang B, Li S, Offringa R, et al. The m6A-related mRNA signature predicts the prognosis of pancreatic cancer patients. Mol Ther Oncolytics. 2020;17:460–70. https://doi.org/10.1016/j.omto.2020.04.011.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Wang H, Diaz AK, Shaw TI, Li Y, Niu M, Cho JH, et al. Deep multiomics profiling of brain tumors identifies signaling networks downstream of cancer driver genes. Nat Commun. 2019;10:3718. https://doi.org/10.1038/s41467-019-11661-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Li T, Li X, Zamani A, Wang W, Lee CN, Li M, et al. C-Rel is a myeloid checkpoint for Cancer immunotherapy. Nat Can. 2020;1:507–17. https://doi.org/10.1038/s43018-020-0061-3.

    Article  CAS  Google Scholar 

  54. Yao W, Rose JL, Wang W, Seth S, Jiang H, Taguchi A, et al. Syndecan 1 is a critical mediator of macropinocytosis in pancreatic cancer. Nature. 2019;568:410–4. https://doi.org/10.1038/s41586-019-1062-1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Andersen NF, Kristensen IB, Preiss BS, Christensen JH, Abildgaard N. Upregulation of Syndecan-1 in the bone marrow microenvironment in multiple myeloma is associated with angiogenesis. Eur J Haematol. 2015;95:211–7. https://doi.org/10.1111/ejh.12473.

    Article  CAS  PubMed  Google Scholar 

  56. Purushothaman A, Uyama T, Kobayashi F, Yamada S, Sugahara K, Rapraeger AC, et al. Heparanase-enhanced shedding of syndecan-1 by myeloma cells promotes endothelial invasion and angiogenesis. Blood. 2010;115:2449–57. https://doi.org/10.1182/blood-2009-07-234757.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Jung YW, Shim JI, Shim SH, Shin YJ, Shim SH, Chang SW, et al. Global gene expression analysis of cell-free RNA in amniotic fluid from women destined to develop preeclampsia. Medicine (Baltimore). 2019;98(3):e13971. https://doi.org/10.1097/MD.0000000000013971.

    Article  CAS  Google Scholar 

  58. Osuna CE, Gonzalez AM, Chang HH, Hung AS, Ehlinger E, Anasti K, et al. TCR affinity associated with functional differences between dominant and subdominant SIV epitope-specific CD8+ T cells in Mamu-a*01+ rhesus monkeys. PLoS Pathog. 2014;10(4):e1004069. https://doi.org/10.1371/journal.ppat.1004069.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Liu Y, Liang G, Xu H, Dong W, Dong Z, Qiu Z, et al. Tumors exploit FTO-mediated regulation of glycolytic metabolism to evade immune surveillance. Cell Metab. 2021;33(6):1221–1233.e11. https://doi.org/10.1016/j.cmet.2021.04.001.

    Article  CAS  PubMed  Google Scholar 

  60. Lopez JA, Noori T, Minson A, Li Jovanoska L, Thia K, Hildebrand MS, et al. Bi-allelic mutations in STXBP2 reveal a complementary role for STXBP1 in cytotoxic lymphocyte killing. Front Immunol. 2018;9:529. https://doi.org/10.3389/fimmu.2018.00529.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Babst M, Katzmann DJ, Snyder WB, Wendland B, Emr SD. Endosome-associated complex, ESCRT-II, recruits transport machinery for protein sorting at the multivesicular body. Dev Cell. 2002;3(2):283–9. https://doi.org/10.1016/s1534-5807(02)00219-8.

    Article  CAS  PubMed  Google Scholar 

  62. Teo H, Perisic O, González B, Williams RL. ESCRT-II, an endosome-associated complex required for protein sorting: crystal structure and interactions with ESCRT-III and membranes. Dev Cell. 2004;7(4):559–69. https://doi.org/10.1016/j.devcel.2004.09.003.

    Article  CAS  PubMed  Google Scholar 

  63. van Niel G, D'Angelo G, Raposo G. Shedding light on the cell biology of extracellular vesicles. Nat Rev Mol Cell Biol. 2018;19(4):213–28. https://doi.org/10.1038/nrm.2017.125.

    Article  CAS  PubMed  Google Scholar 

  64. Chen C, Yu H, Han F, Lai X, Ye K, Lei S, et al. Tumor-suppressive circRHOBTB3 is excreted out of cells via exosome to sustain colorectal cancer cell fitness. Mol Cancer. 2022;21(1):46. https://doi.org/10.1186/s12943-022-01511-1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Shi C, Zhang Y, Yang H, Dong T, Chen Y, Xu Y, et al. Ultrasound-targeted microbubble destruction-mediated Foxp3 knockdown may suppress the tumor growth of HCC mice by relieving immunosuppressive Tregs function. Exp Ther Med. 2018;15:31–8. https://doi.org/10.3892/etm.2017.5421.

    Article  CAS  PubMed  Google Scholar 

  66. Yu A, Snowhite I, Vendrame F, Rosenzwajg M, Klatzmann D, Pugliese A, et al. Selective IL-2 responsiveness of regulatory T cells through multiple intrinsic mechanisms supports the use of low-dose IL-2 therapy in type 1 diabetes. Diabetes. 2015;64:2172–83. https://doi.org/10.2337/db14-1322.

    Article  CAS  PubMed  Google Scholar 

  67. Sato E, Olson SH, Ahn J, Bundy B, Nishikawa H, Qian F, et al. Intraepithelial CD8+ tumor-infiltrating lymphocytes and a high CD8+/regulatory T cell ratio are associated with favorable prognosis in ovarian cancer. Proc Natl Acad Sci U S A. 2005;102:18538–43. https://doi.org/10.1073/pnas.0509182102.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Liu S, Tang Q, Huang J, Zhan M, Zhao W, Yang X, et al. Prognostic analysis of tumor mutation burden and immune infiltration in hepatocellular carcinoma based on TCGA data. Aging (Albany NY). 2021;13:11257–80. https://doi.org/10.18632/aging.202811.

    Article  CAS  Google Scholar 

Download references

Acknowledgments

We would like to thank Editage (www.editage.cn) for English language editing.

Funding

This study was supported by the Natural Science Foundation of Shandong Province, Project No. ZR2019MH089; Jinan Science and Technology Plan (Clinical Medical Science and Technology Innovation Plan), Project no. 202019057.

Author information

Authors and Affiliations

  1. Department of Hepatobiliary Surgery, Shandong Provincial Hospital, Cheeloo College of Medicine, Shandong University, Jinan, 250021, Shandong, China

    Huaqiang Zhu, Jianlu Wang, Chaoqun Ma, Hengjun Gao, Faji Yang, Qingqiang Ni, Hongguang Li, Xu Zhou & Jun Lu

  2. Department of Hepatobiliary Surgery, Shandong Provincial Hospital Affiliated to Shandong First Medical University, Jinan, 250021, Shandong, China

    Huaqiang Zhu, Jianlu Wang, Shuchao Zhao, Chaoqun Ma, Hengjun Gao, Faji Yang, Qingqiang Ni, Hongguang Li, Xu Zhou & Jun Lu

  3. Jinan Health Publicity and Education Center, Jinan, 250021, Shandong, China

    Haini Zhao

  4. ChosenMed Technology (Beijing) Co., Ltd., Beijing, 100176, China

    Dongliang Wang

  5. Department of Gastroenterology, Shandong Provincial Hospital, Cheeloo College of Medicine, Shandong University, Jinan, 250021, Shandong, China

    Chunqing Zhang

Authors
  1. Huaqiang Zhu
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Haini Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Jianlu Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Shuchao Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Chaoqun Ma
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Dongliang Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Hengjun Gao
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Faji Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Qingqiang Ni
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Hongguang Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Xu Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Chunqing Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  13. Jun Lu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

H.Z. conceived, designed, and wrote the manuscript. H.Z., J.W., and S.Z. collected the data, performed bioinformatics analysis, and generated the figures and Tables. H.Z., J.W., S.Z., C.M., and D.W., performed the literature search. H.G., F.Y., Q.N., H.L., and X.Z. revised the images and tables and checked the manuscript. C.Z. and J.L. supervised the study and preparation of the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Chunqing Zhang or Jun Lu.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no conflicts of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Supplementary Table S1.

m6A-related mRNA set (n = 1281).

Additional file 2: Supplementary Fig. S1.

Survival analysis of CHOL patients with m6A-related mRNA expression. A-E) Kaplan–Meier survival of AIP (A), CEBPB (B), SDC1 (C), VPS23 (D) and STXBP2 (E) in CHOL patient. CHOL: Cholangiocarcinoma; AIP: Aryl hydrocarbon receptor interacting protein; CEBPB: CCAAT/enhancer binding protein beta; SDC1:syndecan1; VPS25: vacuolar protein sorting 25 homolog; STXBP2: syntaxin binding protein 2.

Additional file 3: Supplementary Table S2.

Candidata target genes of CCAAT/enhancer binding protein beta (n = 1742)

Additional file 4: Supplementary Table S3.

Comparison of 22 tumor-infiltrating immune cell types between the low- and high-risk groups.

Additional file 5: Supplementary Table S4.

Comparison of 22 tumor-infiltrating immune cell types between the low- and high-risk groups.

Additional file 6: Supplementary Fig. S2.

Correlation between the proportion of tumor-infiltrating cells in the two risk groups with expression of the five m6A-related mRNAs. A) Scatterplot of the relationship between the abundance of the 22 TICs and STXBP2 expression. B) Scatterplot of the relationship between the abundance of the 22 TICs and CEBPB expression. C) Scatterplot of the relationship between the abundance of 22 TICs and VPS25 expression. CEBPB: CCAAT/enhancer binding protein beta; VPS25: vacuolar protein sorting 25 homolog; STXBP2: syntaxin binding protein 2.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhu, H., Zhao, H., Wang, J. et al. Potential prognosis index for m6A-related mRNA in cholangiocarcinoma. BMC Cancer 22, 620 (2022). https://doi.org/10.1186/s12885-022-09665-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12885-022-09665-3

Keywords

  • Cholangiocarcinoma
  • N6-methyladenosine
  • m6A-related mRNA
  • Risk model
  • Prognosis

BMC Cancer

ISSN: 1471-2407

Contact us