Transcriptome profiling data, somatic mutation data, and corresponding clinical data for the TCGA-CORD cohort were downloaded from The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov/), including data from 471 CC and 41 normal case samples. Gene expression profiles were then fully annotated with the Gencode project (Frankish et al., 2019) and distinguished into mRNAs and lncRNAs profiles. The GSE17536 dataset (N = 177) was obtained from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) as an external validation set to better verify the role of target lncRNAs.

Through the review of the latest literature, a total of 50 m6A-, m1A-, m5C-, and m7G-RNA methylation regulators were obtained. Among them, 25 m6A regulators (METTL3, METTL14, METTL16, WTAP, KIAA1429, VIRMA, RBM1, RBM15, RBM15B, and ZC3H13, FTO, ALKBH5, YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3 IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, HNRNPG, RBMX, LRPPRC, and FMR1) (Li et al., 2019b; Hu et al., 2019; An and Duan, 2022), 13 m1A regulators (TRMT6, TRMT61A, TRMT61B, TRMT61C, TRMT10C, BMT2 RRP8, YTHDF1, YTHDF2, YTHDF3, and YTHDC1, ALKBH1, and ALKBH3) (Xie et al., 2020; Song et al., 2021a), 14 m5C regulators (NOP2, NSUN1, NSUN2, NSUN3, NSUN4, NSUN5, NSUN7, DNMT1, TRDMT1, DNMT3A, DNMT3B, TET2, YBX1, and ALYREF) (Meng et al., 2021), and 2 m7G regulators (METTL1 and WDR4) (Tomikawa, 2018) were included. RNA methylation regulators differentially expressed in colon cancer and normal tissues in the TCGA-CORD cohort were identified using the “limma” package. The “maftools” package was used to generate mutation maps of RNA methylation regulators in CC patients. CNV altered positions of RNA methylation regulators on 23 chromosomes were mapped using the “RCircos” package.

Pearson correlation analysis was used to screen for lncRNAs co-expressed with differentially expressed RNA methylation-related genes (|Pearson R|>0.5 and p-value <0.001). Univariate Cox regression analysis was performed to screen for RMlncRNAs significantly associated with OS (p < 0.05), and the Sankey diagram was mapped by the “ggalluvial” R package. The Wilcoxon test was used to detect differences in the expression of prognosis-related RMlncRNAs between tumor tissues and normal tissues.

Based on the expression of RMlncRNAs with prognostic value, unsupervised consensus clustering was performed using “ConsensusClusterPlus” on 433 colon cancer patients to identify potential molecular subtypes (Wilkerson and Hayes, 2010). R packages “ survival” and “survminer” were used to analyze the prognosis of samples with different molecular subtypes. Clinical data were included and analyzed for differences in molecular subtypes by using the “heatmap” R package for distinct clinicopathological features. The proportion of 22 tumor-infiltrating immune cells (TICs) in each sample was quantified using the CIBERSORT algorithm (Newman et al., 2015). The ESTIMATE algorithm was used to calculate the tumor microenvironment (TME) score (including immune score, stromal score, ESTIMATE score, and tumor purity) for each sample (Yoshihara et al., 2013). In addition, we synthesized 38 immune checkpoint genes from the literature and examined the expression of these checkpoint genes among molecular subtypes (Pardoll, 2012; Nirschl and Drake, 2013).

The TCGA-CORD cohort was randomly divided into a training set and a test set (1:1 ratio). A minimum absolute shrinkage and selection operator (LASSO) Cox regression analysis was used to narrow down candidate lncRNAs and develop an RNA methylation-related lncRNA signature (we named it RMlnc-score). The formula is as follows: RMlnc-score = Σ (βi × Expi) (β: coefficients, Exp: lncRNA expression level). Patients were then divided into high RMlnc-score and low RMlnc-score groups based on the median value of RMlnc-score. Kaplan-Meier survival curves were plotted using the R package “survival” to describe the overall survival difference between the high and low score groups. Receiver operating characteristic curves (ROC) analysis was performed to evaluate its sensitivity and accuracy. Heatmaps were generated to reveal differences in signature lncRNA expression in the low and high RMlnc-score groups.

The student’s t-test was used to assess the relationship between RMlnc-score and clinical characteristics. In addition, survival analysis was performed to further elucidate the relationship between RMlnc-score by sex (male and female), age (≤65 and >65 years), T-stage (T1-2 and T3-4), N-stage (N0 and N1-2), M-stage (M0 and M1), and grade (stages I-II and stages III-IV) in each subgroup for prognostic ability. Subsequently, univariate and multivariate Cox regression analyses were used to determine the relationship and independence between clinicopathological characteristics and RMlnc-score. A nomogram and calibration curves were then constructed based on independent prognostic factors from multivariate Cox regression analysis to predict the probability of survival at 1, 3, and 5 years in CC patients. The GSE17536 dataset was used as an external validation cohort to further assess the prognostic value and clinical relevance of model lncRNAs.

