The RNA transcriptome dataset and the corresponding GC clinical information were downloaded from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) database. Genes were grouped into protein coding genes and lncRNA genes based on the human genome annotation data. Expression levels of 13 m6A regulator genes were also determined. We used the Pearson correlation coefficient to evaluate the correlation between m6A regulator genes and lncRNA. The lncRNA with an absolute correlation coefficient >0.3 and a P value <0.05 were considered as a m6A-related lncRNA. Patients were then grouped into two groups: the training dataset and the testing dataset. Retrieved data were used for subsequent bioinformatics analysis.

Univariate cox regression analysis and least absolute shrinkage and selection operator (LASSO)-penalized Cox regression analysis were employed to identify the m6A-related lncRNA prognostic signature in the training dataset. The risk score for each GC patient was calculated based on the following formula:

Risk score = Σ_{Expi *βi}. where Exp_i represents each lncRNA expression and β_i represents the coefficient of each lncRNA. The receiver operating characteristic (ROC) curve analysis was used to evaluate the accuracy of the lncRNA signature in the training dataset, testing dataset, the complete dataset, and the DFS dataset.

The complete dataset and DFS dataset with corresponding clinical information were used for subsequent analysis. To identify the independence of the lncRNA signature, we performed a univariate cox regression and multivariate cox regression analysis between the lncRNA signature and the clinical traits, and the significant independent prognostic factors was selected based on a P value <0.05.

The nomogram prediction model was constructed based on the lncRNA signature risk score and independent clinical factors, using the “rms” R package. The calibration plot was applied to validate the calibration and accuracy of the nomogram (by a bootstrap method with 1,000 replicates).

The Gene Set Enrichment Analyses (GSEA) method was used to explore the potential KEGG pathway implicated in the lncRNA prognostic signature. The reference gene set was retrieved from c2.cp.kegg.v7.1.symbols files, and the significant pathways were screened based on the criterion: P <0.05 and FDR <0.25.

We used non-negative matrix factorization (NMF) clustering methods to explore potential molecular subgroups based on the prognostic lncRNA expression profile. The optimal K cluster was selected on the basis of the cophenetic correlation coefficient. We further used the TIDE algorithm to predict the response of the four subgroups to immunotherapy.

The differentially expression genes (DEGs) between the high-risk group and the low-risk group were identified using the “limma” R package. We then uploaded the top 1,000 DEGs to the CMAP database to identify which target compounds might be useful.

