Identification of N6-Methylandenosine-Related lncRNAs for Subtype Identification and Risk Stratification in Gastric Adenocarcinoma

Objectives The purpose of this study was to investigate the role of m6A-related lncRNAs in gastric adenocarcinoma (STAD) and to determine their prognostic value. Methods Gene expression and clinicopathological data were obtained from The Cancer Genome Atlas (TCGA) database. Correlation analysis and univariate Cox regression analysis were conducted to identify m6A-related prognostic lncRNAs. Subsequently, different clusters of patients with STAD were identified via consensus clustering analysis, and a prognostic signature was established by least absolute shrinkage and selection operator (LASSO) Cox regression analyses. The clinicopathological characteristics, tumor microenvironment (TME), immune checkpoint genes (ICGs) expression, and the response to immune checkpoint inhibitors (ICIs) in different clusters and subgroups were explored. The prognostic value of the prognostic signature was evaluated using the Kaplan-Meier method, receiver operating characteristic curves, and univariate and multivariate regression analyses. Additionally, Gene Set Enrichment Analysis (GSEA), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and Gene Ontology (GO) analysis were performed for biological functional analysis. Results Two clusters based on 19 m6A-related lncRNAs were identified, and a prognostic signature comprising 14 m6A-related lncRNAs was constructed, which had significant value in predicting the OS of patients with STAD, clinicopathological characteristics, TME, ICGs expression, and the response to ICIs. Biological processes and pathways associated with cancer and immune response were identified. Conclusions We revealed the role and prognostic value of m6A-related lncRNAs in STAD. Together, our finding refreshed the understanding of m6A-related lncRNAs and provided novel insights to identify predictive biomarkers and immunotherapy targets for STAD.


INTRODUCTION
Globally, gastric cancer (GC) is the fifth most common cancer and the third most deadly neoplasm (1). Gastric adenocarcinoma (STAD) is the most common pathological type of GC, and despite considerable progress in the diagnosis and therapeutic strategies for STAD, the prognosis of patients with STAD remains poor due to advanced stage and postsurgical recurrence (2,3). Therefore, the identification of novel biomarkers for early detection and effective therapeutic targets for treating patients with STAD is critical and urgent.
Accumulating evidence has shown that long noncoding RNAs (lncRNAs) had various biological functions and played a crucial role in the oncogenesis and progression of GC (4). For example, lncRNA IGF2-AS functions as a competing endogenous RNA (ceRNA) to miR-503 and promotes the pathogenicity of STAD by regulating SHOX2 (5). LINC00707 acts as an oncogene in GC by interacting with RNA-binding protein HuR and increasing the stability of VAV3/F11R mRNAs (6). LncRNA CRNDE could bind to splicing protein SRSF6 to reduce its stability and thus regulate alternative splicing events and affect autophagy regulation in GC (7).
Increasing evidence suggests that RNA modifications play a critical role in tumorigenesis and progression of different cancers, including GC (8,9). N6-methyladenosine (m 6 A), which introduces a methyl group in the nitrogen-6 position of adenosine, is found to be the most frequent internal RNA modification in mammals (10). As a dynamic and reversible process, m 6 A RNA modification is primarily regulated by "writers" (adenosine methyltransferases) and "erasers" (demethylases) and performs different functions by interacting with "readers" (m 6 A-binding proteins). As identified to distribute extensively in a variety of RNAs, such as messenger RNAs (mRNAs), pri-microRNAs (pri-miRNAs), circular RNAs (circRNAs), and lncRNAs, m 6 A is involved in various biological processes related to the occurrence and progression of tumors, including GC (11)(12)(13). For instance, the m 6 A writer METTL3 stimulates m 6 A modification of HDGF mRNA, and the m 6 A reader IGF2BP3 recognizes and binds to the m 6 A site and enhances its stability, which promotes tumor angiogenesis and glycolysis in GC (14). Overexpression METTL3 facilitates the processing of pri-miR-17 into the miR-17 through an m 6 A DGCR8-dependent method, which activates the AKT/mTOR pathway and the progression of GC (15). LINC00470 promotes the degradation of PTEN mRNA to facilitate malignant behavior in GC cells by interacting with METTL3 (16). Additionally, extensive literature has demonstrated that m 6 A plays an important role in immune recognition, immune responses, and tumor microenvironment (TME) (17)(18)(19). However, the specific role and prognostic value of m 6 A-related lncRNAs in STAD remain unclear.
Here, we analyzed the Cancer Genome Atlas (TCGA) database for m 6 A-related lncRNAs involved in STAD, identified two clusters based on m 6 A-related lncRNAs, and constructed an m 6 A-related lncRNA prognostic signature. Then, we estimated its predictive value and diagnostic effectiveness, as well as the correlation of m 6 Arelated lncRNAs with TME and immunotherapy. Furthermore, the molecular mechanisms associated with m 6 A-related prognostic lncRNAs were explored. The finding in this study revealed the critical role of m 6 A-related lncRNAs and shed light on the latent relationship and the underlying mechanism between m 6 A-related lncRNAs and tumor-immune interactions.