The R package “scatterplot3d” was used to perform PCA analysis to explore potential differences between high and low RMlnc-score groups. To analyze the correlation between RMlnc-score and TICs, we used different software (including ssGSEA, xCELL, Timer, Quantiseq, MCPcounter, EPIC, CIBERSORT-ABS, and CIBERSORT) to comprehensively analyze of immune cell infiltration.

The tumor immune dysfunction and exclusion (TIDE) algorithm (Fu et al., 2020) was used to assess the potential response of colon cancer patients in the different RMlnc-score groups to immunotherapy. Data from the Genomics of Drug Sensitivity in Cancer (GDSC) database were used to predict the response of CC patients to chemotherapeutic drug therapy. The “pRRophetic” R package (Geeleher et al., 2014) was used to calculate the half-maximal inhibitory concentration (IC50) of common chemotherapeutic agents.

m6A-Atlas (Tang et al., 2021) and SRAMP(Zhou et al., 2016) were used to predict the m6A site of the lncRNAs; m5C-Atlas (Ma et al., 2022) and RNAm5Cfinder (Li et al., 2018) were used to predict the m5C site of the lncRNAs; m7GHub (Song et al., 2020) and iRNA-m7G (Chen et al., 2019) databases were used to predict the m7G site of the lncRNAs.

All statistical analyses were performed using R software (v4.0.2). p values < 0.05 were considered statistically significant if not explicitly stated.

The workflow of this study is illustrated in Figure 1. First, we investigated the expression of 50 m1A-, m5C-, m6A-, and m7G-RNA methylation regulatory genes in the TCGA-CORD cohort (Figure 2A). The results showed that there were 42 differentially expressed RNA methylation regulatory genes. Among them, 37 regulators were highly expressed in colon cancer tissues, and five were lowly expressed in colon cancer tissues. Next, we investigated the incidence of somatic mutations and copy number variations for 50 regulators in TCGA-CORD. A total of 140 of 399 samples (35.09%) experienced genetic alterations in RNA methylation regulators (Figure 2B). Among them, ZC3H13 (9%) was the gene with the highest mutation frequency, followed by YTHDC2 (6%), and RBM15 (5%.) The investigation of CNV alteration frequency revealed that all RNA methylation regulators were found to show prevalent CNV alterations. Among them, DNMT3B, ALYREF, YTHDF1/3, IGF2BP2/3, YBX1, and HNRNPA2B1 showed significant copy number amplification, while TRMT6, YTHDF2, YTHDC2, and RBM15/15B showed remarkable copy number deletions (Figure 2C). Figure 2D shows the location of CNV changes in RNA methylation regulators on chromosomes. The above analysis revealed a high degree of heterogeneity in the expression and inherited variation status of RNA methylation in CC, demonstrating that RNA methylation-related regulators may play a pivotal position in the occurrence and development of CC.

We identified 1,057 lncRNAs significantly associated with 42 differentially expressed RNA methylation regulators by using Pearson correlation analysis and defined them as RMlncRNAs. Based on the mRNA-lncRNA co-expression pattern, we constructed a Sankey diagram to show their linkage (Figure 3A). After excluding normal tissues or patients lacking survival data, we merged survival information with RMlncRNA expression data of colon cancer patients (final number of patients = 433). Subsequently, we performed univariate Cox regression analysis and found that 23 RMlncRNAs were significantly associated with OS of colon cancer patients (p < 0.05, Figure 3B). Among them, only TNFRSF10A-AS1 was identified as a protective factor with a risk ratio (HR) < 1, while all others were considered as risk factors. The bar graph and heatmap showed significant differences in the expression of these 23 prognosis-related RMlncRNAs between normal and colon cancer tissues (Figures 3C, D).

Based on the expression levels of 23 prognosis-related RMlncRNAs in CC samples, we clustered 433 samples by an unsupervised clustering approach to further elucidate the biological differences between subgroups. Our results showed that K = 2 was the optimal number of clusters with the highest correlation within groups and the least interference between groups (Figures 4A–C). Therefore, CC patients were divided into two subgroups: Cluster1 (n = 170) and Cluster2 (n = 263). The survival analysis results showed a significant survival advantage for Cluster2 patients (p = 0.021, Figure 4D). The heatmap showed differences in prognosis-related RMlncRNA expression between subgroups (Figure 4E), and most RMlncRNAs were highly expressed in Cluster1. In addition, we found that patients with distant metastasis (M1) were more represented in Cluster1 (p < 0.05), while other clinicopathological features were not significantly different between the two subgroups.

