Constructing a novel signature and predicting the immune landscape of colon cancer using N6-methylandenosine-related lncRNAs

Background: Colon cancer (CC) is a prevalent malignant tumor that affects people all around the world. In this study, N6-methylandenosine-related long non-coding RNAs (m6A-related lncRNAs) in 473 colon cancers and 41 adjacent tissues of CC patients from The Cancer Genome Atlas (TCGA) were investigated. Method: The Pearson correlation analysis was conducted to examine the m6A-related lncRNAs, and the univariate Cox regression analysis was performed to screen 38 prognostic m6A-related lncRNAs. The least absolute shrinkage and selection operator (LASSO) regression analysis were carried out on 38 prognostic lncRNAs to develop a 14 m6A-related lncRNAs prognostic signature (m6A-LPS) in CC. The availability of the m6A-LPS was evaluated using the Kaplan–Meier and Receiver Operating Characteristic (ROC) curves. Results: Three m6A modification patterns with significantly different N stages, survival time, and immune landscapes were identified. It has been discovered that the m6A-LPS, which is based on 14 m6A-related lncRNAs (TNFRSF10A-AS1, AC245041.1, AL513550.1, UTAT33, SNHG26, AC092944.1, ITGB1-DT, AL138921.1, AC099850.3, NCBP2-AS1, AL137782.1, AC073896.3, AP006621.2, AC147651.1), may represent a new, promising biomarker with great potential. It was re-evaluated in terms of survival rate, clinical features, tumor infiltration immune cells, biomarkers related to Immune Checkpoint Inhibitors (ICIs), and chemotherapeutic drug efficacy. The m6A-LPS has been revealed to be a novel potential and promising predictor for evaluating the prognosis of CC patients. Conclusion: This study revealed that the risk signature is a promising predictive indicator that may provide more accurate clinical applications in CC therapeutics and enable effective therapy strategies for clinicians.


Introduction
Globally, colon cancer is a common malignant tumor in humans that has a high morbidity and mortality rate (Lu and Zhang, 2020;Xiaoyong et al., 2022). A poor prognosis leads to CC posing a serious threat to the health of humans. Most colon cancers are caused by poor dietary habits, age, lifestyle, lack of exercise, and smoking (Xavier and Podolsky, 2007). In recent years, there have been increased studies focusing on the modification of m6A within the epigenetics field. The m6A modification, which transfers the methyl group to the nitrogen-6 position of the adenosine base in RNA, is the most abundant and reversible mRNA epigenetic modification (Zhang et al., 2021a). The m6A methyltransferase, also known as "writers", installs m6A modification from RNA, while the demethylase, also known as "erasers", removes m6A modification from RNA (Xu et al., 2021a). After m6A modification, the mature mRNA is recognized by the "reader" when it is exported from the nucleus to the cytoplasm (Zaccara et al., 2019). The "reader" is a binding protein capable of recognizing important chemical signals for m6A modification. As a promising biomarker, m6A-related regulatory factors participate in various biological processes in the occurrence and progression of numerous diseases, especially malignant tumors (Tian et al., 2020). The m6A-related regulatory factors are involved in almost every RNA metabolism and processing step that influence RNA function (Lan et al., 2021).
Long non-coding RNA (lncRNA), a type of transcribed RNA, has various important biological functions (Johnsson et al., 2014). As lncRNA is vital in almost all aspects of biological function, it is indispensable for modulation, especially m6A. Furthermore, lncRNA may also regulate m6A methylation through certain pathways. In addition, m6A modification may affect lncRNA by several regulatory mechanisms. Furthermore, m6A modification may alter the structure of sectional RNA, allowing the corresponding RNA-binding proteins to enter the ambient m6A residues. Moreover, modification of m6A may affect the formation of the RNA-DNA triplex, modulating the binding of RNA to target DNA, and thereby regulating target genomic sites (Ma et al., 2019). It has been reported that m6A-related lncRNAs contribute to improving prognostic risk assessment and the development of individualized therapy decisions in various cancers. For example, Feng et al. (2021) demonstrated the prognostic predictive value of the 6A-LPS, which included four m6A-related lncRNAs, and evaluated the correlation with PD-L1 in head and neck squamous cell carcinoma. Xu et al. (2021b) elucidated a risk signature composed of 12 m6A-related lncRNAs as promising prediction targets in lung adenocarcinoma prognosis and immune responses. The m6A-related lncRNA signature has been shown to have significant molecular prognostic value in gastric cancer and may improve the development of personalized immunotherapy strategies . However, the specific mechanism through which m6A-related lncRNAs promote CC occurrence and development remains unknown.
In this study, the expression data of 14,142 lncRNAs and 21 m6A regulators from TCGA were analyzed to identify three m6A modification patterns, as well as develop an m6A-LPS using the LASSO regression analysis. The 14 lncRNAs used to construct the m6A-LPS were reported in CC for the first time. In addition, a nomogram was established to predict the survival time of CC patients. Finally, the relationship with immunotherapy responses was investigated.

