Identification of m6A-Related Biomarkers Associated with Prognosis of Colorectal Cancer

Background

Colorectal cancer, also known as bowel cancer, colon cancer, or rectal cancer, is a cancer formed by uncontrolled cell growth in the colon or the rectum (part of the large intestine), or in the appendix. According to GLOBOCAN 2020 data, colorectal cancer (CRC) is the second most deadly and third most prevalent cancer in the world [1]. The incidence of CRC worldwide has been steadily increasing, accounting for 10% of all cancer diagnoses [2]. By 2020, CRC was the second most deadly cancer in the world, with an estimated 1.9 million new cases and 935 000 deaths. Meanwhile, the global burden of CRC is expected to increase by 60%, with nearly 2.2 million new cases and 1.1 million deaths expected by 2030 [3]. About 20% of CRC patients present with synchronous metastases, most commonly in the liver, and up to 60% of the patients develop distant metastases within 5 years [4]. Therefore, it is urgent to determine genetic and environmental risk factors that may affect the prognosis of CRC patients in clinical practice. Recently, with the help of high-throughput sequencing technology, a breakthrough has been made in the identification of CRC biomarkers at the cellular and molecular levels, which would potentially improve the prognostic prediction accuracy and introduce new therapeutic targets for CRC patients.

Since the discovery of first structurally modified nucleoside, pseudouridine, in the 1950s, more than 150 different chemical modifications have been identified on cellular RNA so far [5]. m6A methylation is the most characteristic mRNA modification and has been the most extensively studied since its discovery [6,7]. The most common methylation modification in eukaryotic mRNA is N6-methyladenosine (m6A), which accounts for more than 80% of all RNA base methylation and exists in various species [8]. The abundance and effects of m6A on RNA depend on the dynamic interplay between its methyltransferase (“writers”, such as METTL3, METTL14, WTAP, KIAA1429, ZC3H13, and METTL16), binding protein (“readers”, such as YTH domain-containing proteins and MRB1) and demethylase (“erasers”, such as FTO and ALKBH5) [9]. Accumulating evidence suggests that RNA methylation modification is also misregulated in human cancers and may be a potential ideal target for cancer treatment [10]. m6A methylation modification affects multiple aspects of RNA metabolism, ranging from RNA processing, nuclear export, RNA translation, to decay [11]. In addition, research proved that m6A methylation modification of mRNA and non-coding RNA plays an important role in a variety of common cancers, including solid tumors and non-solid tumors, and regulates cell proliferation and migration in cancer by affecting the biological functions of cells, tumor cell differentiation, and homeostasis [12,13]. Therefore, by finding m6A-related genes and m6A RNA methylation modification sites in cancer, new therapeutic targets could be provided for cancer treatment. Recently, studies focused on the role of m6A-related genes and their methylation regulators have revealed that METTL3 interacts with the microprocessor protein DGCR8 and actively regulates the pri-miR221/222 process in an m6A-dependent manner, which possibly has a carcinogenic effect in bladder cancer [14]. The reduction of RNA m6A methylation can activate oncogenic Wnt/PI3K-Akt signaling and promote malignant phenotypes of gastric cancer cells [15].

In the study, signature analysis with clinical information was performed for determining the m6A-related genes expression profile of patients with COAD (colon adenocarcinoma) and READ (rectum adenocarcinoma esophageal carcinoma) from TCGA database and Gene Expression Omnibus (GEO) databases. LASSO and Cox regression aided identification of potential m6A-related genes to predict the survival of patients with colorectal cancer. We determined 9 prognostic m6A-related genes from TCGA dataset and further validated the model in GEO dataset. The results obtained in this study would help to predict the prognosis of CRC patients and improve personalized treatment and management.

Material and Methods

DATASETS ACQUISITION:

