The Role of N6-Methyladenosine-Associated lncRNAs in the Immune Microenvironment and Prognosis of Colorectal Cancer

Background The role of N6-methyladenosine long noncoding RNAs (lncRNAs) in colorectal cancer (CRC) is elusive. Materials and Methods We identified m6A-associated lncRNAs by using the data gathered from The Cancer Genome Atlas (TCGA) and stratified CRC patients into different subgroups. Cox regression analysis was performed to construct an m6A-associated lncRNA signature. The role of this signature in the immune microenvironment and prognosis was dissected subsequently. Finally, a gene set enrichment analysis (GSEA) was conducted to predict the possible mechanisms based on the signature. Results Three m6A-associated clusters were constructed from 866 differentially expressed lncRNAs. Cluster 2 had poor prognosis and low immune cell infiltration. An m6A-associated lncRNA signature consisting of 14 lncRNAs was constructed and recognized as an independent prognostic indicator of CRC by using survival analysis and receiver operating characteristic (ROC) curves. The clinical features and immune cell infiltration status were significantly different in patients stratified by the risk score. Furthermore, GSEA showed that the P53 pathway and natural killer cell-mediated cytotoxicity were more enriched in the low-risk group. Conclusion Our data revealed that m6A-associated lncRNAs could be potential prognostic indicators of immunogenicity in CRC.


Introduction
Colorectal cancer (CRC) is the third most prevalent gastrointestinal malignancy worldwide [1]. A significant number of CRC patients will ultimately relapse after curative treatments [2]. Hence, there is an urgent need to investigate prognostic markers for CRC.
Immune microenvironment has been found to be closely associated with the clinical outcome of immunotherapy and tumor development. In the present study, the coexpression network of the m6A-associated lncRNAs was investigated to obtain 68 m6A-associated prognostic lncRNAs. en, we established three m6A-associated clusters in CRC, analyzed the characteristics of immune cell infiltration among tumor cells, and investigated whether m6A-associated lncRNAs clusters have prognosis values in CRC patients. Furthermore, we constructed a signature using 14 m6A-associated lncRNAs which could predict the prognosis of CRC patients.

2.2.
Identification of m6A-Associated lncRNAs with a Prognostic Value. As shown in Figure 1(a), we annotated m6Aassociated lncRNAs and clinical characteristics, then investigated the role of each lncRNA on the prognostic outcome of the patients with CRC. A total of 68 m6A-associated lncRNAs with obvious prognostic values were detected and used for further study.

Establishment of m6A-Associated lncRNA Clusters.
To classify different m6A clusters based on lncRNAs, we mapped these 68 m6A-associated lncRNAs to expression profile of CRC samples to perform clustering using the Consensus Cluster Plus (CCP) tool. As shown in Figure 1(b), the number of clusters was sequentially set from 1 to 9, and CCP analysis indicated that the results were most stable when these m6A-related lncRNAs were separated into three clusters using the Consensus Cluster Plus R package (Figures 1(c), 1(d)). e OS data of each cluster was calculated using the Kaplan-Meier method, and the results displayed that there was significant difference among the survival of CRC patients in these three clusters (Figure 1(e)).

Clinical
Characteristics and the Immune Score of Each Cluster in CRC. As compared with cluster 1 and cluster 3, cluster 2 had the highest N stage, M stage, and TNM stage (Figure 2(a)). e ESTIMATE algorithm was employed to evaluate the accurate estimate score (tumor purity), immune score, and stromal score in accordance with the gene expression profiles of CRC patients. Our findings showed that compared with clusters 1 and 3, cluster 2 had the lowest estimate score, immune score, and stromal score (Figure 2(b)).
2.5. m6A-Associated lncRNAs Signature Construction. As shown in Figure 3(a), a total of 14 m6A-associated lncRNAs that had a coexpression relationship with 8 m6Aassociated genes were recognized as effective independent prognostic factors. Among them, AC137932. 3 (Figures 3(b)-3(d)). Our findings showed that the mortality was closely associated with the risk score. Moreover, the area under the curve (AUC) is measured, and the value for the prognostic risk score was 0.764 which is higher than AUCs of the other clinicopathological factors (Figure 4(e)). e AUC values corresponding to 1-, 3-, and 5-year of OS were 0.764, 0.743, and 0.753, respectively ( Figure 3(f )). ese data indicated the good prediction accuracy of this model.