Materials and methods
Datasets and screening of m6A-related lncRNAs The transcriptome sequencing results, relevant clinical details, and somatic mutation data of CC patients were collated via TCGA (https:// cancergenome.nih.gov/), which embodied mRNA and lncRNA expression data from 473 colon cancers and 41 adjacent tissues of CC patients. The GTF file was downloaded via Ensembl for annotation to distinguish between mRNA and lncRNAs. Furthermore, a list of 21 m6Arelated genes was screened based on previous literature, including 8 writers (RBM15, RBM15B, WTAP, METTL3, METTL14, KIAA1429, CBLL1, and ZC3H13), 2 erasers (FTO and ALKBH5), and 11 readers (YTHDC1, YTHDC2, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPC, HNRNPA2B1, YTHDF1, YTHDF2, YTHDF3, and RBMX). If the expression level of a lncRNA was related to one or more of these genes, it was identified as an m6A-related lncRNA. The correlation between lncRNA and m6a regulators in the dataset was determined by the Pearson method. The square of correlation coefficient|R 2 |>0.5 and p < 0. 001 were used as the identification criteria for m6A-related lncRNAs.

Consensus clustering and comparison of immune cell infiltration
Consensus clustering analysis was used to identify distinct modification patterns and cluster the TCGA-COAD samples based on the expression levels of prognostic lncRNAs. Consensus clustering was used to determine the optimal number of clusters. The ConsensusClusterPlus R package, which was the unsupervised clustering method, provided stable visual evidence to estimate the number of unsupervised clusters in a dataset (Wilkerson et al., 2010). In order to ensure the classification was stable, 1,000 repetitions were performed utilizing ConsensuClusterPlus package (Wilkerson et al., 2010). In order to analyze the differences in immune infiltration levels among multiple clusters, 22 immune cell types infiltrating the TCGA-COAD samples were identified using the CIBERSORT package Y et al., 2022). The p-value corresponding to the output result is of < 0.05 was considered statistically significant.

Construction of the m6A-LPS
The single-factor Cox regression was used in conjunction with survival data, and the prognosis-related lncRNAs with p < 0.05 were screened using the log-rank test according to the TCGA database. Subsequently, the Lasso regression analysis was carried out on these lncRNAs to develop the m6A-LPS. The risk score for each case was determined based on the expression of the predicted lncRNA multiplied by the coefficient from the LASSO algorithm. The Kaplan-Meier (K-M) method was used for assessing the difference Frontiers in Genetics frontiersin.org in overall survival (OS). The receiver operating characteristic (ROC) curve was drawn using the R package "survivalROC". The sensitivity and specificity of the m6A-LPS were examined using the area under the curve (AUC). In order to determine if the m6A-LPS and clinicopathological features served as independent factors for OS, the univariate and multivariate Cox regression analyses were carried out. The principal component analysis (PCA) was used to compare the distribution of high-and low-risk score patients in CC by the "prcomp" function of the "stats" R package.
Estimation of immune landscape and immunosuppressive molecules with the m6A-LPS In order to assess the various immune landscapes in the different risk subgroups, the currently acknowledged methods were utilized to identify the infiltration of 22 immune cell types in the TCGA-COAD samples, including TIMER (Li et al., 2017;Wang et al., 2020), CIBER SORT (Chen et al., 2018;Zhang et al., 2020), XCELL (Aran et al., 2017;Aran, 2020), QUANTISEQ (Plattner et al., 2020), MCPcounter (Dienstmann et al., 2019), EPIC (Racle et al., 2017), and CIBERSORT-ABS (Tamminga et al., 2020). A detailed Spearman correlation analysis was performed using the R ggplot2 packages. A p-value of < 0.05 was considered statistically significant. The differences of 22 immune cell types were also analyzed by the Wilcoxon test using the R ggpubr package. The differences in the expression levels of biomarkers related to Immune Checkpoint Inhibitors (ICIs) between different subgroups were examined using a violin plot obtained with the R ggpubr package.

