m6A-related lncRNAs in GC
We defined m6A-related lncRNAs as those that were significantly correlated (P < 0.0001, |Cor|> 0.4) with m6A-related genes. Ultimately, 491 m6A-related lncRNAs were obtained. Furthermore, co-expression network of the m6A-related genes and lncRNAs was plotted (Fig. 2a). According to univariate Cox regression analysis (Fig. 2b), we obtained 23 m6A-related prognostic lncRNAs (P < 0.05) (Additional file 3: Table S3). Based on 32 normal samples and 375 GC samples from TCGA dataset, differential expression analysis of these m6A-related lncRNAs was made. Amidst 23 m6A-related lncRNAs of tumor samples, 15 lncRNAs exhibited significantly higher expression degrees (P < 0.05) while 8 lncRNAs demonstrated relatively lower expression levels (P < 0.05) (Fig. 2c, d) when compared with normal samples.
Cluster analysis of m6A-related prognostic lncRNAs
According to the expression profiles of m6A-related prognostic lncRNAs, we conducted unsupervised clustering to analyze the GC samples from TCGA dataset and divided samples into different subtypes. As exhibited in Fig. 3a, k = 2 was the most optimized selection. In order to further explore the relation between clustering result and clinical prognosis, we made survival analysis to compare the OS of GC patients between two subtypes. The consequence demonstrated that OS rate of cluster1 was inferior to that of cluster2 (P = 0.001) (Fig. 3b). Moreover, correlations between the cluster analysis and other clinical parameters, including age, gender, TNM stage, tumor stage and tumor grade, were plotted in the heatmap (Fig. 3c).
Functional enrichment analysis
In consideration of the favorable clustering result of OS for GC patients, we conducted GSEA between cluster1 and cluster2 to explore the potential biofunction of m6A-related lncRNAs. GSEA results indicated that several different tumor hallmarks were significantly enriched in two clusters, such as cell cycle, P53 signaling pathway, ECM receptor interaction and MAPK signaling pathway (P < 0.05) (Additional file 4: Figure S1a–d). Meanwhile, we found that certain immunity pathways were intimately associated with mast cells, Dendritic cells (DC), Natural Killer (NK) cells and T cells (P < 0.05) (Additional file 4: Figure S1e–h). Consequently, these results revealed that m6A-related prognostic lncRNAs were highly relevant to tumorigenesis and immune pathways.
Immune analysis of m6A related prognostic lncRNAs
Given the fact that several immune-related signaling pathways were enriched in both clusters, we further made an analysis on immunity which was comprised of immune checkpoints expression, TIC abundance profile and TME scores, so as to investigate the difference between the two clusters. First of all, we analyzed differentially expressed immune checkpoints, including PD‐1, PD‐L1, PD‐L2, IDO1, CTLA‐4, TIM‐3, LAG3 and TIGIT. Compared to cluster2, expression levels of both PD‐L2 and TIM-3 were significantly upregulated in cluster1 (P < 0.05) (Fig. 4a, b).
Subsequently, the co-expression analysis between immune checkpoints and m6A-related prognostic lncRNAs was implemented in R software. Amid 23 m6A-related prognostic lncRNAs, 16 lncRNAs was significantly correlated with PD‐L2 (7 positive correlations and 9 negative correlations; P < 0.05) while 15 lncRNAs was closely relevant to TIM (4 positive correlations and 11 negative correlations; P < 0.05) (Fig. 4c, d). In comparison to cluster2, cluster1 had more mast cells resting (P = 0.0003), monocytes (P = 0.0096) and T cells CD4 memory resting (P = 0.00044), but less macrophages M1 (P = 0.021), T cells CD4 memory activated (P = 0.0012) and T cells follicular helper (P = 1.1e-05) amidst the TICs with differential profiles (Fig. 4e, Additional file 5: Figure S2). Besides, Immune, stromal, and ESTIMATE scores, being relevant to the ratio of the immune/stromal components, were higher in cluster1 than cluster2 (P < 0.05) (Fig. 4f–h). Altogether, these results suggested that m6A-related prognostic lncRNAs were intimately associated with tumor immunity.
Construction of the m6A-LPS
To investigate the prognostic role of m6A-related lncRNAs in GC, we employed the Lasso algorithm to construct a m6A-LPS and further combined it with 23 m6A-related prognostic lncRNAs obtained from previous univariate Cox regression analysis. Based on the optimal value of λ (λ = 9), we ultimately screened out nine m6A-related prognostic lncRNAs (P < 0.05) (Fig. 5a, b). Taking advantage of both regression coefficients and expression levels of nine m6A-related prognostic lncRNAs (AC026691.1, Coefficient = 0.4785; AL139147.1, Coefficient = 0.4706; AL590705.3, Coefficient = 0.3874; TYMSOS, Coefficient = -0.0586; AL355574.1, Coefficient = -0.1085; AL390961.2, Coefficient = -0.2289; AC005586.1, Coefficient = -0.2724; and AP000873.4, Coefficient = -0.3635), we estimated the risk score of the m6A-LPS (Fig. 5c; Additional file 6: Table S4).
In order to appraise the prognostic role of m6A-LPS, we divided GC patients into train and test sets randomly (Additional file 7: Table S5). Furthermore, based on the median value of the risk score, we stratified GC patients into high- and low-risk groups. K-M survival curves demonstrated that GC patients in the high‐risk group had lower OS rates than their counterparts in both train and test sets (P < 0.001) (Fig. 5d, g). Afterwards, the ROC curves were plotted to evaluate the accuracy of m6A-LPS in predicting survival of GC at 1, 3, and 5 years. (Train set: AUC at 1, 3, 5 year is 0.705, 0.801, and 0.807, respectively; Test set: 0.682, 0.681 and 0.678, separately; Fig. 5e, h). Besides, we assessed risk scores of each GC case amidst train and test sets, which implied that GC patients in the low‐risk group had better survival status and shorter dead status than high-risk group (Fig. 5f, i). In general, above results indicated that m6A-LPS had a promising capacity to predict the survival of GC patients.
Independent prognostic analysis and stratification analysis
To determine whether the m6A-LPS could function as an independent prognostic indicator, we performed both univariate and multivariate Cox analysis on signature-based risk scores in train and test sets. The outcome of univariate Cox analysis suggested that m6A-LPS-based risk score was closely associated with unfavorable OS in both train set [hazard ratio (HR): 10.409, 95% CI 5.025–21.564, P < 0.001, Additional file 8: Figure S3a] and test set (HR: 4.490, 95% CI 1.823 − 11.061, P < 0.001, Additional file 8: Figure S3c). Moreover, the result of multivariate Cox regression analysis also demonstrated that corresponding risk score was intimately related to low OS in both train set (HR: 11.097, 95% CI 4.830 − 25.493, P = 0.001, Additional file 8: Figure S3b) and test set (HR: 6.411, 95% CI 2.394 − 17.172, P < 0.001, Additional file 8: Figure S3d). Meanwhile, age and stage were also verified to be closely relevant with the OS in univariate and multivariate Cox analysis. Generally speaking, these results confirmed that our m6A-LPS-based risk score might be an independent factor for prognostic evaluation.
Additionally, based on clinicopathological parameters containing age, gender, tumor stage and tumor grade, we employed stratification analysis to examine the predictive value on the survival of m6A-LPS in each section. The results demonstrated that the high-risk group had an intimate correlation with worse survival (P < 0.001) in both older and (≥ 65 years) younger group (< 65 years), male and female, advanced- (Stage III–IV) and early‐stage (Stage I–II) patients and tumor grade (G1–2 or G3) (Additional file 8: Figure S3e–l). The outcomes suggested that the m6A-LPS was adept in stably discriminating patients with undesirable prognosis. Besides, the overall correlation between risk score and clinical factors was plotted in Fig. 6a.
Correlation between the risk score and ICB therapy response
To further explore the underlying immune-related mechanism in GC patients, the differential expression of immune checkpoints analysis and TICs risk score were identified. Intriguingly, we found that expression levels of PD-1 and CTLA4 in the low-risk group were higher than in the high-risk one (Fig. 6b, c), which indicated that there were more immune escape and more protein expression of immune checkpoints in the low-risk group. On account of above findings, we further investigated the role of these lncRNAs in ICBs therapy. It was implied that the patients who responded to ICBs therapy had a lower risk score than patients who did not (Fig. 6d), suggesting that low-risk patients were underlying beneficiaries of ICBs therapy. Furthermore, we found that the risk score had significant negative correlations (P < 0.05) with B cells memory (R = − 0.18), Macrophages M0 (R = − 0.18), T cells CD4 memory activated (R = − 0.23), and T cells follicular helper (R = − 0.28). Meanwhile, the risk score was positively related to DCs resting (R = 0.21), Macrophages M2 (R = 0.28), Mast cells resting (R = 0.22), NK cells activated (R = 0.15), T cells CD4 memory resting (R = 0.24), and Monocytes (R = 0.31) (Additional file 9: Figure S4).
Validation of the m6A-LPS in GC tissues
We collected ten pairs of tumor tissues and para tumor tissues from Xijing Hospital. Via performing qRT-PCR, we found that two m6A-related lncRNAs (AC022031.2 and AL590705.3) were upregulated, while seven lncRNAs (TYMSOS, AC026691.1, AL355574.1, AP000873.4, AL390961.2, AC005586.1, and AL139147.1) were downregulated in cancer tissues (Additional file 10: Figure S5).
Reduced lncRNA AC026691.1 could promote proliferation and migration of GC cells
In view of the above findings, lncRNA AC026691.1 had the highest regression coefficients. To validate the tumorigenic role m6A-related lncRNA played in GC, siRNA was employed to silence the lncRNA AC026691.1 in SGC-7901 and BGC-803 cells. Moreover, we detected the transfection efficiency with qRT-PCR (Fig. 7a). As indicated in the results of CCK-8 assay, via silencing with siRNA, the ability of cell proliferation was dramatically enhanced in SGC-7901 and BGC-803 cells (Fig. 7b). Additionally, outcomes of the wound-healing assay demonstrated that the migration capability was significantly promoted (Fig. 7c–f). Taken above results together, it was natural to summarize that lncRNA AC026691.1 could play a vital role in suppressing GC.
lncRNA AC026691.1 had a strong interaction with FTO in m6A modification of GC
To explore the underlying mechanism of lncRNA AC026691.1 in m6A modification, we made correlation analysis between AC026691.1 and m6A-related genes. Corresponding result suggested that lncRNA AC026691.1 had the most positive relation with FTO (Coefficient = 0.5080, P = 5.43e−26) (Additional file 11: Table S6). After silencing lncRNA AC026691.1, the expression level of FTO comparatively decreased in SGC-7901 and BGC-803 cells (P < 0.05) (Fig. 8a, b). Additionally, the m6A level exhibited a downward tendency in GC cell lines after downregulating the lncRNA AC026691.1 (P < 0.05) (Fig. 8c, d). Furthermore, by downregulation of FTO in GC cell lines, the expression level of LncRNA AC026691.1 dramatically declined in qRT-PCR results (P < 0.05) (Fig. 8e). And the underlying regulatory site of FTO was further predicted and demonstrated in Additional file 12: Figure S6. In general, the above results indicated that the lncRNA AC026691.1 interacted closely with FTO, via regulation of which, the expression level of m6A got reduced.