Computational and statistical analyses were performed using R software (https://www.r-project.org/). The Kruskal-Wallis test was used to compare the four subgroups and to explore their response to immunotherapy. Differences between the high-risk group and the low-risk group were determined using the Kaplan-Meier curve and log-rank test. Clinical data were analyzed using the chi-square test or Fisher’s exact test. For all analyses, a P value <0.05 was considered statistically significant.

m6A RNA methylation regulators play a crucial role in progression of malignant tumors; therefore, we explored the expression profiles of 13 m6A RNA methylation regulators in GC. High expression levels of the 13 genes were observed in the tumor tissue compared with normal tissue ( Figure S1A ). Notably, the tumor tissue showed significantly higher expression levels of m6A regulators, except for FTO and ALKBH5, compared with the normal tissue with a P value <0.05. Further, the correlation analysis showed that HNRNPC and YDHDF2 methylation regulators were highly correlated with GC progression (Cor = 0.57) ( Figure S1B ).

To identify potential m6A-related lncRNAs, the Pearson correlation coefficient was used to assess the relationship between m6A methylation regulators and lncRNAs. As a result, we identified 1,351 interactions and 800 m6A-related lncRNAs with an absolute correlation coefficient >0.3 and a P value <0.05. Among them, 117 interactions were identified between lncRNAs and m6A methylation regulators ( Table S1 ).

A total of 339 patient samples were grouped into a training dataset (N = 136) and a testing dataset (N =203). We then performed a univariate cox regression analysis and LASSO-penalized Cox regression analysis to construct a 11-lncRNA signature model in the training dataset ( Figure 1 ). The risk score for each patient in the training dataset, testing dataset, and the complete dataset was calculated based on the risk formula: risk score = AL049840.3 * 0.599866058 + AC008770.3 * (-1.237087957) + AL355312.3 * (-0.19130367) + AC108693.2 * (-0.956067535) + BACE1-AS * (-0.362760192) + AP001528.1 * 0.528553101 + AP001033.2 * 0.594102051 + AC092574.1 * (-0.618599189) + AC010719.1 * (-0.337936563) + AC009090.3 * 0.770122519 + SAMD12-AS1 * (-0.766369919). Patients in the training dataset, testing dataset, and complete dataset were further grouped into a high-risk group and low-risk group based on the median risk score. We found that patients in the high-risk group corresponded to a greater number of deaths than the low-risk group patients in the training dataset and testing dataset, respectively ( Figure S2 ). The Kaplan-Meier (KM) curve analysis results showed that the low-risk group have a better prognosis than the high-risk group in the testing, training, and the complete dataset (P <0.05) ( Figure 2 ). In addition, the area under the curve (AUC) for 5-year overall survival (OS) was 0.94, 0.85, and 0.89 in the training dataset, testing dataset, and complete dataset, respectively, implying that the lncRNA signature has a good accuracy in the prognostic prediction of GC ( Figure 3 ).

The association between the lncRNA signature and clinical factors were evaluated by performing a KM curves analysis. As shown in Figure 4A , we found that the risk score significantly increased from the younger group to older group (P <0.05). Moreover, the stage also presented a significant divergence from stage I to stage IV (P < 0.05) ( Figure 4B ). These results indicate that the greater the progression is, the higher the risk score of GC is. Moreover, we also investigated the prognostic value of the m6A-reated lncRNA signature in GC patients stratified by clinical pathological variables, including age, gender, cancer grade, and stage. For all different stratifications, patients in the high-risk group tended to have a lower overall survival rate compared to the low-risk group ( Figures 5A–H ). These results suggest that the m6A-related lncRNA signature can predict the prognosis of GC without the need to consider clinical factors.

To determine the independence of the lncRNA signature in the clinical factors, we conducted a univariate cox regression analysis and multivariate cox regression analysis between the lncRNA signature and the clinical factors in the complete dataset. As shown in Figures 6A, B , we found that the lncRNA signature can act as an independent prognostic factor for the prognosis of GC.

The potential pathways or functions of the m6A-related lncRNA signature were explored by performing a Gene Set Enrichment analysis. As shown in Figures 7A–C , we discovered that patients in the high-risk group were mainly involved in the cytokine receptor interaction, focal adhesion, and ECM receptor interaction pathway, while basal transcription factors pathway were enriched in the low-risk group ( Figure 7D ).

Considering the significance of the disease-free survival in the prognosis of GC, we further calculated the risk score for each GC patient based on the previous risk formula. The KM curve analysis result suggests that patients in the high-risk group have a poor survival compared to the low-risk group ( Figure S3A ). The ROC curve analysis results suggest that the lncRNA signature in the DFS have a good accuracy of prognostic prediction ( Figure S3B ). Moreover, the univariate cox regression and multivariate cox regression analysis suggest that the lncRNA signature can serve as an independent prognostic factor for the prognostic prediction of GC ( Figures 8A, B ). Moreover, we also compared our risk model with other reported lncRNA signatures, such as Wu et al. (15) and Chen et al. (16). As shown in Figure S4 , we observed that our risk signature has the highest AUC value when compared with other risk models, suggesting that our signature is a reliable prognostic model.

To further establish an individualized prediction model in the OS and DFS, we further constructed a nomogram for GC patients based on the independent prognostic factors, respectively ( Figures 9A , 10A ). The calibration plots show that the performance of the nomogram have a good concordance with the prediction of 1-, 3-, and 5-year OS ( Figures 9B–D ) and DFS ( Figures 10B–D ).

To explore the potential molecular drugs of GC, we first performed a differentially expression analysis between the high-risk and low risk group. As a result, a total of 150 DEGs were identified and further uploaded to the CMAP drug database. In total, 64 drugs with 45 mechanisms of action (MOA) were shared by the above drugs ( Figure 11 and Table S2 ). According to the drug result, metixene, dicycloverine, cyclopentolate, and procyclidine shared the mechanism of Acetylcholine receptor antagonist; dobutamine, orciprenaline, and etilefrine shared the mechanism of Adrenergic receptor agonist; nimesulide, nabumetone, and LM-1685 shared the mechanism of cyclooxygenase inhibitor; thioperamide, doxepin, and doxylamine shared the mechanism of Histamine receptor antagonist. Our study identified drugs targeting the m6A-related lncRNA signature and might provide therapeutic targets for further analysis.

The NMF algorithm, using the m6A-related lncRNAs, was used to explore the molecular subtypes of GC. To ensure a robust and reliable subtype, low expression level lncRNAs were filtered out and retained the lncRNAs that were associated with survival of GC for univariate cox regression analysis, to ensure a robust and reliable subtype. The data was then used in a NMF clustering analysis. The cophenetic correlation coefficients were calculated to determine the optimal k value, and k = 4 was eventually selected as the optimal cutoff after comprehensive consideration ( Figure 12A , the four subtypes were named C1, C2, C3, and C4). When k = 4, the consensus matrix heatmap still presented sharp and crisp boundaries, suggesting a stable and robust clustering for the samples ( Figure 12B ). Notably, the KM curve analysis result showed a significant survival difference among subgroups, of which C1 and C2 showed a better survival outcome compared to C3 and C4 (P = 0.031) ( Figure 12C ). Moreover, the different expression patterns of the lncRNAs in the four subgroups suggest that these lncRNAs might play an important role in the subgroups ( Figure 12D ). Further, we used the TIDE algorithm to predict the response of the four subtypes to immunotherapy. Subtype C1 and C2 responded better to immunotherapy compared with subtype C3 and C4 (P = 0.0072). In addition, we also evaluated the SNP alteration among these subtypes. We observed that subtype C3 has the highest SNP alteration, while C1 corresponded to the lowest SNP alteration ( Figure 13 ).