GSEA enrichment analysis
We used GSEA to examine immune function and biological pathways. Functional enrichment analyses via Gene Ontology (GO) and the Kyoto Gene and Genomic Encyclopedia (KEGG) pathway analyses, were conducted by gene set enrichment analysis (GSEA) 4.1.0 (Subramanian et al., 2005;Kanehisa et al., 2017). p < 0.05 was considered a significant category.

Quantitative real-time PCR
In accordance with the instructions of the manufacturer, cDNA synthesis was performed by using the SYBR-Green Mix (Vazyme Biotech Co.,Ltd.). The Real-time PCR (qRT-PCR) analysis was conducted on The LightCycler 480 Real-Time PCR System. The primers used for qRT-PCR were purchased from Servicebio (Wuhan China). The related GAPDH mRNA expression was used as an endogenous control. The primer sequences used in our study were summarized in Table 1.

Statistical analysis
Gene expression and risk scores between subgroups were compared using the Wilcox test. The Kaplan-Meier (K-M) method was used for assessing the difference in OS among three subgroups. An assessment of linear correlation between two random variables was performed using the Pearson correlation coefficient. All statistical analyses were performed using R (version 3.6.2), GraphPad Prism 9 and Perl software in the present experiment. All analyses performed were bilateral, with p < 0.05 being statistically significant.

Determining m6A-related lncRNAs in CC patients
The workflow for the m6A-LPS analysis is illustrated in Figure 1. Initially, the expression data of 14,142 lncRNAs and 21 m6A regulators from TCGA, as well as the survival status, and clinical

Gene
The primer sequences (5′-3′) Frontiers in Genetics frontiersin.org data, were identified and obtained. A total of 1,234 lncRNAs were identified by determining the Pearson coefficient of the lncRNA-mRNA co-expression analysis ( Figure 2A; Supplementary Table S1). As shown in Figure 2B; Supplementary Table S2, 38 prognostic m6A-related lncRNAs were filtered through a univariate Cox regression analysis. The differences in expression of the 38 prognostic lncRNAs in CC and normal samples were analyzed using the Wilcoxon test. The differential expression of 38 lncRNA heatmaps ( Figure 2C) and boxplots ( Figure 2D) were plotted using the "pheatmap" and "ggpot2" packages, respectively. The correlations between the 38 prognostic lncRNAs and PD-L1 in TCGA were examined using Spearman's method, as displayed in Figure 2E Consensus clustering and comparison of immune cell infiltration Consensus cluster analysis was used to further explore the expression characteristics of m6A-related lncRNAs in the TCGA-COAD samples. The TCGA-COAD cohort was clustered into k subtypes (k = 2-9) based on the 38 prognostic lncRNA expression levels with the R package ConsensusClusterPlus (Figures 3A-D; Supplementary Figure S1). The proportion of ambiguous clustering measurements and the similarity in expression levels we measured for the m6A-related lncRNAs shared by TCGA ultimately determined k = 3 to have the best cluster stability. There was also the least crossover between CC samples, when the consistency matrix with a k value of 3 was selected.
In order to explore the clinical application value among the three clusters, the relationship between cluster and clinical features was evaluated using the Chi-square test, and the obtained heatmap is displayed in Figure 3G Among the three clusters, there was a significant difference in the N stage and, according to the Kaplan-Meier method (p = 0.010) ( Figure 3E), a clear difference in survival time. The intracluster proportions for the three clusters were then analyzed based on age, gender, T, N, M, and clinical stage as shown in Figure 3H Subsequently, the CIBERSORT algorithm was used to analyze the infiltration of 22 immune cell types among three clusters. It was found that the Tregs differed significantly among the three clusters (Supplementary Figure S2). The T cells regulatory (Tregs) differed significantly among three clusters ( Figure 3F). Furthermore, cluster3 displayed a greater number of Mast cells resting compared with cluster1. The NK cells resting, CD4 T cells memory activated, CD8 T cells, T cells gamma delta, and T cells regulatory differed significantly between cluster1 and cluster2. The T cells CD4 memory activated, NK cells resting, T cells gamma delta, and T cells regulatory differed significantly between cluster3 and cluster2.

