Identification of m6A-Related lncRNAs in PDAC Patients
The work flow was shown in Fig. 1. Cumulatively, 140 PDAC patients from TCGA cohort and 63 PDAC patients from ICGC cohort were included in our study; for whom, the baseline clinical features were presented in Additional file 1: Table S1. Firstly, we identified 14830 lncRNAs in the TCGA cohort and 12559 lncRNAs in the ICGC cohort. Next, we extracted the expression data of lncRNAs and 24 m6A-related genes from the TCGA and the ICGC cohorts. Pearson correlation analyses were performed to identify m6A-related lncRNAs in two cohorts (|Pearson R| > 0.4 and p < 0.001). We obtained the 2672 correlations in the TCGA cohort and 21008 correlations in ICGC cohorts, then the corresponding relationship between m6A-related genes and lncRNAs in two cohorts were intersected. Finally, 585 shared correlations were obtained in both two cohorts, then 262 m6A-related lncRNAs were extracted from the shared correlations.
Identification of prognostic m6A-Related lncRNAs
Univariate Cox regression was applied to identify prognostic m6A-related lncRNAs from the 262 m6A-related lncRNAs in the TCGA traing cohort (p < 0.05). The result showed that 28 m6A-related lncRNAs were significantly associated with the OS (Fig. 2A), including 2 risky lncRNAs and 26 protective lncRNAs, and the correlations between the 28 lncRNAs and the m6A-related genes in the TCGA cohort are shown in Fig. 2B.
Establishment of the m6A-LPR in the TCGA Training Cohort
To construct the m6A-LPR for predicting the OS of PDAC patients, LASSO Cox analysis was performed to the 28 prognostic m6A-related lncRNAs in the TCGA cohort and 12 m6A-related lncRNAs were screened (Fig. 3A). Next, multivariate Cox proportional regression analysis was performed for analyzing these 12 m6A-related lncRNAs in TCGA cohort to construct the m6A-LPR. The m6A-LPR comprising four m6A-related lncRNAs, including MIR4435-1HG, RP5-1112D6.4, RP11-582J16.5 and RP11-999E24.3, was developed through the summary of the expression values of these four m6A-related lncRNAs multiplied by corresponding coefficients derived from the above multivariable Cox regression analysis (Fig. 3B, C). The downregulated RP5-1112D6.4, RP11-582J16.5, and RP11-999E24.3 with HR < 1 were considered to be tumor suppressors, while the upregulated MIR4435-1HG with HR > 1 was considered to be oncogenes. The Kaplan–Meier survival curves showed that higher expression of RP5-1112D6.4, RP11-582J16.5 and RP11-999E24.3 and lower expression of MIR4435-1HG were correlated with impoved OS in the TCGA cohort (Fig. 3D-G). The expression of four m6A-related lncRNAs were also associated with the clinicopathological and immune signatures of PDAC, such as WHO grade, TP53 mutation status, KRAS mutation status, ESTIMATE score, immune score, stromal score, and tumor mutation burden (TMB) (Fig. 3H). The genomic information of these four m6A-related lncRNAs and the corresponding correlation with m6A regulators were presented in Table 1. The risk score for each patient was calculated by following formula based on m6A-LPR: [(0.205)*expression value of MIR4435-1HG] - [(0.525)*expression value of RP11-582J16.5] - [(0.658)*expression value of RP11-999E24.3] - [(0.275)*expression value of RP5-1112D6.4]. Based on the median value of risk scores, patients were categorized into low-risk and high-risk groups. Kaplan–Meier survival curves showed that PDAC patients with lower risk scores had better OS (p = 0.0011) (Fig. 4A). Survival status and risk score distributions were illustrated in Fig. 4B. The ROC curves demonstrated that m6A-LPR had a good performance for predicting OS in the TCGA cohort (1-year AUC = 0.760, 2-year AUC = 0.722; Fig. 4C).
Validation of the m6A-LPR in the ICGC Cohort
In order to further validate the robustness of m6A-LPR, the risk scores were calculated in the ICGC cohort using the above equation. Based on the median risk score, PDAC patients in the ICGC cohort were also divided into low- and high-risk subgroups. Consistent with the results in the TCGA training cohort: PDAC patients with lower risk scores had better OS in the ICGC validation cohort (p = 0.020) (Fig. 4D). Survival status and risk score distributions were shown in Fig. 4E. The ROC curves also demonstrated that m6A-LPR had a prognostic value for PDAC patients in the ICGC cohort (1-year AUC = 0.657, 2-year AUC = 0.729; Fig. 4F). These results showed that the m6A-LPR based prognostic signature had a robust and stable ability in prognosis prediction for PDAC.
Principal Component Analysis
Principal component analysis (PCA) was applied to evaluate the discrepancies between the low- and high-risk subgroups based on the expression of the four m6A-related lncRNAs in m6A-LPR (Fig. 5A, B). The results showed that the samples screened by the four m6A-related lncRNAs could clearly divide the whole patients into a low-risk and high-risk group in both the TCGA and ICGC cohorts.
Stratification Analysis of the m6A-LPR in clinicopathological features
PDAC patients with WHO grade III-IV, mutant TP53 and mutant KRAS (Fig. 6A–C) had higher risk scores, whereas the risk scores were not correlated with stage, T category and N category (Fig. 6D-F). To evaluate whether m6A-LPR was an independent prognostic factor for PDAC patients, univariate and multivariate Cox analyses were performed. In the TCGA cohort, univariate Cox analysis showed that m6A-LPR was significantly associated with OS (HR: 2.72, 95% CI: 1.94–3.81, p < 0.001) and multivariate Cox analysis indicated that m6A-LPR was an independent predictor of OS (HR: 2.77, 95% CI: 1.93–3.96, p < 0.001; Fig. 6G, H). In the ICGC validation cohort, univariate and multivariate Cox analyses also indicated that m6A-LPR was an independent predictor of OS for PDAC patients (univariate: HR: 1.83, 95% CI: 1.12–2.66, p = 0.012; multivariate: HR: 1.75, 95% CI: 1.12–2.50, p = 0.020; Fig. 6I, J). These results demonstrated that m6A-LPR might be helpful for clinical prognosis evaluation as an independent prognostic indicator.
Stratification Analysis of the m6A-LPR in immune features
The relationship between m6A-LPR and tumor immunity was further evaluated. Tumor purity, TMB, Immune Checkpoint Molecules and the infiltration level of immune cell were estimated. PDAC patients in the high-risk group had remarkably lower stromal, immune, and ESTIMATE scores, indicating a lower level of stroma, immune cell infiltration, and tumor purity (Fig. 7A). Furthermore, PDAC patients in the high-risk group had significantly higher levels of TMB (Fig. 7B) and lower levels of CTLA4 expression (Fig. 7C). To further investigate the underlying molecular mechanisms of the m6A-LPR and its relevance to tumor immunity, the relative abundance of 22 tumor-infiltrating immune cells was assessed for each patient using CIBERSORT. The high-risk group was associated with significantly higher levels of plasma B cells and resting NK cells infiltration, and lower levels of infiltrating resting memory CD4 T cells, monocytes and resting mast cells (Fig. 7D).
Functional and pathway enrichment analysis
To explore the potential biological processes and pathways of the molecular discrepancy between the low-risk and high-risk groups, 927 differential expression genes (DEGs) were identified between the low-risk and high-risk groups in the TCGA cohort (|log2 (fold change) | > 1 and p < 0.05). Functional and pathway enrichment analysis indicated these DEGs were mainly enriched in these aspects: digestion, neuronal system and peptide hormone metabolism (Reactome Gene Sets); NABA matrisome-associated (Canonical Pathways); pancreatic secretion, neuroactive ligand-receptor interaction and cytokine-cytokine receptor interaction (KEGG Pathways); regulation of ion transport, chemical synaptic transmission, regulation of system process, signal release, neuropeptide signaling pathway and second-messenger-mediate signaling (GO Biological Processes) (Fig. 8A–C). These results could give us some insights into the potential molecular mechanisms of the m6A-LPR.
Explore Potential Drugs that have a Therapeutic Effect on PDAC
From the Drug-LncRNA Module of the LncMAP database, we obtained 304304 drug-lncRNA interaction pairs. A total of 28 prognostic m6A-related lncRNAs were then imported into the database to predict the potential drugs of the genes, and 75 drug-lncRNA interactions were extracted when FDR < 0.05. The network including 18 prognostic m6A-related lncRNAs and 18 drugs was identified (Fig. 9). The five most interaction with prognostic m6A-related lncRNAs drugs were Panobinostat, L-685458, Palbociclib, Crizotinib and TAE684.