e Validation of the Signature in CRC.
e prognostic value of the m6A-associated lncRNA signature was investigated in CRC patients from TCGA dataset. e patients were classified by various clinical parameters, consisting of gender, age, T, N, M, and TNM stage. In almost all subgroups, the patients with a low-risk score trended to have a higher OS rate than that of the high-risk group ( Figure 4).
Next, we evaluated the independence and effectiveness of this model in predicting prognosis of CRC patients. Our findings showed that this m6A-associated lncRNA signature could be an effective and independent factor for predicting the outcome of the patients with CRC ( Figures 5(a), 5(b)). en, a nomogram was conducted to predict 1-, 3-and 5year OS of the patients with CRC based on the results of univariate and multivariate Cox regression analyses, including age, TNM stage, and the risk score ( Figure 5(c)). e calibration curves demonstrated well-prediction accuracy of this nomogram in CRC patients (Figures 5(d)-5(f )).

Gene Set Enrichment
Analysis. Finally, we evaluated the potential biological mechanisms associated with the risk model by GSEA. As shown in Figure 6, the P53 signaling pathway (NOM p-val � 0.0019, FDR q-val � 0.155) and natural killer cell-mediated cytotoxicity (NOM p-val � 0.0172, FDR q-val � 0.195) were more enriched in the low-risk group. Our study suggested that this risk-related model could be used for the personalized treatment of CRC patients.

Discussion
Previous studies have demonstrated the pivotal roles of m6A modification in various cancers including CRC [14][15][16]. Investigating the potential prognostic role of m6A-associated lncRNAs will facilitate understanding of the molecular mechanisms of CRC. In our work, 68 prognostic m6A-associated lncRNAs were identified, then three m6A-  associated lncRNAs cluster groups were constructed using 426 CRC samples from TCGA database. Compared with cluster 1 and cluster 3, cluster 2 had the worst OS time and high pathological stage. In addition, ESTIMATE analyses revealed that the immune score was remarkably reduced in cluster 2. Our data suggested that m6A-associated lncRNAs might be used as a predictive biomarker. It is generally known that there are currently some CRC prognostic indicators, including the TNM stage and tumor grade. However, more accurate prognostic factors are           TNM stage and identified that the signature was closely associated with the progression of CRC. Meanwhile, the GSEA analysis preliminary displayed that these lncRNAs were closely involved in the P53 pathway and NK cellmediated cytotoxicity. Further studies are needed to demonstrate the mechanisms involved in this lncRNA signature.

Conclusion
In summary, our work defined a m6A-associated lncRNA signature which could predict the prognosis of CRC patients.
is m6A-associated lncRNA signature will provide guidance for individualized treatment.

Data Acquisition and Processing of the CRC Dataset.
e public RNA sequencing (RNA-seq) data from 512 patients' CRC were downloaded from TCGA (https://portal. gdc.cancer.gov/). Patients without survival information were removed.

Identification of m6A-Associated lncRNAs in CRC.
e m6A-associated genes were gathered from TCGA database and selected based on the previously published articles [20,21]. e m6A-associated lncRNAs were screened by Spearman correlation coefficient formula with|R| value > 0.6 and p value <0.001.

Consensus Clustering of m6A-Associated lncRNAs.
On the basis of the expression levels of m6A-associated lncRNAs, the CRC patients were separately divided into three groups (clusters 1, 2, and 3) according to optimal kmeans clustering. Cluster analysis was performed with the Consensus Cluster Plus R package. e overall survival (OS) data of each cluster was calculated using the Kaplan-Meier method. e correlation between m6A-associated lncRNAs and clinical characteristics was analyzed according to TCGA database. e ESTIMATE algorithm was employed to estimate the tumor immune microenvironment.

m6A-Associated lncRNA Signature Construction.
e prognostic m6A-associated lncRNAs were identified via univariate cox regression analysis. e prognostic signature was established via multivariate cox regression analysis. e risk scores of CRC patients were calculated by the following formula: Risk score � Expi * βi, where Expi represents the expression, and βi represents the coefficient of m6A-associated lncRNAs. e accuracy of the m6Aassociated lncRNAs was assessed via the ROC curve analysis.

Statistical Analysis.
All data were analyzed via by using R statistical software version 4.0.3. A p value less than 0.05 was statistically significant.    14 Journal of Oncology Data Availability e data that support our findings are openly available in TCGA (https://portal.gdc.cancer.gov/) repository.