Construction of the m6A-LPS
The LASSO regression analysis was performed on 38 prognostic lncRNAs to develop a 14 m6A-LPS in TCGA dataset (Figures 4A, B;     The m6A-LPS predicts overall survival in patients with CC. LASSO regression was performed, calculating the minimum criteria (A,B). The patients in the high-risk cohort had significantly shorter OS than those in the low-risk cohort in the training (C) and test sets (F); ROC curves of the m6A-LPS and clinical features for prediction of 3-year OS in the training (D) and test sets (E); risk score distribution, the distribution of survival time and survival status, and the heatmap of the expression of 14 m6A-related prognostic lncRNAs in the training (G) and test sets (H).
Frontiers in Genetics frontiersin.org Supplementary Table 3). The correlations between the 14 prognostic lncRNAs and the m6A genes in TCGA were displayed in Figure 2F As shown in Figure 2G; Table 2, forest maps corresponding to these 14 prognostic lncRNAs were drawn, with the prognostic favorable factors TNFRSF10A-AS1, AC099850.3, AL137782.1, and AC073896.3 among them. The rest, on the other hand, were classified as prognostic unfavorable factors, and their overexpression might reduce survival. The Kaplan-Meier method was used to determine the prognostic value of the 14 hub lncRNAs in CC, with the median expression being the cutoff value. As displayed in Supplementary Figure ). In order to validate the accuracy of the m6A-LPS, all cases were randomly separated into either a training set or test set. Furthermore, all cases in the TCGA-COAD cohort were separated into high-and low-risk subgroups in every set based on the median risk score. The K-M analysis revealed that the OS of high-risk patients declined compared to the low-risk patients in the training set (p < 0.001) ( Figure 4C), which was further validated by the results of the test set (p = 0.001) ( Figure 4F). Finally, the ROC curve revealed that the m6A-LPS exhibited a relatively decent predictive value in both the training and test sets in predicting a 3-year OS of CC in comparison to other clinical factors, as shown in Figures 4D, E. Figures 4G, H displays the risk score distribution, survival status, survival time, and expression heatmap based on the m6A-LPS in both the training and test sets. The ROC analysis revealed that the AUCs of the m6A-LPS at 1-, 3-, 5-year was 0.714, 0.777, 0.799 in the training set, respectively ( Figure 5A); and that the AUCs of the m6A-LPS at 1-, 3-, 5-year was 0.686, 0.768, 0.774 in the test set ( Figure 5B), respectively. The calibration curves for the probability of OS at 1, 3, 5 year showed good consistency between the actual observation and the m6A-LPS prediction (Figures 5C, D). According to Figures 5E, F, the mortality rate of the high-risk patients was higher in the training and test sets than in the low-risk patients. Following that, PCA indicated that the training ( Figure 5G) and test sets ( Figure 5H) should be distributed in two directions.

An independent prognostic analysis and construction of an m6A-LPS-based nomogram
The univariate and multivariate Cox regression analyses were conducted to confirm whether the m6A-LPS and clinicopathological features had independent prognostic characteristics for patients with CC in the TCGA-COAD cohort. The univariate Cox regression analyses revealed that the m6A-LPS, age, clinical stage, T, M, and N were all significantly related to OS ( Figure 5I). Furthermore, the findings of the multivariate Cox analysis indicated that both m6A-LPS, and T were independent predictors for determining the OS of CC patients in the TCGA-COAD cohort ( Figure 5J). As an applicable quantitative tool in the clinic, a nomogram containing m6A-LPS with clinical features was established to evaluate the life expectancy of CC patients in the TCGA-COAD cohort ( Figure 6A). The total score for each case was obtained by calculating the points of every parameter to predict the 1-year, 3-year, and 5-year OS among CC patients. The calibration curves for the probability of OS at 1, 3, 5 year showed good consistency between the actual observation and the nomogram prediction ( Figure 6B). Finally, as shown in Figure 6C, the ROC curve revealed that the m6A-LPS-based nomogram had relatively decent predictive value in predicting 1-year (AUCs = 0.779), 3-year (AUCs = 0.812) and 5-year (AUCs = 0.820) OS.

