Expression of m6A RNA methylation regulators and m6A-regulated lncRNAs in PDAC
The m6A RNA methylation regulators play an indispensable role in tumor occurrence, proliferation, invasion and metastasis. We probe and systematically download the 23 m6A regulator genes expression and lncRNA expression of 178 tumor samples and 4 adjacent paracancerous samples from TCGA database and 167 normal samples from GTEx database. All the expressions of 23 m6A RNA methylation regulator genes differed between tumor and normal samples (P<0.01).
From the heatmap and boxplot (Figure S1A and Figure S1B) of 23 m6A RNA methylation regulator genes expression, it can be observed that “writers” (ZC3H13, RBM15, RBM15B), “readers” (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGFBP1, IGFBP3), and “erasers” (ALKBH5) were up-regulated in tumor tissues compared to normal tissues. The remaining m6A regulator genes are highly expressed in normal tissues compared with tumor tissues. To further identify the relationship between m6A methylation regulators and lncRNAs, Pearson correlation coefficient and RNA-sequence data was used to select the m6Aregulated lncRNAs. Finally, a total of 1411 lncRNAs with absolute correlation coefficient >0.4 and P<0.05 was identified as lncRNAs coexpressed with m6A. (Table S1). Through univariate Cox analyses (Figure 2A), 28 m6A-regulated lncRNAswere obtained.. Among that, 18 m6A-regulated lncRNAs were protective genes, 10 m6A-regulated lncRNAs were risk genes. . After that, we compared the expression differences of 28 m6A-regulated lncRNAs between tumor and control groups in TCGA database combined with GTEx database, and the visualized heatmaps and boxplots are shown in Figure 2B. The LINC01559, LINC00941, LINC01705, CASC8, AP000695.2 and SH3PXD2A−AS1 tend to beupregulated sharply in tumor tissues.
Correlation between consensus clustering analysis of m6A-regulated lncRNAs and survival of PDAC patients
Concordance clustering analysis was performed on PDAC samples, which were identified in the cumulative distribution function (CDF), k takes values ranging from 2-9 and k denotes cluster counting (Figure S2). According to similarity of m6A-regulated lncRNAs expression and ambiguous proportion for clustering measurements, the optimal clustering parameter was k = 2 (Figure 3A) We merged the survival information and target m6A-regulated lncRNAs expression data of PDAC patients. A total of 177 PDAC patients were grouped into Cluster 1 (n=87) and Cluster 2 (n=90) based on the m6A-regulated lncRNAs expression. Comparing the clinical variables and demographic variables between Cluster 1 and Cluster 2 according to 28 m6A-regulated lncRNAs, it was found that there was no difference, but there was a significant difference in the overall survival between the two groups. The overall survival of Cluster 1 was superior to that of Cluster 2 (P<0.001) (Figure 3B).
Correlation between PD-L1 and other immune checkpoints gene expression and m6A-related lncRNAs in Cluster 1 and Cluster 2
To explore the relationship between m6A-regulated lncRNAs and immune checkpoints (ICPs) gene expression, we compared the distribution of ICPs gene between tumor group and paracancerous cohort and between Cluster 1 and Cluster 2(Figure 7A-H). There was no difference in PD-L1 gene expression between Cluster 1 and Cluster 2 as well as tumor tissues and normal tissues. (Figure 7A, B). However, the expression of the CD47 (P<0.05) (Figure 7C), CD276 (P<0.001) (Figure 7E), IDO-1 (P<0.001) (Figure 7G) genes that potential be ICPs was statistically different between the tumor and normal groups, showing significantly upregulated expression in the tumor tissue compared to the normal tissue. Meanwhile, the CD276 (B7-H3) gene also showed expression discrepancy in Cluster 1 and Cluster 2 (P<0.001) (Figure 7F). Co-expression analysis of the genes showed a positive correlation between PD-L1 expression and m6A-regulated lncRNA U62317.1 and LINC01705 (Figure 7I).
Stratification analysis of the m6A‑related lncRNAs in immune infiltration and tumor microenvironment
The tumor microenvironment (TME) is complex and constantly evolving. In addition to stromal cells, fibroblasts, and endothelial cells, the TME contains immune cells [24]. Stromal cells in the TME may secrete growth factors, which may not only promote tumor cell growth, but may also chemoattract other cells to migrate into the TME. Immune cells like natural killer and T cells, which can eliminate tumor cells, portend a promising outcome. Researchers have endeavored to understand the characteristics of the TME and the status of immune cells and their interactions with tumor cells. Certain antibody therapies have demonstrated having the ability to activate a patient's own immune system to kill tumor cells [25]. We further evaluate the association between m6A-regulated lncRNAs and tumor immune infiltration. Stromal and infiltration level of 22 immune cells was evaluated by calculating the stromal, immune and ESTIMATE score using CIBERSORT. In the subtyping of PDAC samples based on 28-m6A regulated lncRNAs, it can be seen that B cells naive, B cells memory, Plasma cells, CD8 T cells, T cells CD4 memory resting, Monocytes, Macrophages M0, and Dendritic cells activated in Cluster 1 and Cluster 2 have differential distributions (Figure 8A). In Cluster 1, B cells naive, plasma cells, CD8 T cells and Monocytes were infiltrated at higher levels in the TME, whereas in Cluster 2, expression levels of B cells memory, T cells CD4 memory resting and Macrophages M0 were upregulated (Figure 8A). The Cluster 1 had higher stromal score, immune score and estimate score, implicated lower tumor purity (Figure 8B, C, D). All of these scores have statistically significant (P<0.001).
Gene Set Enrichment Analyses (GSEA)
The Gene Set Enrichment Analyses (GSEA) was performed to explore the signature of m6A regulated lncRNAs in Cluster 1 and Cluster 2. Filtering criteria were FDR < 0.05 to find notable enriched pathways. As shown in Figure 9, the obvious enriched signature pathway in Cluster 2 was “cell cycle” in Gene Set Enrichment Anaylses (GSEA). Unfortunately, no significant signature enrichment was found in Cluster 1.
Construction of m6A-regulated lncRNAs prognostic nomogram
The 177 PDAC samples with survival time and status in TCGA were divided into a train set (N=108) and a test set (N=69) according to 6:4 proportion. Based on the 28 m6A-regulated lncRNAs as well as the LASSO regression analysis (Figure S3), we finally obtained 11 m6A-regulatedlated lncRNAs correlated with prognosis in train set. The risk score for each patient in train and sets was calculated as follows: risk score = U62317.1 * 0.340993840 + LINC00941 * 0.120489300 + AL139287.1 * (-0.173219714) + AC012615.1 * (-0.075084474) + CASC8 * 0.229254017 + Z97832.2 * (-0.279993778) + AC099778.1 * (-0.037908532) + AC245041.2 * 0.058969635 + AC091057.1 * 0.580745126 + TSPOAP1-AS1 * (-1.587783887) + AC005332.3 * (-0.324442618). After that, the respective patients were divided into high and low risk groups based on the median risk score of the training set in both the training and testing sets. Multivariate Cox analysis in the train set showed that age, tumor grade, and the risk score were independent risk factors associated with prognosis (Figure 4A). Multivariate Cox analysis in the test set also indicated that the risk score was an independent risk factor associated with prognosis (Figure 4B). The distributions of the risk score and survival status of training set and testing set are presented in Figure 5(training set: A, C, E and G; testing set: B, D, F and H). The Kaplan-Meier curve of OS and log-rank t test was obtained in training set and testing set, patient’s survival probabilities in the two datasets were statistically differed (P<0.001 in training set and P=0.008 in testing set) (Figure 5E, F). The ROC predicting 5-year survival probabilities in training set and testing set were presented in Figure 5G and Figure 5H. The AUC for 1-year, 3-year and 5-year OS in training set was 0.861, 0.829 and 0.885 and in testing set was 0.708, 0.785 and 0.818. Therefore, we constructed a prognostic risk nomogram based on the risk score (Figure 6).
The clinicopathological factors, clusters and immune-scores associated with risk scores in PDAC
Based on the median score of 11 prognostic m6A-regulated lncRNAs, 177 PDAC patients were divided into low-risk and high-risk groups. We analyzed the clinicopathological variables, cluster analysis, gene mutations and immune scores in the low-risk and high-risk groups in both the training and testing sets. The visualized clinical relevant heatmap based on 11-m6A regulated lncRNAs between low-risk and high-risk group is shown in Figure 10. With a difference in distribution between the low - and high-risk groups were the immune scores (P<0.05) and cluster analysis (P<0.001)(Figure 10A), which were higher in the low-risk group, and cluster analysis, which showed a denser distribution of Cluster 2 in the high-risk group. Patients with a low immune score, have an elevated risk and it follows that the immune score is favorable for prognosis. In terms of the distribution of clinicopathological characteristics and risk scores, patients with T3-T4 stage had an elevated risk compared with T1-T2 (P=0.016) (Figure 10B). PD-L1 expression was upregulated in the high-risk compared with the low-risk group (P=0.037).
We parallel the Kaplan-Meier curve of OS among different age, gender, clinicopathological grade, TNM stage. As shown in Figure11, the OS curves showed distribution discrepancy in terms of all clinical variables and demographic characteristics except M1 and stage III-IV variables. Patients in the low-risk group of these variables all had a superior overall survival than those in the high-risk group (Figure 11A-F) (P<0.01).
Association of m6A-related lncRNAs with immunocytes
We performed correlation analysis between the risk score and the 22 immune cells infiltrated in the TME, deleting samples with P < 0.05 for immune cells in TME. A total of 8 immune cells associated with the risk score were obtained, which were B cells naive, Plasma cells, CD8 T cells, T cells CD4 memory resting, Monocytes, Macrophages M0, Macrophages M1, and Dendritic cells activated. Positively correlated with the risk score were Dendritic cells activated, Macrophages M0, Macrophages M1, CD4 T cells (Figure 12A-D). However, B cells naive, plasma cells, CD8 T cells, Monocytes were negatively correlated with the risk score (Figure 12E-H).