We further explored the differences in immune microenvironment characteristics between distinct clusters to understand the interactions between RNA methylation-related lncRNAs and the immune microenvironment (TME). The results of CIBERSORT analysis (Figure 5A) showed that 8 of the 22 immune infiltrating cells differed between clusters, with T cells CD8, T cells regulatory (Tregs), NK cells resting, NK cells activated, monocytes, dendritic cells resting, and neutrophils showed more infiltration in Cluster2, while only T cells CD4 memory activated were highly enriched in Cluster1. The percentages of 22 immune cell types in GC patients between the two clusters are shown in Figure 5B. ESTIMATE analysis showed (Figure 5C) that the immune score (p < 0.001), stromal score (p = 0.0062), and ESTIMATE score (p < 0.001) were significantly higher in Cluster2 than Cluster1, while the tumor purity in Cluster1 (p < 0.001) was considerably higher than Cluster2. In addition, we tried to determine the correlation between subgroups and some immune checkpoints. We found remarkable differences in the expression levels of 18 immune checkpoint genes between the two subtypes (p < 0.05). The expression levels of PD-1, PD-L1, HAVCR2, CTLA4, LDHA, LGALS9, TNFRSF18, YTHDF1, LAG3, CD40, TNFRSF4, TNFRSF9, CD86, B2M, and CD8A were higher in Cluster2, whereas PDCD1LG2, IL12A, PVR, and JAK1 were higher in Cluster 2 (Figure 5D). Previous studies have shown that high immune scores and activation of suppressive immune checkpoints (like HAVCR2, PD-L1, CTLA-4) play a crucial role in “hot tumors” (Zhan et al., 2021). “Hot tumors” are more likely to benefit from immune checkpoint blockade (ICB) therapy, whereas “cold tumors” with low levels of immune infiltration are more likely to become resistant to immunotherapy (Galon and Bruni, 2019). Therefore, we may consider cluster 1 as the “cold tumor” and cluster 2 as the “hot tumor”, which may predict different immunotherapy responses.

The 433 colon cancer patients were randomly divided into a training set (n = 217) and a test set (n = 216). To avoid overfitting, we screened the seven most powerful prognostic RMlncRNAs by LASSO regression analysis, which were used to construct the RNA methylation-related lncRNA prognostic signature (RMlnc-score) (Figures 6A, B). The correlation coefficients are shown in Table 1. Patients were classified into low RMlnc-score and high RMlnc-score groups according to the cut-off values of RMlnc-score. The RMlnc-score for each patient was calculated as follows:RMlnc-score=(0.0645*ALMS1-IT1 expression) + (−0.1268*TNFRSF10A-AS1 expression) + (0.6464*FRMD6-AS1 expression) + (0.6173*STARD7-AS1) + (0.4430*LINC02257 expression) + (0.2254*AP001505.1 expression) + (0.2329*AC019205.1 expression). The Kaplan-Meier curves showed that in the training set (p < 0.001, Figure 6C) and test set (p = 0.002, Figure 6D), patients in the high RMlnc-score group had a worse prognosis compared to the low RMlnc-score group. The area under the curve (AUC) for 5-years overall survival (OS) was 0.741 and 0.734 for the training and test sets, respectively (Figures 6E, F). In the overall cohort (Figure 6G), the RMlnc-score (our study) had an AUC of 0.737 at 5-years overall survival, which was substantially higher than ChaiLncSig (AUC = 0.653), YunLncSig (AUC = 0.658), and ZhangLncSig (AUC = 0.659). This suggests that the RMlnc-score has higher accuracy in predicting survival compared to three recently published lncRNA signatures for colon cancer (Chai et al., 2021; Yun and Yang, 2021; Zhang et al., 2021d). The survival status and RMlnc-score score curves for the training and test sets showed (Figure 6H, I, 7I) that RMlnc-score was proportional to the number of deaths in CC patients. The heatmaps showed (Figures 6J, K) that the expression of ALMS1-IT1, FRMD6-AS1, STARD7-AS1, LINC02257, AP001505.1, and AC019205.1 was upregulated in the high RMlnc-score group, while TNFRSF10A-AS1 was upregulated in the low RMlnc -score group was up-regulated. In addition, we performed a validation analysis of the signature lncRNA in the GSE17536 cohort. However, due to fewer non-coding genes in the microarray data, we only detected ALMS1-IT1 and FRMD6-AS1. Our results showed that high expression of ALMS1-IT1 (p = 0.044) and FRMD6-AS1 (p = 0.034) was significantly associated with poor prognosis of patients (Figures 7A, B). High expression of ALMS1-IT1 was associated with high grade (p = 0.0017, Figure 7C) and high expression of FRMD6-AS1 was associated with high stage (p = 0.021, Figure 7D) and high grade (p = 0.028, Figure 7E).