PCA confirms the ability of the m6A-LPS to group
To validate risk models, PCA was also applied, and the results were displayed using the R software's "scatterplot3D" tools. PCA was used to identify the distinct patterns of m6A distribution on expression profiles of the whole gene ( Figure 6E), 21 m6A genes ( Figure 6F), m6A-related lncRNAs( Figure 6G), and 14 hub lncRNAs of the m6A-LPS ( Figure 6H). It was clear from the m6A-LPS results that there were more noticeable differences in distributions between the low-and high-risk groups than from the other three approaches. Based on the m6A-LPS, we intuitively observed that CC patients were effectively divided into two subgroups. The risk score's concordance index increased over time, outpacing the concordance index of other clinical indicators ( Figure 6D). In light of these findings, the m6A-LPS may be a better predictor of the prognosis of CC.

Clinical correlation analysis
In order to explore the clinical practice of the m6A-LPS, the chisquare test was used to investigate the relationship between the m6A-LPS and clinical features. The heatmap from the resulting diagram is displayed in Figure 8A. The cluster, T, N, and clinical stages were all significantly different in the various risk-subgroups.  Frontiers in Genetics frontiersin.org AC245041.1 was significantly linked to age; AL137782.1 was significantly associated linked to the M stage.

Estimation of the immune landscape and immunosuppressive molecules with the m6A-LPS
A detailed correlation analysis using the Spearman method was performed to determine the correlation between the m6A-LPS and tumor-infiltrating immune cells based on the cohort data from TCGA-COAD, with the results shown in a bubble plot ( Figure 8A). The results indicated that the m6A-LPS had a negative relation with neutrophils, plasmacytoid dendritic cells, and resting mast cells, and a positive relation with regulatory T cells, monocytes, B cells, and CD8 + T cells (Supplementary Figure 5; Supplementary Table 5). The results were similar to consensus clustering analysis. The regulatory T cells are highly expressed in both cluster 2 and high-risk groups, and both have the worse prognosis ( Figures 7N, O). As ICIs were recommended for treating CC in clinics, the differences in ICI-related biomarkers between the two

Investigation of the correlation between the m6A-LPS and chemotherapeutics
In order to evaluate the m6A-LPS in clinical practice for treating CC, the therapeutic response was estimated by calculating the half-maximal inhibitory concentration (IC 50 ) of the standard chemotherapeutic drugs for each sample. IC 50 of chemotherapeutics was obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) website in TCGA. The differences in the IC 50 between the high-and low-risk subgroups were examined using the Wilcoxon signed-rank test. The results were displayed in a box plot obtained utilizing the R packages ggpubr, pRRophetic, and ggplot2. The box plot revealed a positive association between a highrisk score and a higher IC 50 of mitomycin (p = 0.0091) ( Figure 8G) and gemcitabine (p = 0.0038) ( Figure 8H), suggesting that the m6A-LPS has potential as a predictor in CC patients for chemosensitivity. , and TBX2 (F); *p < 0.05, ** p < 0.01, and *** p < 0.001. Boxplots evaluating the response to the mitomycin (G) and gemcitabine (H) chemotherapeutic between high-and low-risk patients.
Frontiers in Genetics frontiersin.org

Exploration of the correlation between the m6A-LPS and somatic variants
The maftools R package was utilized to explore the differences between the two subgroups in tumor mutation burden (TMB) in the TCGA-COAD cohort. Figures 9A, B showed the first 20 most frequently mutated genes in the high-and low-risk subgroups. It indicates that the low-risk subgroup has a more extensive TMB than the high-risk subgroup. Then, all cases in the TCGA-COAD cohort were separated into high TMB and low TMB subgroups according to the cutoff of TMB. The K-M analysis revealed that the OS of the high TMB subgroup declined compared to the low TMB subgroup (p = 0.036) ( Figure 9C). Following that, the synergistic effect of TMB and risk score in the prognostic stratification of CC patients was assessed. The survival curve of TMB combined with risk score showed that there was a clear difference in survival time among the four subgroups (p < 0.001), including high TMB & high-risk score, high TMB & low-risk score, low TMB & low-risk score, and Low TMB & high-risk score ( Figure 9D).