RNA expression data and clinical information of primary CRC tissues were downloaded from TCGA database (https://portal.gdc.cancer.gov/). The dataset included survival data and FPKM expression value from 644 CRC tumor samples (478 COAD+166 READ). Normalized array data (GSE39582, GPL570 Affymetrix Human Genome U133 Plus 2.0 Array, France) and sample annotation files for 579 primary CRC tissues were obtained from Gene Expression Omnibus (GEO).

IDENTIFICATION OF M6A-RELATED GENE SET:

Currently known m6A RNA methylation regulators include methyltransferase (METTL3, METTL14, METTL16, WTAP, KIAA1429, RBM15, and ZC3H13), binding protein (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, and HNRNPC), and demethylase (ALKBH5 and FTO). There wer 3412 m6A-related genes correlated with CRC identified from the m6Avar database (http://rmvar.renlab.org/) [16]. The candidate m6A-related gene set was obtained by removing duplicate genes and the genes with no expression value or expression value less than 80% of total expression value in the samples. After refinement, 3161 candidate genes were kept for further analysis.

M6A-RELATED DEGS FROM PATIENTS WITH DIFFERENT PATHOLOGICAL STAGES:

To investigate the expression difference of m6A-related genes in different pathological stages, one-way ANOVA was used for identifying differentially expressed genes (DEGs) among 4 pathological stages. After removing samples without specific pathological stages, 620 tumor samples of CRC patients were left for further evaluation. Subsequently, all the differential expression analysis was executed with the threshold of P value <0.05. The heatmap which showed the expression differences of m6A-related genes in 4 different pathological stages was plotted using R package pheatmap.

SURVIVAL ANALYSIS:

Univariate Cox regression analysis was employed to identify m6A-related genes associated with the prognosis of CRC patients. According to the expression of genes, we classified the patients into high or low expression groups. With the survival time and survival status (live or dead) as input files, the overall survival (OS) probability between high or low expression groups was calculated and compared. OS of every CRC tumor sample was estimated by Kaplan-Meier method. The log-rank test was used to determine the significance of OS probability between different subgroups.

After that, the OS related genes were used for LASSO-Cox regression analysis to construct the prognostic model of these m6A-related genes with R package glmnet [23–25]. Ten-fold cross-validation minimum criteria were used to select the least-squares minimum value (min) with the minimum mean across validation error. Finally, we used to prognostic model to screen the gene set and calculate each patient’s risk score by a standard formula, which combines the expression levels of m6A-related genes with LASSO-Cox regression coefficients. To confirm whether the risk score independently affected the patients’ OS, the multivariate Cox regression model was used to assessed the association of pathoclinical features with the OS.

CONSTRUCTION OF PPI NETWORK:

The search tool for the retrieval of interacting genes (STRING database, V11; http://string-db.org/) was employed to predict the protein-protein interactions network of prognostic m6A-related DEGs [17]. On the STRING website, after put prognostic m6A-related DEGs, the website will show PPI network files according to its internal database. Subsequently, Cytoscape software (http://cytoscape.org/) was applied to visualize and analyze biological networks and node degrees of the 146 candidate genes based on a confidence score >0.4 [18].

GENE ONTOLOGY AND KEGG PATHWAY ENRICHMENT ANALYSIS:

Gene ontology (GO) is a tool for gene annotation using a dynamic, controlled vocabulary that classifies genes into 3 categories: biological process, molecular function, and cellular component [19,20]. GO analysis uses different genes to annotate gene functions based on the GO database. After obtaining all the functions involved in the genes, Fisher’s exact test and multiple comparison test are used to calculate the significance level (P value) and false positive rate (FDR) of each function, to screen out the significant functions embodied by the differentially expressed genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) database is used to assign gene sets to specific pathway maps of molecular interactions, reactions, and relation networks [21]. At present, KEGG Pathway is divided into 8 categories: overall network, metabolic processes, genetic information transmission, environmental information transmission, intracellular biological processes, biological systems, human diseases, and drug development. Pathway analysis was based on the KEGG database, using Fisher’s exact test and chi-square test for differentially expressed genes to analyze the significance of the pathway participated by the target gene.

GO functional annotation and KEGG pathway enrichment analyses of DEGs were performed by R package clusterProfiler [22] with P value <0.05 as statistically significant to further explore the functions and involved pathways of differentially expressed m6A-related genes.

STATISTICAL ANALYSIS:

All statistical analyses were conducted in R program (version 4.0.3; https://www.r-project.org/). Survival analysis was performed using Kaplan-Meier method and log-rank test with R package “survival” [26], while survival curves were plotted by “SurvMiner” package [27]. Cox proportional risk regression model was used for multivariate analysis. For all statistical tests, a P value of 0.05 was considered significant.

Results

IDENTIFICATION OF DIFFERENTIALLY EXPRESSED M6A-RELATED GENES FROM PATIENTS WITH DIFFERENT PATHOLOGICAL STAGES:

We selected 15 m6A RNA methylation regulators and 3412 m6A RNA methylation-related genes. After removing duplicate genes and the genes with no expression value or expression values less than 80% of total expression value in the samples, 3161 candidate genes were kept for further analysis. The workflow was shown in Figure 1. Subsequently, a total of 673 differentially expressed genes (DEGs) were identified by one-way ANOVA analysis among different stages (Figure 2). Notably, 4 m6A RNA methylation regulatory factors, including methyltransferase like 14 (METTL14), YTH domain containing 2 (YTHDC2), YTH N6-methyladenosine RNA binding protein 2 (YTHDF2), and zinc finger CCCH-type containing 13 (ZC3H13), were significantly differentially expressed among the 4 pathological subgroups. The expression profile of METTL14, YTHDC2, and YTHDF2 exhibited sustained decreasing with the progression of CRC, while the expression level of ZC3H13 was steadily elevated during CRC development (Supplementary Figure 1).

SCREENING OF PROGNOSTIC M6A-RELATED GENES:

To obtain prognostic m6A-related genes, univariate Cox regression analysis was performed with the threshold of P value <0.05. As a result, a total of 146 m6A-related genes were found to be associated with the overall survival of CRC patients (Supplementary Table 1).

PROTEIN–PROTEIN INTERACTION NETWORK OF THE M6A-RELATED DEGS:

The PPI networks for m6A-related genes were constructed according to the STRING database. The NetworkAnalyzer module was used to analyze and compare the network, hub genes were identified by analyzing the network topology. The top 10 hub genes were catenin beta 1 (CTNNB1), tripartite motif containing 37 (TRIM37), member RAS oncogene family (RAB7A), kinetochore scaffold 1 (CASC5/KNL1), centromere protein E (CENPE), G2/mitotic-specific cyclin-B1 (CCNB1), ubiquitin conjugating enzyme E2 H (UBE2H), heat shock protein family A (Hsp70) member 8 (HSPA8), kinesin family member 1A (KIF1A), and F-box and WD repeat domain containing 4 (FBXW4) (Figure 3).

FUNCTIONAL ANNOTATION OF M6A-RELATED GENES:

To further explore the biological functions and involved pathways of prognostic m6A-related genes, we used R package clusterProfiler to perform GO function annotation and KEGG pathway enrichment analysis. The most significant (P value <0.05) GO terms and KEGG pathways are shown in Figure 4. Results indicated that cerebral cortex neuron differentiation, forebrain development, digestive tract development, digestive system development, and odontogenesis of dentin-containing teeth were the top 5 significantly enriched GO terms in biological process category, while pathways of neurodegeneration-multiple diseases, small cell lung cancer, Amoebiasis, Tuberculosis and HIF-1 signaling pathway were the top 5 significantly enriched KEGG pathways.

CONSTRUCTION OF THE PROGNOSTIC RISK MODEL:

After 1000 resamples by Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression analysis through “glmnet” and “survival” R package, 9 m6A-related genes were selected to construct the prognostic model (Figure 5). Based on the gene signature with LASSO coefficients, the risk score=(0.02327*LRRC17)+ (0.8170*NFKB1)+ (0.07971*NOS2)+ (0.1193* PCDHB2)+ (0.006866*RAB7A)+ (0.07122*RPS6KA1)+ (−0.0897*RRNAD1)+ (0.3128*TLE6)+ (−0.06165*UBE2H). The TCGA cohort was then split into the high- and low-risk groups according to the median value of prognostic risk score. Kaplan-Meier survival analysis showed that patients with lower risk scores had significantly better overall survival than those with higher risk scores (P<0.0001) (Figure 5A). As shown in Figure 5B, the expression of RRNAD1 and UBE2H was negatively associated with the risk score of colorectal cancer patients, while the expression of other genes was positively associated with the risk score of patients with colorectal cancer. We also visualized the risk score distribution according to the length of follow-up months (Figure 5C). In the TCGA datasets, a significant positive correlation was found between the expression of NOS2 and TLE6 (Figure 5D). Figure 6A and 6B showed screening process of 9 m6A-related prognostic genes by LASSO-Cox regression analysis and random permutation. In Figure 6A, each line represents one gene, and the gene with non-zero coefficients was kept. In Figure 6B, when the line trends to flat, the lamda value would be chosen to conduct regression model. To verify whether these candidate prognostic gene models were influenced by clinical factors such as age, TNM stage and sex, we performed multivariate Cox regression analysis. The results showed that TNM stage, the status of vascular invasion, and the risk score were all independently related to the OS of TCGA CRC patients (Figure 6C).

VALIDATION OF THE PROGNOSTIC MODEL:

Concurrently, the prognostic risk score of the m6A-related signature was also validated in the GSE39582 dataset. In the exploration, the median risk score of 573 patients in GSE39582 was calculated as 2.188 according to the risk score formula, among which 286 patients with a score <2.188 were assigned as the low-risk group, while other 287 patients with a score ≥2.188 specified as the high-risk group. Patients with lower risk score had significantly longer overall survival than those with higher risk score in the validation set (P value=0.0082, Figure 7). Multivariate Cox regression analysis also showed that risk score was a prognostic factor for CRC patients’ OS (Figure 6D).

Discussion

Colorectal cancer is one of the most lethal solid tumors, with complex molecular and cellular heterogeneity. In the past few decades, there has been a great deal of research focusing on the molecular mechanisms of colorectal cancer, but most of them have concentrated on the aberration of protein-coding genes, leaving post-transcriptional processes mysterious. However, post-transcriptional alterations play a significant role in the preservation of tumor cells by modulating every hallmark in cancer [28,29]. RNA methylation modifications compose over 60% of all RNA modifications, and the most common type of RNA methylation modification is N6-methyladenosine (m6A) RNA methylation [30]. Accumulating evidence suggests that m6A modification plays a critical role not only in hypertension and cardiovascular disease, but also in tumor genesis and metastasis [31]. Therefore, the identification of m6A-related genes and m6A RNA methylation regulators abnormal expression may improve our understanding of colorectal cancer and provide us with valuable therapeutic targets.

Previous study about CRC and m6A methylation focused on m6A regulators [32], one subtype of CRC, such as COAD [33], or the comparison of tumor and normal tissues [34]. In the study, we systematically analyzed RNA sequencing data of CRC (including COAD and READ) patients from TCGA and GEO database. The m6A RNA methylation-related genes and m6A regulators were both brought into this investigation. We found that 673 candidate m6A RNA methylation-related genes were abnormally expressed among different pathological stages. Moreover, based on TCGA dataset, 146 m6A-related genes were related to patients’ overall survival based on the result of univariate Cox regression analysis. By constructing the PPI network, we identified 10 hub genes, namely CTNNB1, TRIM37, RAB7A, CASC5/KNL1, CENPE, CCNB1, UBE2H, HSPA8, KIF1A and FBXW4. CTNNB1 is a key component of the Wnt/β-catenin signaling pathway and was identified as the regulator of m6A modification in hepatoblastoma induced by METTL3 [35]. Liu et al found that Sec62 is upregulated by the METTL3-mediated m6A modification and promotes the stemness and chemoresistance of CRC by binding to β-catenin and enhancing Wnt signaling [36]. That might indicate that the WNT/β-catenin pathway should receive more attention in further study. A case-control study of breast cancer demonstrated that KIF1A promoter methylation can distinguish breast cancer (BC) cases from controls in plasma and was inversely associated with DNA repair ability (DRC) levels [37]. Studies confirmed that CCNB1 silencing can activate the p53 signaling pathway, further inhibit cell proliferation, and promote cell senescence in pancreatic cancer [38]. Meanwhile, in a systematic analysis of melanoma, it was proved to be positively correlated with either YTHDF1 or HNRNPA2B1, suggesting that both genes may affect m6A modification by CCNB1 gene [39].

To further investigate the influence of m6A RNA methylation regulatory factors on the prognosis of colorectal cancer, we used LASSO-Cox regression to establish a prognostic risk model based on 9 m6A RNA methylation-related genes. Survival analysis showed that the high-risk and low-risk subgroups classified by the model did have different prognostic destination in both the TCGA training group and the validation set. Although there has been no earlier research to prove that the biomarkers, we found in the prognostic model a close relationship with RNA modification process, but some of them (GSR and S1PR3) were correlated with the apoptosis and the adaption of acidic microenvironment of CRC cells [40,41]. Further experimental validation of the relationship between those biomarker genes and m6A RNA methylation regulators is needed to test the feasibility of our prognostic model.

Conclusions

In summary, our study systematically analyzed the expression profile of m6A RNA methylation-related genes, their prognostic significance, potential functions and pathways, and protein-protein interactions from CRC patients with the help of TCGA and GEO databases. Nine candidate m6A-related mRNA biomarkers (LRRC17, NFKB1, NOS2, PCDHB2, RAB7A, RPS6KA1, RRNAD1, TLE6, and UBE2H) were found to be closely related to the clinicopathology and prognosis of colorectal cancer. This study not only suggests the potential value of m6A-related genes as novel prognostic biomarkers for colorectal cancer, but also provides important clues for the diagnosis and treatment of colorectal cancer patients.