Quantification of m6A RNA methylation modulators pattern was a potential biomarker for prognosis and associated with tumor immune microenvironment of pancreatic adenocarcinoma

Abstract

Background: m6A is the most prevalent and abundant form of mRNA modifications and is closely related to tumor proliferation, differentiation, and tumorigenesis. In this study, we try to conduct an effective prediction model to investigated the function of m6A RNA methylation modulators in pancreatic adenocarcinoma and estimated the potential association between m6A RNA methylation modulators and tumor microenvironment infiltration for optimization of treatment.

Methods: Expression of 28 m6A RNA methylation modulators and clinical data of patients with pancreatic adenocarcinoma and normal samples were obtained from TCGA and GTEx database. Differences in the expression of 28 m6A RNA methylation modulators between tumour (n = 40) and healthy (n = 167) samples were compared by Wilcoxon test. LASSO Cox regression was used to select m6A RNA methylation modulators to analyze the relationship between expression and clinical characteristics by univariate and multivariate regression. A risk score prognosis model was conducted based on the expression of select m6A RNA methylation modulators. Bioinformatics analysis was used to explore the association between the m6Ascore and the composition of infiltrating immune cells between high and low m6Ascore group by CIBERSORT algorithm. Evaluation of m6Ascore for immunotherapy was analyzed via the IPS and three immunotherapy cohort. Besides, the biological signaling pathways of the m6A RNA methylation modulators were examined by gene set enrichment analysis (GSEA).

Results: Expression of 28 m6A RNA methylation modulators were upregulated in patients with PAAD except for MTEEL3. An m6Ascore prognosis model was established, including KIAA1429, IGF2BP2, IGF2BP3, METTL3, EIF3H and LRPPRC was used to predict the prognosis of patients with PAAD, the high risk score was an independent prognostic indicator for pancreatic adenocarcinoma, and a high risk score presented a lower overall survival. In addition, m6Ascore was related with the immune cell infiltration of PAAD. Patients with a high m6Ascore had lower infiltration of Tregs and CD8⁺T cells but a higher resting CD4⁺ T infiltration. Patients with a low m6Ascore displayed a low abundance of PD-1, CTLA-4 and TIGIT, however, the IPS showed no difference between the two groups. The m6Ascore applied in three immunotherapy cohort (GSE78220, TCGA-SKCM, and IMvigor210) did not exhibit a good prediction for estimating the patients' response to immunotherapy, so it may need more researches to figure out whether the m6A modulator prognosis model would benefit the prediction of pancreatic patients' response to immunotherapy.

Conclusion: Modulators involved in m6A RNA methylation were associated with the development of pancreatic cancer. An m6Ascore based on the expression of IGF2BP2, IGF2BP3, KIAA1429, METTL3, EIF3H and LRPPRC is proposed as an indicator of TME status and is instrumental in predicting the prognosis of pancreatic cancer patients.

Keywords: Immune cell infiltration; Immunotherapy; Pancreatic adenocarcinoma; Prognostic; m6A RNA methylation modulators; m6Ascore.

Fig. 1

Expression of m6A RNA methylation modulators in 140 patients with PAAD and 167 normal samples based on TCGA and GTEx database. The expression levels of 28 m6A RNA methylation modulators in patients with PAAD presented in A a heatmap (green means high expression while blue means low expression and blue means normal samples while pink means tumour samples) and B a boxplot (blue means normal samples while red means tumour samples). The association between the expression of the modulators and survival time with statistically significant C The CNV variation frequency of 28 m6A regulators in 140 patients with PAAD. D The location of CNV alteration of m6A regulators on 23 chromosomes. E Spearman’s correlation analysis of the 28 m6A RNA methylation modulators in PAAD patients. F Interaction associations among 28 m6A RNA methylation modulators are visualized in the network and the relationship between m6A regulators and prognosis. *p < 0.05, **p < 0.01, ***p < 0.001

Fig. 2

Prognostic value of risk score from m6A RNA methylation modulators expression in PAAD patients. A Cox univariate analysis of m6A RNA methylation modulators. B and C Selection of genes by LASSO Cox regression. D Nomogram predicting 1-, 2-, and 3-year survival of patients with PAAD based on the expression of selected m6A RNA methylation modulators. The top row shows the point value for each variable, rows 2–6 indicate the expression of selected m6A RNA methylation modulators, the sum of these values is located on the Total Points axis, and the line drawn downward to the survival axes is used to determine the likelihood of 1-, 2- or 3-year survival. Calibration of nomograms predicting. Calibration of nomograms predicting E 1-year, F 2-year, and G 3-year survival

Fig. 3

Relationship between m6A risk score and OS of PAAD patients. A landscape of risk grouped by risk score. B Survival analyses for low and high m6Ascore patient groups using Kaplan-Meier curves. C Time-dependent ROC analysis of risk score in predicting prognoses. D The difference of clinical features between high and low risk group

Fig. 4

The risk score is an independent prognostic indicator of PAAD patients. A Univariate and B multivariate Cox regression analysis of the association between clinical and pathological factors, including the risk score and OS of PAAD patients. C Nomogram based on age, gender, tumour grade and stage, and risk score. Calibration of nomograms predicting D 1-year, E 2-year, and F 3-year survival

Fig. 5

Identifying differentially expressed immune-related genes and enrichment analysis. A Differentially expressed genes between high and low m6Ascore group and visualized in volcano plot. B Intersection between differentially expressed genes from TCGA database and immune-related genes from Immport database showed in venn diagram. C Heatmap of the top 50 intersection immune-related genes with higher logFC. D Gene ontology enrichment analysis of intersection immune-related genes, BP biology process; CC cellular component; MF molecular function. E KEGG pathway enrichment analysis of intersection immune-related genes

Fig. 6

The associated KEGG pathways of m6A RNA methylation genes by GSEA. A EIF3H. B IGF2BP2. C IGF2BP3. D LRPPRC. E KIAA1429

Fig. 7

The association between tumor microenvironment and m6Ascore. Discrepancy of stromalscore (A), immunescore (B) and estimatescore (C) in two groups. Spearman’s correlation analysis of m6Ascore with tumor microenvironment. D stromalscore. E immunescore. F estimatedscore

Fig. 8

The different tumor-infiltrating immune cells between high and low m6Ascore patients with PAAD. A Profiles of 22 type of tumor-infiltrating immune cells in two group. The relationship between OS and A resting memory CD4+ T cells, B CD8+ T cells, C Tregs cells

Fig. 9

IPS and immunotherapy gene expression analysis. The levels of immune checkpoints molecules including PD- L1 (A), PD-1 (B), CTLA-4 (C), TIM-3 (D), LAG-3 (E), and TIGIT (F) in low-risk and high-risk groups. The association between IPS and the m6Ascore in PAAD patients based on TCIA database, H CTLA4⁻ PD1⁻ I CTLA4⁻ PD1⁺ J CTLA4⁺ PD1⁻ CTLA4⁺PD1⁺

Fig. 10

IHC analysis of six m6A modulators in pancreatic cancer tissues and adjacent normal tissues