GSEA enrichment analysis
GSEA was applied to predict the potential GO and KEGG pathways, and the top 5 most relevant items were selected. Figures 10A-F showed the potential functions and access among the three clusters. Following that, GSEA was applied to predict the potential GO and KEGG pathways in the high-and low-risk subgroups. With regard to KEGG pathway analysis as shown in Figure 10G, we found that the following pathways were active in Frontiers in Genetics frontiersin.org high-risk subgroups: "oocyte meiosis", "pathways in cancer", "regulation of actin cytoskeleton"; and that the following pathways were active in low-risk: "drug metabolism cytochrome p450", "retinol metabolism", "glycine serine and threonine metabolism", "tryptophan metabolism", "fatty acid metabolism". In addition, the results of GO analysis ( Figure 10H) indicated that

FIGURE 10
Significantly enriched GO terms and KEGG pathways via GSEA. The top five most relevant items GO terms (A) and KEGG pathways (B) between the cluster2 and cluster1; The top five most relevant items GO terms (C) and KEGG pathways (D) between the cluster3 and cluster1; The top five most relevant items GO terms (E) and KEGG pathways (F) between the cluster3 and cluster2; The top five most relevant items GO terms (G) and KEGG pathways (H) between the high-and low-risk subgroups.
Frontiers in Genetics frontiersin.org 14 "negative regulation of cellular amide metabolic process," and "intrinsic apoptotic signaling pathway," were mainly enriched in high-risk subgroups, and that "monocarboxylic acid catabolic process," "organic acid catabolic process," "aromatase activity," "fatty acid catabolic process," and "fatty acid beta-oxidation" were mainly enriched in low-risk subgroups. The results revealed the differentially enriched pathway and biological processes were mostly related to tumor progression and metabolism.

Validation of survival predictive ability of the m6A-LPS
The data of 452 CC patients in the TCGA-COAD cohort were examined. In addition, a K-M curve was drawn to demonstrate the correlation between clinicopathologic characteristics and prognosis. The results indicated that age, T, N, M, and clinical stage were statistically significant for prognosis (Supplementary Figure S4). Following that, a stratified analysis for CC patients was performed, which revealed that the m6A-LPS, as a dependable prognostic indicator, had a robust predictive ability for the survival of CC patients ( Figure 11A-F). However, the m6A-LPS was not able to independently estimate the survival of CC patients in the T1-2 ( Figure 11D) subgroups. Frontiers in Genetics frontiersin.org 16

