Evaluated the Role of N6-Methyladenosine-associated lncRNAs in Prognosis and Immune Microenvironment of Colorectal Cancer

Background: The role of N6-methyladenosine medicated long non-coding RNAs (lncRNAs) is elusive in colorectal cancer (CRC). Materials and Methods: We identi�ed m6A-associated lncRNAs by using the data gathered from The Cancer Genome Atlas (TCGA) and strati�ed CRC patients into different subgroups. Cox-regression analysis were performed to construct a m6A-associated lncRNA signature. And the role of this signature in prognosis and immune microenvironment was dissected subsequently. Finally, a gene set enrichment analysis (GSEA) were executed to predict the possible bio-mechanisms on the basis of the signature. Results: Three m6A-associated clusters were constructed from 866 differentially expressed lncRNAs. Cluster 2 had poor prognosis and low immune cell in�ltration. A m6A-associated lncRNA signature consisting of 14 lncRNAs was constructed, and recognized as an independent prognostic indicator in CRC by using survival analysis and receiver operating characteristic (ROC) curves. The clinical features and immune cell in�ltration status were signi�cantly different in the patients strati�ed by risk score. Furthermore, GSEA showed that 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 and immunogenicity indicators in CRC.


Introduction
Colorectal cancer is the third prevalent gastrointestinal malignancies worldwide 1 . A signi cant number of CRC patients will ultimately relapse after curative treatments [2][3] . Hence, there is an urgent requirement to investigate alternative prognostic markers in CRC patients.
Recent studies have highlighted that m6A RNA modi cation plays important roles in many biological processes including cancer pathogenesis 5 . Aberrant expressions of m6A regulators (e.g. METTL14, METTL3, KIA, ALKBH5, FTO and YTHDF1/2/3) have been identi ed in numerous tumors 6-10 . Varieties of biological functions, ranging from tumor initiation, invasion, metastasis to tumor stem cell pluripotency, could be mediated by m6A methylation [11][12] . Long non-coding RNAs (lncRNAs) are important epigenetic regulators, which play critical roles in diverse physiological and pathological processes [13][14] . Studies have reported that many lncRNAs participate in tumor initiation and progression [15][16][17] . Despite extensive efforts to de ne the pathogenesis of lncRNAs, the roles of lncRNAs in the m6A modi cation remains largely elusive in CRC. Immune microenvironment has been found to be closely implicated in the clinical outcome of immunotherapy and tumor development. In the present study, the co-expression network of the m6A-associated lncRNAs to obtain 68 m6A-associated prognostic lncRNAs was generated. Then we established three m6A-associated clusters in CRC, analyzed the characteristics of immune cell in ltration among tumor cells and investigated whether m6A-associated lncRNAs clusters have prognosis values in CRC patients. Furthermore, we constructed 14 m6A-associated lncRNAs signature which could predict the prognosis of CRC patients.

Results
The differential expressions of m6A-associated lncRNAs.
Identi cation of m6A-associated lncRNAs with prognostic value As shown in Figure1A, we annotated m6A-associated lncRNAs and clinical characteristics, then investigated the role of each lncRNA on the prognostic outcome of the patients with CRC. And there are 68 m6A-associated lncRNAs with obvious prognostic value were detected for further study.

Establishment of m6A-associated lncRNA clusters
To classify the different tumor clusters with m6A phenotypes on the basis of lncRNAs, we mapped 68 m6A-associated lncRNAs to expression pro le of CRC samples to perform consistent clustering with Consensus Cluster Plus (CCP) tool. As shown in Figure1B, the number of clusters was sequential set from 1 to 9 and CCP analysis indicated that the results were most stable when these m6A-related lncRNAs were separate into three tumor clusters ( Figure 1C, D). The OS data of each cluster was calculated using Kaplan-Meier method, and the results displayed that there were signi cant difference among the survival of the three clusters ( Figure 1E).
Clinical characteristics and Immune Score of each cluster in CRC As compared with cluster1 and cluster3, cluster 2 had the highest N stage, M stage and TNM stage ( Figure 2A). ESTIMATE-algorithm was employed to accurate estimate score (tumor purity), immune score and stromal score in accordance with the gene expression pro les of CRC patients. Our ndings showed that compared with cluser 1 and 3, the cluser 2 had the lowest estimate score, immune score and stromal score ( Figure 2B). m6A-associated lncRNAs signature construction As shown in Figure 3A, a total of 14 m6A-associated lncRNAs that had a co-expression relationship with 8 m6A-associated genes were recognized as effective independent prognostic factors. Among them, Our ndings showed that the morality occurrence was closely associated with the risk score.
Moreover, the area under the curve (AUC) is measured and the value for prognostic risk score was 0.764 which is higher than AUCs of the other clinicopathological factors ( Figure 4E). And the AUC values corresponding to 1, 3, and 5 years of OS were 0.764, 0.743 and 0.753, respectively ( Figure 3F). These data accomplished good prediction accuracy of this model.