Stratified survival analysis in combination with clinical characteristics (Figures 8A–L) showed that in age>65 (p < 0.001), age ≤ 65 (p < 0.001), male (p < 0.001), female (p = 0.004), stage III-IV (p = 0.002), T3-4 (p < 0.001), M0 (p < 0.001), and N1-2 (p < 0.001) subgroups of patients, survival was significantly lower in the high RMlnc-score group than in the low RMlnc-score group. By comparing the RMlnc-score of patients in different groups, we found that RMlnc-score increased with increasing T-stage, N-stage, M-stage, and clinical stage, while no significant differences were seen for age and gender (Figures 8M–R). Univariate Cox regression analysis showed that age, stage, T-stage, N-stage, M-stage, and RMlnc-score (all p < 0.001) were strongly associated with prognosis (Figure 9A). Multivariate Cox regression analysis confirmed that age, T-stage, and RMlnc-score were independent prognostic factors for CC patients (Figure 9B). Based on the three independent prognostic factors in the multivariate Cox regression analysis, we created a nomogram capable of predicting the incidence of OS in CC patients at 1, 3, and 5 years (Figure 9C). The calibration curve demonstrated the high accuracy and sensitivity of this nomogram (Figure 9D).

The results of principal component analysis (PCA) showed no significant differences between the high RMlnc-score group and the low RMlnc-score group in the expression of all genes (Figure 10A), RNA methylation-related genes (Figure 10B), and RNA methylation-related lncRNAs (Figure 10C). However, in the expression of the seven lncRNAs used in the prognostic model (Figure 10D), there was a significant difference between the high RMlnc-score and low RMlnc-score groups. We also explored whether our model could predict immune cell infiltration in CC. The bubble plot (Figure 10E) showed that RMlnc-score was positively correlated with CD4⁺ T cells, cancer-associated fibroblast (CAFs), myeloid dendritic cell, macrophage M0, NK cell activated, hematopoietic stem cell while negative correlation with CD4+8 cell, monocyte, neutrophil, and B cell plasma. The ssGSEA results (Figure 10F) showed that some immune cells, including dendritic cells (DCs), activated dendritic cells (aDCs), immature dendritic cells (iDCs), mast cells, neutrophils, NK cells, and type 2 T helper were significantly increased in the low RMlnc-score group, and some pathways associated with immune function, namely APC co-stimulation, C-C chemokine receptor, and cytolytic activity, were significantly activated in the low RMlnc-score group.

TMB and MSI have been reported to be predictive biomarkers of immunotherapeutic response (Ock et al., 2017; Vandekerkhove et al., 2021). Therefore, we first compared somatic mutations in high RMlnc-score and low RMlnc-score and visualized the top 20 genes with the highest mutation frequency (Figures 11A, B). However, there was no significant difference in tumor mutational load between the high RMlnc-score and low RMlnc-score groups (Figure 11C). We then compared the differences in MSI distribution between the different scoring groups and found that the low RMlnc-score group was associated with higher microsatellite instability (MSI) (Figure 11D). TIDE, a novel predictive marker of immunotherapy, was better than known immunotherapy biomarkers (including TMB and PD-L1 expression) for response to immunotherapy in certain tumors (Wang et al., 2019). Higher TIDE scores indicate that tumor cells are more likely to induce immune escape, thus indicating a lower response rate to immunotherapy. Surprisingly, we found that patients in the low RMlnc-score group had significantly lower TIDE scores (including T cell dysfunction and exclusion scores) than those in the high RMlnc-score group (Figures 11E–G). The above findings suggested that RMlnc-score correlates with the response of CC patients to immunotherapy and may help predict the efficacy of ICB immunotherapy.

To explore the effect of RMlnc-score on drug response, we compared the half-maximal inhibitory concentration (IC50) of the commonly used drugs in both groups. The results showed that the IC50 values of bicalutamide, lapatinib, sorafenib, metformin, and temsirolimus were higher in the high RMlnc-score group, indicating that patients in the low-scoring group were more sensitive to these five drugs. In contrast, axitinib, bexarotene, bosutinib, elesclomol, embelin, etoposide, imatinib, lenalidomide, methotrexate, midostaurin, nilotinib, pazopanib, shikonin, vinblastine, vinorelbine, and vorinostat had higher IC50 in patients with low RMlnc-score, implying that patients in the high RMlnc-score group were more sensitive to these drugs (Figure 12).

After scanning the m6A-Atlas, m5C-Atlas, and m7GHub databases, we eventually obtained six m6A, nine m5C, and one m7G modification sites on STARD7-AS1 and five m5C modification sites on FRMD6-AS1, which have been experimentally validated (Supplementary Table S1). Then, we also utilized the widely used bioinformatics tools SRAMP, RNAm5Cfinder, and iRNA-m7G to predict potential m6A, m5C, and m7G modification sites on our seven lncRNAs. Some meaningful results showed that all seven lncRNAs were potentially methylated (Supplementary Table S2).