Discussion
LncRNAs have been shown to promote or inhibit proteincoding genes in a number of physiological and pathological processes (Wilusz et al., 2009;Ma et al., 2017;Zhang et al., 2019). Increasing evidence suggests that m6A-related genes can regulate the incidence and progression of a variety of malignant tumors by modifying specific lncRNAs. Furthermore, lncRNAs could be used as competitive endogenous RNAs to regulate tumor invasive progression by targeting m6A-related genes.
Recently, several studies reported that the m6A-associated regulators were dysregulated in CC, as well as being involved in the initiation and progression of the disease. Ni et al. (2019) identified a novel negative functional loop, the lncRNA GAS5-YAP-YTHDF3 axis, in a clinical study. Compared to the adjacent tissues, the expression of GAS5 was lower in CRC tissues. In CRC patients, the high expression levels of YAP and YTHDF3 were linked to the low expression of GAS5. These findings offer a novel therapeutic option for CRC patients. According to , METTL14 was downregulated in CRC and the low expression of METTL14 may be associated with the proliferation and invasion in CRC. Moreover, it was discovered that lncRNA XIST was a downstream target of METTL14, and that the expression of XIST was negatively correlated with METTL14 and YTHDF2. The lncRNA RP11 is upregulated in clinical settings and is regulated by m6A methylation. Furthermore, m6A-induced lncRNA RP11 can promote CRC cell metastasis and proliferation through the upregulation of Zeb1 (Wu et al., 2019).
The majority of previously published studies focused on the intrinsic carcinogenic pathways of CC. However, how lncRNA interacts with m6A regulators in carcinogenesis and progression of CC, as well as the overall patterns of the m6A-related lncRNA in CC, remain unknown. Therefore, in-depth studies of the m6Arelated lncRNA will contribute to the identification of therapeutic biomarkers that have clinical prognostic value and the development of more effective therapeutic strategies. In the present study, we identified the m6A-related lncRNAs using correlation analysis implemented by the Pearson method from TCGA. The TCGA confirmed 38 m6A-related lncRNAs related to the prognosis, and three m6A modification patterns with significantly different N stages, survival time, and immune landscape were identified according to these lncRNAs.
Additionally, 14 of 38 prognostic lncRNAs were applied to develop an m6A-related lncRNA model to predict the OS of CC patients. In further analysis, the overall survival, cluster, T, N, and clinical stage significantly differed in the different risksubgroups in further analysis. Notably, the ROC curve revealed that the m6A-LPS exhibited a relatively decent predictive value in both the training and test sets in predicting a 3-year OS of CC in comparison to other clinical factors. In addition, a nomogram was established to predict the survival time of CC patients. Finally, the relationship with immunotherapy responses was investigated. Surprisingly, there was a significant correlation between the m6A-LPS and the tumor-infiltrating immune cells. The m6A-LPS has potential as a predictor of ICI treatment response and chemosensitivity in CC patients. Chang et al. (2021) identified a positive feedback loop "ITGB1-DT/ITGB1/Wnt/β-catenin/MYC" among the 14 hub lncRNAs, significantly promoting the proliferation, migration, and invasive ability of lung adenocarcinoma cells. The enhanced expression of ITGB1-DT, an oncogenic lncRNA, has been linked to the poor OS and disease-free survival (Chang et al., 2021). NCBP2-AS1 has been implicated in a model of colon adenocarcinoma recurrence prognosis based on competing endogenous RNAs (Jin et al., 2020). In renal cell carcinoma, AP006621.2 has been identified as a redox-related prognostic factor (Qi-Dong et al., 2020). In bladder cancer, SNHG26 was discovered to be an epithelial-mesenchymal transition-related prognostic factor (Tong et al., 2021). In addition, AC073896.3 and TNFRSF10A-AS1 have been used as components of an autophagy-related lncRNA signature to improve colorectal cancer prognosis (Wei et al., 2020;Zhou et al., 2020). Meanwhile, AC073896.3 and TNFRSF10A-AS1 are prognostic favorable factors for colon cancer, which is consistent with our study (Wei et al., 2020;Zhou et al., 2020). AL137782.1 and AC073896.3 are protective factors for colon cancer, which is consistent with our study (Zhang et al., 2021b). Recently, AC099850.3 was proven to be involved in the migration and proliferation of hepatocellular carcinoma by regulating the expression of cell-cycle-related molecules (Wu et al., 2021). In addition, AC099850.3 was highly expressed in HCC and squamous cell carcinoma of the tongue and increased the predicted value of the signature through coexpression and the ceRNA mechanism Jia et al., 2020;Wu et al., 2020;Jiang et al., 2021). However, there have been few reports on the other lncRNAs.
Compared with the similar study by Chen S , our model has a better predictive value for OS of cancer patients at 1-, 3-, 5-year. Furthermore, we also comprehensively investigated the relationship between the model and consensus cluster, tumor infiltration immune cells, immunosuppressed biomarkers, somatic variants, and chemotherapeutic drug efficacy. Inevitably, there were several limitations to this study. As a retrospective study, this study demonstrated a certain degree of heterogeneity among patients, hence, future validation studies are required. In addition, the m6A-LPS, as well as its correlation with TME, did not undergo external verification because of insufficient data from other independent cohorts. Therefore, several methods were utilized to verify this novel prognostic signature. Finally, the data of TCGA released publicly was mainly mined and analyzed as this study required further external validation in multicenter cohorts.
Frontiers in Genetics frontiersin.org

Conclusion
This is the first study that comprehensively identified and systematically analyzed the expression data of m6A-related lncRNAs in CC in the TCGA database. The m6A-LPS, which is based on 14 m6Arelated lncRNAs, has been revealed to be a novel potential and promising biomarker for evaluating the prognosis of CC patients. The survival rate, clinical features, tumor infiltration immune cells, immunosuppressed biomarkers, and chemotherapeutic drug efficacy were all re-evaluated. This study revealed that the risk signature is a promising predictive indicator that may provide more accurate clinical applications in CC therapeutics and enable effective therapy strategies for clinicians.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions
YWa and HC conceived and designed the study, and revised the manuscript. DZ, YL, YWu, HM, XL, HW, XJ, GZ, and LF conducted all data collection and analysis and compiled charts. All authors read and approved the final manuscript. All authors contributed to the article and approved the submitted version.