The validation of the signature in CRC
The prognostic value of the m6A-associated lncRNA signature was investigated in CRC patients classi ed by various clinical parameters, consisting of gender, age, T, N, M and TNM stage. In almost all subgroups, the patents with low risk score trended to have higher OS rate than high-risk group (Figure 4).
Next we evaluated the independence and effectiveness of this model in predicting prognosis of CRC patients. Our ndings showed that this m6A-associated lncRNA signature could be an effective and independent factor for predicting the outcome of the patients with CRC ( Figure 5 A-B). Then, a nomogram was conducted to predict 1, 3 and 5-years OS of the patients with CRC based on the results of univariate and mutivariate Cox-regression analyses, including age, TNM stage and risk score ( Figure 5C). The calibration curves proved well prediction accuracy of this nomogram in CRC patients ( Figure 5 D-F).

Gene set enrichment analysis
Finally, we tried to detect 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 Nature Killer cells mediated cytotoxicity (NOM p-val=0.0172, FDR q-val=0.195) were more enriched in low risk group. Our study suggested that the risk-related model might contribute to organize the personalized treatment for CRC patients.

Discussion
Previous advances have demonstrated the pivotal roles of m6A modi cation in various cancers including CRC 18-20 . Investigating the potential prognostic role of m6A-associated lncRNAs will facilitate understanding the molecular mechanisms of CRC. In our work, 68 prognostic m6A-associated lncRNAs were obtained, then three m6A-associated lncRNAs cluster groups were constructed using 426 CRC samples from TCGA database. Compared with cluster1 and cluster3, cluster2 had worst OS time and later pathological stage. In addition, ESTIMATE analyses further revealed that the immune score was remarkably reduced in Cluster 2. Our data suggested that m6A-associated lncRNAs might serve as a predictive biomarker.
It is generally known that there are currently some CRC prognostic indicators, including TNM stage and tumor grade. However, More accurate prognostic factors are required to predict and analyze the OS rate in

Conclusion
In summary, our work de ned an innovative m6A-associated lncRNA signature which could provide a new perspective on predicting the prognosis of CRC patients. Furthermore, our m6A-associated lncRNA signature provides an important clue for the development of individual treatment.

Data processing of the CRC dataset
The public RNA sequencing (RNA-seq) were downloaded from TCGA. Patients without survival information were removed.
Identi cation of m6A-associated lncRNAs in CRC The m6A-associated genes were gathered from TCGA database and selected based on previously published articles. And the m6A-associated lncRNAs were screened by spearman correlation coe cient formula with /R/_value >0.6 and p_value <0.001.

Consensus clustering of m6A-associated lncRNAs
On the basis of m6A-associated lncRNAs expression levels, the CRC patients were separately divided into three groups according to the optimal k-means clustering. Cluster analysis was performed with Consensus Cluster Plus R package. The OS data of each cluster was calculated using Kaplan-Meier method. The correlation between m6A-associated lncRNAs and clinical characteristics was analyzed according to the TCGA database. And ESTIMATE-algorithm was employed to estimate the tumor immune microenvironment.
m6A-associated lncRNA signature construction The prognostic m6A-associated lncRNAs were indentied via univariate cox regression analysis. And the prognostic signature was established via mutivariate cox regression analysis. The risk scores of CRC patients were calculated by following formula, Risk score = ∑Expi*βi, here Expi represents the expression and βi represents the coe cient of m6A-associated lncRNAs. And the accuracy of the m6A-associated lncRNAs was assessed via the ROC curve analysis.

Statistical analysis
All data were analyzed via by using R statistical soft-ware version 4.0.3. A p_value less than 0.05 was statistically signi cant.

Consent for publication
Not applicable.

Availability of data and materials
The data that support our ndings are openly available in TCGA (https://portal.gdc.cancer.gov/) repository.

Competing interests
No author has con ict of interest.       The independence and effectiveness of this model in predicting prognosis of CRC patients. Forest plots of univariate (A) and mutivariate (B) Cox regression analysis in CRC. Nomogram model (C) to predict 1, 3 and 5-years survival rates of CRC patients. Calibration graph showed that predicted 1 (D), 3(E) and 5years (F) survival rates were close to actual survival rates.