The gene expression data from the Cancer Genome Atlas (TCGA) (Cancer Genome Atlas Network, 2012) and the protein expression data from Clinical Proteomic Tumor Analysis Consortium (CPTAC) (Vasaikar et al., 2019) were downloaded from TCGA data portal (https://portal.gdc.cancer.gov/) and LinkedOmics (http://linkedomics.org), respectively. The expression levels of genes and proteins data were pre-normalized into the form of log2 (1 + Fragment Per Kilo-Million (FPKM)) and log2 intensity.

Prior to differential expression analysis, the genes with low abundance were excluded if their expression were below one FPKM in more than 85% of the total number of samples. The differential expression analysis was conducted using R limma method (Ritchie et al., 2015). The student t test and log2 fold change (FC) were applied to measure the differences. The genes with an adjusted p-value < 0.05 and a log2(FC) > 0.5 were considered as differentially expressed genes/proteins.

The consensus clustering algorithm was applied to determine the cancer subtypes using NMF method in R CancerSubtypes package (Xu et al., 2017). The optimal clustering number was determined according to the delta area plot.

The pairwise correlation between lncRNAs and m6A proteins were conducted using Spearman’s correlation analysis. Moreover, the physical interactions between the lncRNAs and m6A proteins were predicted by LncADeep (Yang et al., 2018b).

Two pathway enrichment methods were used for the data analysis, as they used different statistical testing approaches to analyze different types of data. The first is overrepresentation enrichment analysis (ORA), which used a given gene set as input, and identified Reactome pathways enriched by this gene set, while such gene sets consisted of significantly upregulated genes across the subtypes. The second is the gene set enrichment analysis (GSEA), which used all genes sorted by the statistics of interest in descending order, such as the correlation coefficient, and identified a gene subset significantly enriched in the head or tail of the ordered genes, along with significantly enriched pathways in this gene subset. These two analyses were implemented in R clusterProfiler package (Yu et al., 2012). The pathway activities were estimated by single-sample enrichment analysis (ssGSEA), which were implemented in R GSVA package (Hänzelmann et al., 2013).

The Cox proportional hazard regression model was employed in the survival analysis. The response variable is the survival time and status, and the predictive variables were the identified subtypes and expression levels of certain genes (high vs low, as compared to the median of gene expression). The survival analysis and visualization were implemented in R survival (https://cran.r-project.org/web/packages/survival/index.html) and survminer (https://cran.r-project.org/web/packages/survminer/index.html) packages.

The two-sample or pairwise mean comparison was conducted using the Mann-Whitney test. The two-sample or pairwise proportion comparison was conducted using Pearson’s chi-squared test. The p-value < 0.05, 0.01, and 0.001 were represented by the symbols ∗, ∗∗, and ∗∗∗.

To explore the expression patterns of m6A regulators and lncRNAs in COAD, we collected the gene expression data from Cancer Genome Atlas (TCGA) and protein expression data from the Clinical Proteomic Tumor Analysis Consortium (CPTAC). Moreover, we also collected a total of 23 m6A regulators from a previous study (Shen et al., 2021), including 8 writers (METTL3, METTL14, METTL16, RBM15, RBM15B, WTAP, KIAA1429, and ZC3H13), 3 erasers (ALKBH5, ALKBH3, and FTO) and 12 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, RBMX, and EIF3A) (Supplementary Data S1). Specifically, we identified that 10 m6A regulators including 3 writers (KIAA1429, ZC3H13, and RBM15) and 7 readers (YTHDF3, IGF2BP2, HNRNPA2B1, HNRNPC, RBMX, YTHDF1, and IGF2BP3) were upregulated in COAD samples at both mRNA and protein levels (Figures 1A,B, adjusted p-value < 0.05), as compared with the adjacent normal tissues. Moreover, we also identified 4060 differentially expressed lncRNAs (3,233 upregulated and 827 downregulated) (Supplementary Data S1). The correlation analysis between the DE-lncRNAs and DE-m6A regulators revealed 2,479 m6A-related lncRNAs (Spearman’s correlation >0.3 or < −0.3, Supplementary Data S1).

As tumor samples often exhibited high heterogeneity, we then investigated whether tumor samples could be divided into multiple subtypes based on those m6A-related lncRNAs. The consensus clustering analysis was conducted on the m6A-related gene expression data, which divided the tumor samples into 2–10 clusters, respectively. The delta area plot indicated that the optimal cluster number was 4, as cluster stability decreased significantly thereafter (Figure 2A). The correlation analysis across the tumor samples revealed that inter-cluster similarity was significantly lower, while the intra-cluster similarity was high. (Figure 2B). Notably, the COAD samples could be divided into four subtypes, each with 11, 200, 154, and 83 samples, respectively. Notably, the subtype 0 contained a small fraction of samples, all of which were formalin-fixed and paraffin-embedded (FFPE). Considering that the FFPE samples had different histological characteristics from the fresh samples, we excluded the subtype 0 in further analysis. Collectively, the m6A-related lncRNAs could clearly stratify the COAD samples into three subtypes.

With the m6A-lncRNA-related subtypes, we examined its association with the clinical characteristics of COAD patients. Particularly, the subtype-3 had the shortest overall survival when compared with the subtypes 1 and 2 (Figure 3A, log-rank test, p-value < 0.05). Consistently, the subtype-3 had a higher proportion of patients with advanced stage than the other two subtypes (Figure 3B). Moreover, we also found that nearly 40% of samples from the subtype 2 were with microsatellite instability (MSI), and such proportion was significantly higher than that observed in the two other subtypes (Figure 3C). In addition, other clinical factors such as patient age, KRAS mutation, perineural invasion, and lymphatic invasion were not significantly correlated to this subtyping results (Figures 3D–G). These results indicated that the m6A-lncRNA-related subtypes were prognostically relevant.

To further characterize the functionalities of the m6A-lncRNA-related subtypes, we conducted differential expression analysis and gene set enrichment analysis (GSEA) to identify subtype-specific genes and pathways. Specifically, we identified 4, 210, and 240 genes that were specifically upregulated in subtypes 1, 2, and 3, respectively. It should be noted that those genes were expressed higher in a specified subtype when compared with the other two subtypes Considering that there were few upregulated genes in subtype 1, we also identified 516 genes simultaneously upregulated in subtypes 1 and 3 when compared with subgroup 2. As shown in Figure 4A, those genes exhibited significantly different expression patterns across the subtypes.

The GSEA revealed that the genes upregulated in subtypes 1 and 3 were primarily enriched in NOTCH-related signaling pathways, while those specifically upregulated in subtype 3 were involved in energy metabolism such as oxidative phosphorylation, the citric acid (TCA) cycle and respiratory electron transport, respiratory electron transport, ATP synthesis by chemiosmotic coupling, and heat production by uncoupling proteins (Figure 4B). Notably, the subtype 2 was characterized by the inflammatory response-related pathways like interleukin-4 and interleukin-13 signaling, signaling by interleukins, interleukin-10 signaling, and neutrophil degranulation (Figure 4B). Accordingly, the components involved in NOTCH signaling pathways, such as NOTCH1 and NCOR2, inflammatory factors like IL18, IL1RN, IL1R2 and IFNGR1, and energy metabolism-related genes, such as D2HGDH, ATP6V1C2, RBL1 and LPIN3, were specifically upregulated in the corresponding subtypes (Figure 4C). These results indicated that the m6A-lncRNA-related subtypes had significantly different gene expression patterns and dysfunctional pathways.

As the m6A regulators were upregulated in COAD, we then investigated whether these genes were upregulated in any specific subtype. Among the ten m6A regulators upregulated in COAD, 8 genes including one writer (RBM15) and 7 readers (YTHDF3, IGF2BP2, HNRNPA2B1, HNRNPC, RBMX, YTHDF1, and IGF2BP3), were specifically upregulated in subtype 3 (Figure 5A, adjusted p-value < 0.05). Moreover, we also identified 21 m6A-related lncRNAs as subtype 3 specific, indicating that the m6A readers were closely associated with the biological characteristics of subtype 3 (Figure 5B, adjusted p-value < 0.05). More importantly, among the 21 m6A-related lncRNAs in subtype 3, four lncRNAs including CTD-3184A7.4, RP11-458F8.4, ANKRD10-IT1, and RP11-108L7.15 were found to be associated with poor prognosis (Figure 6), suggesting that the m6A regulators and the related lncRNAs were prognostically relevant. In addition, to test whether the m6A proteins could potentially regulate lncRNAs through RNA methylation, we predicted the physical interaction between the lncRNAs and m6a proteins using a deep learning method, LncADeep (Yang et al., 2018b). Notably, RP11-108L7.15 was predicted to interact with HNRNPA2B1 and RBMX, while CTD-3184A7.4 and RP11-458F8.4 might be regulated by YTHDF1 (Table 1), suggesting that the m6A proteins might act as the upstream regulators of the lncRNAs.

As the biological functions of the lncRNAs were usually unknown, to shed light on how the four prognostically relevant lncRNAs promoted the tumor progression, we conducted correlation analysis and GSEA. As the four lncRNAs were specifically upregulated in subtype 3, and this subtype was characterized by abnormal activity of energy metabolism, the GSEA further revealed that three out of these lncRNAs, including CTD-3184A7.4, RP11-458F8.4, and RP11-108L7.15 were positively correlated with the energy metabolism-related pathways such as oxidative phosphorylation, respiratory electron transport, the citric acid (TCA) cycle and respiratory electron transport (Figures 7A–C). In contrast, ANKRD10-IT1 was predicted to be involved in G alpha (12/13) signaling events and Rho GTPase cycle (Figure 7D).

Moreover, to further reveal the expression pattern of genes enriched in these pathways, we estimated the activities of these pathways based on the gene expression data using single-sample gene set enrichment analysis (ssGSEA). The survival analysis revealed that enhanced activities of these pathways were associated with shorter survival time (Figure 8), suggesting that high expression of the lncRNAs associated with these pathways might result in poor prognosis. Taken together, these results suggested that the four prognostically relevant m6A-related lncRNAs might participate in energy metabolism-related pathways, G alpha (12/13) signaling, or Rho GTPase cycle, thereby resulting in unfavorable outcomes.