Acquisition of Datasets
The RNA-seq transcriptome data [fragments per kilobase million (FPKM)] (20) from 373 samples and clinical information from 406 patients with STAD in the TCGA database (http://cancergenome.nih.gov/) were downloaded for our study. Patients with complete clinicopathological and survival information were included for further assessment. different clinicopathological characteristics and OS were compared between clusters I and II. Furthermore, the differences in the content of immune infiltrating cells, TME scores, and the expression of immune checkpoint genes (ICGs) between different clusters were explored (25). The content of immune infiltrating cells was identified by CIBERSORT (26), and the immune/stromal scores and tumor purity were calculated through the "ESTIMATE" package in R (27).
Then, we randomly divided the patients with STAD into two groups: the training group and the testing group. Subsequently, based on m 6 A-related prognostic lncRNAs identified by univariate Cox regression analysis, the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm was used to identify m 6 A-related lncRNAs with powerful prognostic significance and construct the prognostic risk model from the training group data. According to the best penalty parameter l, the coefficients of the m 6 A-related lncRNAs were calculated. The risk score (RS) was estimated using the following formula: where Coef(i) is the coefficient and X(i) represents the expression levels of m 6 A-related lncRNAs. Using the median RS obtained as the demarcation value, patients with STAD were classified into two groups: high-risk and low-risk groups. Kaplan-Meier analysis and the receiver operating characteristic (ROC) curves were used to validate the predictive efficiency (28). Then, the accuracy of the model was validated from the test group and the combined group using the same method. Furthermore, the differences in clinicopathological features, the content of immune infiltrating cells, TME scores, and ICGs expression between high-risk and low-risk groups were also explored. Moreover, we further predicted the response to immune checkpoint inhibitors (ICIs) in subgroups based on the immunophenoscores (IPS) of patients with STAD obtained from The Cancer Immunome Atlas (TCIA) (https://tcia.at/ home). Additionally, the prognostic value of the RS was verified using univariate and multivariate Cox regression analyses. The hazard ratio (HR) with 95% confidence intervals and log-rank p-value were calculated using the "glmnet" and "survival" R packages (29).
To explore the biological functions associated with m 6 Arelated lncRNAs, Gene Set Enrichment Analysis (GSEA), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and Gene Ontology (GO) analysis were performed. Genes that were significantly upregulated (fold change >1 and p < 0.05) or downregulated (fold change <−1 and p < 0.05) between clusters I and II or between the high-risk and low-risk groups were identified using the "edgeR" package in R, which were used for GO and KEGG pathway analysis. Additionally, genes in different clusters and different risk groups were functionally annotated using GSEA. Based on the m 6 A-related prognostic lncRNAs, the target miRNAs were predicted via miRcode database and target mRNAs of these miRNAs were found in different databases, such as TargetScan, miRTarBase, and miRDB. Target mRNAs in the ceRNA network were also functionally annotated using GO and KEGG pathway analyses. The flow chart of bioinformatic analysis was shown in Figure 1.

Cell Culture
The GC cell line MGC-803 and the normal human gastric epithelial cell line GES-1 were purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA). All cells were cultured in RPMI-1640 medium (Life Technologies, Grand Island, NY, USA) supplemented with 10% fetal bovine serum (Life Technologies) at 37°C in a humidified atmosphere with 5% CO 2 .

Statistical Analysis
The expression data of m 6 A-related regulators and all lncRNAs in tumor tissues and adjacent mucosa of STAD obtained from TCGA was compared using one-way analysis of variance (ANOVA); the clinical characteristics of different groups were compared using the Chi-square test; the Kaplan-Meier method was used to perform a bilateral logarithmic rank test in overall survival analysis; p-value <0.05 were regarded as statistically significant. All statistical analyses were implemented using R v4.0.3 (https://www.r-project.org/) or GraphPad Prism software (Version 8.0).

Identification of m 6 A-Related Prognostic lncRNAs
Firstly, the expression levels of 24 m 6 A-related genes and all lncRNAs from the TCGA were extracted respectively. Through coexpression analysis, we identified 471 m 6 A-related lncRNAs (|cor| > 0.4, p-value <0.05). The gene co-expression network of 24 m 6 A-related genes and 471 m 6 A-related lncRNAs was shown in Figure 2A. After conducting univariate Cox analysis, 19 candidate lncRNAs that were highly correlated with OS were identified (p < 0.05) ( Figure 2B). The expression of 19 m 6 A-related prognostic lncRNAs was compared between tumor tissues and adjacent mucosa ( Figures 2C, D). Among these lncRNAs, seven (AL139147.1, AC022031.

Consensus Clustering of m 6 A-Related Prognostic lncRNAs Identified Two Clusters of STAD With Different Clinicopathological Features and Immune Landscape
Based on the expression levels of 19 m 6 A-related prognostic lncRNAs, consistent clustering analysis of patients with STAD was implemented. Patients were clustered into two clusters due to the minimal interference between the two subgroups ( Figures 3A-D).
The distribution of the clinicopathological characteristics in clusters I and II were displayed as a heat map ( Figure 3E). Evident differences between the two clusters according to tumor grade (p < 0.01) were observed. Notably, the OS rate of clusters I and II were significantly different based on the Kaplan-Meier method and cluster II was associated with poorer OS ( Figure 3F).
In terms of TME scores, it was found that ESTIMATE, immune, and stromal scores significantly increased in cluster II, while tumor purity increased in cluster I. In addition, the content of 22 immune cells in clusters I and II was compared ( Figure 4A). As a result, cluster I contained more follicular helper T cells (p < 0.001) and M0 macrophages (p < 0.01), and cluster II had more monocytes (p < 0.01), M2 macrophages (p < 0.05), resting dendritic cells (DC) (p < 0.001) and resting Regarding the expression of ICGs, we investigated the distribution of 38 ICGs obtained from previous studies in different clusters (30)(31)(32)(33)(34). Through differential expression analysis, we found that 22 ICGs in the tumor tissues differentially expressed compared with the adjacent mucosa ( Figure 5A), and 11 of 22 ICGs in cluster I differentially expressed compared with cluster II ( Figure 5B). Moreover, we observed three ICGs (namely, YTHDF1, IL23A, LDHC) were significantly overexpressed in cluster I and eight ICGs (namely, PDCD1LG2, CD86, HAVCR2, LAMA3, TNFSF4, IL12B, LDHA, ICOS) were significantly overexpressed in cluster II.

Construction and Verification of the m 6 A-Related lncRNAs Prognostic Signature
Based on 19 candidate lncRNAs that were highly correlated with OS, we used the LASSO method in the training group to construct an m 6 A-related lncRNA signature for evaluating the prognosis of patients with STAD. Finally, 14 lncRNAs were chosen to establish a prognostic signature and the risk score was calculated ( Figures 6A, B). Using the median risk score as the demarcation value, the patients in the training group (n = 170) were classified into two groups, namely, the high-risk and lowrisk groups. To test the efficacy of the prognostic model, survival and ROC curve analyses were conducted. Kaplan-Meier analysis showed that the low-risk group had a significantly longer survival time than the high-risk group (p < 0.001) ( Figure 6C). The value of the area under the curve (AUC) in the time-dependent ROC curve of 1, 3, and 5 years was 0.718, 0.808, and 0.873 severally ( Figure 6D), suggesting good prediction performance of the survival model. To further validate this 14 lncRNA prognostic signature, verification analysis in the test group (n = 168) and the combined group (n = 338) was implemented. As a result, the high-risk group in the test group and the combined group had significantly shorter survival time compared with the low-risk group, which was previously observed in the training group ( Figure 6C). The time-dependent ROC curve of the test group and the combined group also had well-prediction performances, and the AUC value of 1, 3, and 5 years is shown in Figure 6D.
The distribution plot of the risk score and survival status showed that the higher the risk score, the more deaths of patients with STAD ( Figures 6E, F).The expression of 14 m 6 A-related prognostic lncRNAs in the training group, the test group and the combined group, was displayed as a heat map ( Figures 6G).
To examine whether the risk score was an independent prognostic factor, univariate and multivariate Cox regression analyses were conducted. In the training group, the risk score was significantly associated with OS both in univariate and multivariate Cox regression analyses (p < 0.001), in addition to age at diagnosis and pathological stage (p < 0.05). It should be noted that the risk score was also closely related with OS in the test group (p < 0.05) and the combined group (p < 0.001) by the same analysis, which indicated that the risk score was an independent powerful prognostic factor for the prognosis of OS in STAD (Figures 7A, B).

Subgroup Analysis With Different Clinicopathological Features and Immune Landscape
The expression of 14 m 6 A-related prognostic lncRNAs and the distribution of clinicopathological characteristics, the immune scores of TME, and the clustering of patients in the high-risk and low-risk groups were displayed as a heat map ( Figure 7C). Evident differences between the two groups according to different clusters (p < 0.001) were observed. Significant differences of risk score were found between (1): different tumor grades (p < 0.01) (2), different immune scores of TME (p < 0.05), and (3) different clusters (p < 0.01) ( Figure 7D).
To evaluate whether the m 6 A-related lncRNAs prognostic model could serve as a prognostic indicator for OS in subgroups of patients with different clinical characteristics, we stratified subgroups by age (age ≤65 and age >65), gender (female and male), grade (G1-G2 and G3), clinical stage (stages I-II and III-IV), stage T (T1-T2 and T3-T4), stage M (M0 and M1), and stage N (N0 and N1). As the result shown in Figure 8, the OS of the lowrisk patients based on age (p = 0.004 in age ≤65 and p < 0.001 in age >65), sex (p = 0.004 in female and p < 0.001 in male), grade (p = 0.006 in G1-G2 and p < 0.001 in G3), clinical stage (p = 0.017 in stages I-II and p < 0.001 in stages III-IV), stage T3-T4 (p < 0.001), stage M0 (p < 0.001), and stage N (p = 0.014 in N0 and p < 0.001 in N1) was significantly higher than those of the high-risk patients.
Furthermore, we analyzed the association between immune cells and ICGs and subgroups. As a result, we found that the high-risk group had significant positive correlations with infiltrating levels of resting DC (r = 0.17, p = 0.025), eosinophils (r = 0. 23 (Figures 9G-J). Next, we observed that four ICGs (namely, JAK2, LAMA3, PDCD1LG2, JAK1) were significantly high expression in high-risk group and six ICGs (namely, YTHDF1, TNFRSF18, TNFRSF4, PVR, PDCD1, IFNG) were significantly high expression in low-risk group ( Figure 5C). Furthermore, the difference in the response of ICIs between highrisk and low-risk groups was explored. As a result, the low-risk group was more likely to respond to immunotherapy than the highrisk group whether they were compared in IPS-PD1 (

Construction of the ceRNA Network and Functional Enrichment Analysis
To explore the biological function of 19 m 6 A-related prognostic lncRNAs, a ceRNA network was constructed based on the mechanism of lncRNAs regulating mRNAs expression by sponging miRNAs. Two lncRNAs were extracted from the miRcode database and 28 pairs of interaction between the Two lncRNAs and 24 miRNAs were identified. Based on three mRNA predicting databases mentioned previously and differentially expressed mRNA between the normal and tumor groups of STAD in TCGA, we identified 70 target mRNA. Finally, two lncRNAs, 24 miRNAs, and 70 mRNAs were included to construct ceRNA by Cytoscape software 3.7.1 ( Figure 11A). KEGG pathways and GO analysis were performed to annotate the function of 70 target mRNAs, and we found that these target mRNAs were enriched in transcription coactivator activity and tubulin binding (GO analysis); microRNAs in cancer, MAPK signaling pathway, proteoglycans in cancer, PI3K-Akt signaling pathway, focal adhesion, Rap1 signaling pathway, and Ras signaling pathway (KEGG pathways) (Figures 11B, C).
As we stratified the patients with STAD into clusters I and II or high-risk and low-risk groups, genes that were significantly upregulated (fold change >1 and p < 0.05) or downregulated (fold change <−1 and p < 0.05) between the different clusters or different subgroups were identified using the "edgeR" package in R. GO and KEGG pathways analysis were used for biological functional analysis. Concerning the differentially expressed genes between different clusters, these genes were associated with immune-related biological processes, such as "antigen binding" and "immunoglobulin receptor binding," and malignancyassociated pathways, including "focal adhesion" and "Wnt signaling pathway." The differentially expressed genes between the high-risk and the low-risk groups were enriched in malignancy-associated biological processes and pathways, containing "chemokine receptor binding", "chemokine activity (GO analysis)", and "focal adhesion" and "Wnt signaling pathway (KEGG pathway) (Figures 11B, C)".
Furthermore, we used GSEA to predict the functional difference between clusters I and II or between the high-risk and the low-risk groups. The results showed that cluster II had a worse OS and a lower 5-year survival rate was closely related with (C) The heat map showed the expression of 14 m 6 A-related prognostic lncRNAs, the distribution of clinicopathological characteristics, the immune scores of TME, and the clustering of patients in the high-risk and low-risk groups. ***p < 0.001. (D) The risk score was significantly higher in patients with higher-grade tumors (p < 0.01), higher immune scores of TME (p < 0.05), and cluster II (p < 0.01).
the malignant hallmarks of cancer, such as "focal adhesion", "ECM-receptor interaction", and "cell adhesion molecules (CAMs)" and immune-related pathway, such as, "complement and coagulation cascades". The high-risk group was also associated with cancer-related pathways, which were similar to cluster II (Figure 12).

Validation of the Expression Levels of the m 6 A-Related lncRNA in Cell Lines
For validating the expression levels of the m 6 A-related prognostic lncRNAs from prognostic signature, we detected six m 6 A-related prognostic lncRNA expression levels in GC cell line MGC-803 and normal human gastric epithelial cell line GES-1. Our results showed that LINC01537 and MAGI2-AS3 were significantly downregulated in MGC-803 compared with GES-1 while LINC00106 and STARD4-AS1 were significantly upregulated. SREBF2-AS1 was highly expressed in MGC-803 and SENCR was lowly expressed in MGC-803, but they had no significant difference in cells ( Figure 13).

DISCUSSION
Recently, an increasing number of studies focusing on the role of lncRNAs in GC proved that lncRNAs exerted a critical oncogenic role based on their dysregulated expression and localization (35). Additionally, as the most abundant posttranscriptional modification in eukaryotic noncoding RNAs (ncRNAs), m 6 A has a huge effect on its stability and transport (36)(37)(38). Previous studies have shown that m 6 A "writers" and "erasers" could adjust the levels of m 6 A modification in mRNAs and ncRNAs to regulate binding sites to m 6 A "reader" proteins. Different m 6 A reader proteins recognize and bind to methylated ncRNAs to perform different functions. For instance, YTHDF3 recognizes and binds to m 6 A-modified lncRNA GAS5 and promotes its degradation, which elevates YAP expression and promotes colorectal cancer (CRC) progression (39). The m 6 A mark increases the stability of lncRNA FAM225A, which promotes the progression of nasopharyngeal carcinoma by acting as ceRNA to sponge miR-590-3p/miR-1275 (40). IGF2BP2 recognizes and binds to m 6 A-modified circRNA NSUN2 and increases its export to the cytoplasm (41). Overexpression METTL3 can significantly increase the nuclear localization of lncRNA RP11 in CRC cells (42). Additionally, copious studies have shown that m 6 A RNA modification affects the maturation and response function of immune cells in tumor immunity and could remodel TME (43,44). For example, durable neoantigen-specific immunity is regulated by m 6 A-modified mRNA through YTHDF1 (45). METTL3mediated m 6 A of CD40, CD80, and TLR4 signaling adaptor Tirap transcripts could enhance their translation in DC for stimulating T-cell activation (46). Moreover, some studies have indicated that m 6 A-modified lncRNAs participate in the immune response. As DC migration is critical for the protective immunity and immune homeostasis, m 6 A modification promotes the degradation of lncRNA Dpf3, which is associated with DC migration (47). Thus, in consideration of the crucial role of lncRNAs, m 6 A RNA modification, and immunity response in STAD, these researches call our attention to investigate the gene profile of m 6 A-related lncRNAs in STAD, explore whether m 6 Arelated lncRNAs could serve as ideal biomarkers for STAD prognosis and participate in STAD initiation, progression, TME, and immunotherapy. In our study, a total of 373 samples, 406 patients with STAD, 24 m 6 A-related regulators, and 13,162 lncRNAs were included to exploit the specific role of m 6 A-related lncRNAs in STAD. Nineteen candidate lncRNAs that were highly correlated with OS were identified. Based on the expression of 19 m 6 A-related prognostic lncRNAs, patients with STAD were clustered into two clusters (clusters I and II), which had significant differences in the OS rate and tumor grade. Then, 14 of 19 m 6 A-related prognostic lncRNAs were used to establish a prognostic signature with the LASSO method in the training group.
Kaplan-Meier analysis showed that the OS of patients with low-risk scores was longer than those of patients with high-risk scores. Additionally, the result of ROC curve analysis indicated that the 14-lncRNAs signature could serve as a highly specific and sensitive prognostic survival model in STAD. Moreover, the results were further validated in the test groups and in the combined group. This signature can be used as an independent prognostic factor for STAD, suggesting that these 14 lncRNAs may be vital m 6 A-related lncRNAs and significant prognostic factors for patients with STAD. Furthermore, this m 6 A-related lncRNA prognostic model could serve as a prognostic indicator for OS in subgroups of patients with different clinical characteristics, especially age, sex, tumor grade, clinical stage, stages T3-T4, stage M0, and stage N. Moreover, different clusters or subgroups were correlated with diverse TME scores, immune infiltration cells, and ICGs A B FIGURE 12 | GSEA analysis of the cluster II and the high-risk group. (A) The cluster II was mainly enriched in "focal adhesion", "ECM-receptor interaction", and "cell adhesion molecules (CAMs)", and "complement and coagulation cascades". (B) The high-risk group was mainly enriched in "focal adhesion" and "ECM-receptor interaction". expression, which suggested that m 6 A-related lncRNAs may play a critical regulatory role in the TME and immune exhaustion of STAD. The tumor purity of cluster II that was associated with poorer OS and the high-risk group decreased, which indicated the promotion role of stromal cells and the dysfunction of immune infiltration cells in STAD. Fibroblasts, the most abundant stromal cell, could promote tumor proliferation and metastasis by secreting cytokines, growth factors, and chemokines and the remodeling of extracellular matrix (ECM) (48). Concerning immune infiltration cells, there is a significant positive correlation between M2 macrophages and cluster II and the high-risk group, which implicated the potential regulatory role of m 6 A-related lncRNAs in the polarization of tumor-associated macrophages. M2 macrophages could facilitate tumor occurrence and metastasis through the inhibition of T-cell-mediated antitumor immune response and the promotion of tumor angiogenesis (49). Although abundant immune cells infiltrated in TME, they could not lead to the antitumor immune response because of immune checkpoints blockade and the effects of hypoxia and metabolic alternation. To be more specific, hypoxia is strongly linked to immunosuppression via immune trafficking, angiogenesis, and alteration of molecular markers, which has hugely detrimental effects on T-cell function (50). Similarly, metabolic alternation also can lead to T-cell suppression, which was related to a host of metabolic byproducts (51). Moreover, as our results showed, some prominent ICGs, such as PDCD1LG2, HAVCR2, and ICOS, were significantly overexpressed in cluster II and the high-risk group, which resulted in immune exhaustion in TME. In addition, we discovered the difference in the response to ICIs between subgroups. Together, these findings suggest that the m 6 A-related lncRNAs may play a potential influence on the dysfunction of immune infiltration cells and the response to immunotherapy in STAD.
To provide a comprehensive analysis of m 6 A-related lncRNAs, a ceRNA network consisting of two lncRNAs, 24 miRNAs, and 70 mRNAs was constructed, and differentially expressed genes between different clusters or subgroups were identified for viewing the latent functions of m 6 A-related lncRNAs. With the two key m 6 A-related lncRNAs found in STAD, lncRNA MAGI2-AS3 has been reported to be associated with patient clinical characteristics in different tumors and can regulate a variety of biological processes (52). MAGI2-AS3 significantly promoted GC progression and migration via maintaining ZEB1 overexpression by sponging miR-141/200a (53). KEGG pathway and GO analysis shows that 70 target mRNAs and differentially expressed genes between different clusters or subgroups were enriched in several biological processes and pathways associated with the occurrence and progression of STAD (54)(55)(56), including "MAPK signaling pathway," "PI3K-Akt signaling pathway," "Wnt signaling pathway," "focal adhesion," and so on. Genes in different clusters or subgroups were functionally annotated using GSEA. Genes in cluster II and the high-risk group were also enriched in the cancer-related pathways. In addition, it should be noted that several biological processes and pathways associated with immune response were identified, such as "antigen binding", "immunoglobulin receptor binding", and "complement and coagulation cascades". Previous studies have demonstrated that the activation of the complement system and coagulation cascades plays an important role in cancer malignant biological behavior (57)(58)(59).

CONCLUSION
In this study, we identified two clusters based on 19 m 6 A-related lncRNAs and constructed a prognostic signature comprising 14 m 6 A-related lncRNAs in STAD, which had significant value in predicting the OS of patients with STAD, clinicopathological characteristics, TME, ICGs expression, and the response to ICIs. Additionally, biological processes and pathways associated with m 6 A-related lncRNAs were identified, which improved our understanding of the role of m 6 A-related lncRNAs in the occurrence and progression of STAD. This work also provides important evidence for the development of predictive biomarkers and immunotherapy for STAD.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study, which can be found in the Cancer Genome Atlas (TCGA) database.

ETHICS STATEMENT
This study met the publication guidelines stated by TCGA (https:// cancergenome.nih.Gov/publications/publicationguidelines). All data used in the study were obtained from TCGA. Ethics approval and informed consent were not required.

AUTHOR CONTRIBUTIONS
YH, PL, and FL conceived and designed the study. YH, ZY, CH, XJ, and YY organized the database and performed statistical analyses. YH and YW wrote the first draft of the manuscript. ZY, CH, and XJ prepared the figures and tables and were involved in manuscript writing. KZ, YW, PL, and FL revised and proofread the manuscript. All authors contributed to the article and approved the submitted version.