Characterizing m6A Modification Factors and Their Interactions in Colorectal Cancer: Implications for Tumor Subtypes and Clinical Outcomes

doi:10.21203/rs.3.rs-4712713/v1

Download PDF

Research Article

Characterizing m6A Modification Factors and Their Interactions in Colorectal Cancer: Implications for Tumor Subtypes and Clinical Outcomes

https://doi.org/10.21203/rs.3.rs-4712713/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background

The study aims to comprehensively combine colorectal cancer data cohorts in order to analyze the effects of various DNA methylation-coding genes on colorectal cancer patients. The annual incidence and mortality of colorectal cancer are very high, and there are no effective treatments for advanced colorectal cancer. DNA methylation is a method widely used to regulate epigenetics in the molecular mechanism study of tumors.

Method

Three single-cell cohorts GSE166555, GSE146771, and EMTAB8107, and five transcriptome cohorts GSE17536, GSE39582, GSE72970, and TCGA-CRC (TCGA-COAD and TCGA-READ) were applied in this study. 2 erasers (ALKBH5 and FTO), There are 7 writers (METTL3, METTL14, WTAP, VIRMA, RBM15, RBM15B, and ZC3H13) and 11 readers (YTHDC1, IGF2BP1, IGF2BP2, IGF2BP3, YTHDF1, YTHDF3, YTHDC2, and HNRNPA2) B1, YTHDF2, HNRNPC and RBMX, a total of 20 M6A regulators, were used as the basis of the dataset in this study and were applied to the construction of molecular typing and prognostic models. Drugs that are differentially sensitive in methylation-regulated gene-related prognostic models were identified using the ConsensusClusterPlus package, which was also used to identify distinct methylation regulatory expression patterns in colorectal cancer and to model the relationship between tissue gene expression profiles and drug IC50 values. Finally, TISCH2 assessed which immune cells were significantly expressed with M6A scores. The immunosuppression of M6A methylation is spatially explained.

Results

This study used data from 583 CRC patients in the TCGA-CRC cohort. Firstly, the mutation frequency and CNV variation frequency of 20 m6A modification-related factors were analyzed, and the corresponding histogram and heat map were drawn. The study next analyzed the expression variations between mutant and wild forms of the VIRMA gene and explored differences in the expression of these variables in tumor and normal tissues. In addition, the samples were divided into different subgroups by molecular clustering method based on m6A modification, and each subgroup's expression and clinicopathological characteristics were analyzed. Finally, we compared prognostic differences, tumor microenvironment (TME) characteristics, immune cell infiltration, and gene function enrichment among different subpopulations. We also developed a colorectal cancer m6A-associated gene signature and validated its prognostic effects across multiple cohorts. Finally, using single-cell RNA sequencing data, we confirmed that tumor cells show elevated expression of m6A-related gene signatures.

Discussion

This study explored the mutation frequency, expression differences, interactions, molecular clustering, prognostic effect, and association with tumor characteristics of m6A modification-related factors in CRC and validated them at the single-cell level. These results clarify the association between m6A alteration and colorectal cancer (CRC) and offer important insights into the molecular recognition and management of cancer.

Colorectal cancer (CRC)

m6A modification

Mutation frequency

CNV variation frequency

Differential expression

One of the major cancers that pose a threat to Chinese citizens' lives and health and have created a significant societal burden is colorectal cancer (CRC). As per the most recent data given by the National Cancer Center, in China in 2015, there were 387,600 new cases of colorectal cancer, which accounted for 9.87% of all malignant tumor incidences; also, colorectal cancer caused 187,100 fatalities, which made up 8.01% of all malignant tumor deaths[1]. The typical progression time from precancerous lesions to cancer is 5 to 10 years, which presents a critical window of opportunity for early diagnosis and therapeutic intervention [2]. The stage at diagnosis also has a significant bearing on the prognosis of colorectal cancer. Only 14% of patients with stage IV colorectal cancer who have distant metastases will survive for five years, compared to 90% of those with stage I colorectal cancer. Molecular immunotherapy, chemical drug therapy, surgery, drug treatment, and other techniques are currently the main forms of treatment for colon cancer [3]. However, we still do not have a better treatment for advanced colon cancer.

Several genetic and molecular alterations have occurred in the development and course of colorectal cancer. The two most significant ones are the ongoing activation of oncogenes and the ongoing deletion of tumor suppressor genes [4, 5]. Classical tumor suppressor genes such as P53, APC, PTEN, SMAD4, ACVR2, and BAX are found to be absent to varying degrees in colorectal tumors[6, 7]. For example, 70 ~ 80% of colorectal tumors have corresponding functional inactivation caused by somatic APC gene mutations, and the inactivation of the APC gene is an important initiation event for intestinal tumors to become cancerous. 40% of colorectal adenomas (> 4 cm) have KRAS gene mutations[8]. The RAS gene downregulates many signaling pathways, such as MAPK and PI3K. About 5–10% of colorectal tumors have B-raf mutations, and the RAS gene activates the B-raf gene to activate the MAPK effectors MEK1 and MEK2 to perform their function. About 60 ~ 70% of colorectal tumors have P53 gene mutations, and point mutations in this gene will seriously interfere with the normal cell cycle, thereby accelerating the growth and invasion of colorectal cancer cells[9, 10]. Therefore, gene sequencing technology is widely used to detect biomarkers in solid tumors and predict immunotherapy effects.

With no alteration to the DNA sequence, DNA methylation is a significant epigenetic change that covalently bonds a methyl group at carbon position 5 of cytosine of CpG dinucleotides in the genome. Organizations to control the expression of genes [11, 12]. This alteration may impact the gene's promoter region's accessibility, which may control the gene's transcriptional activity. Aberrant alterations in DNA methylation, such as the hypomethylation of tumor suppressor genes and the hypermethylation of oncogenes, are frequently seen in cancer. The development, occurrence, and prognosis of tumors are closely associated with these alterations. Understanding different methylation patterns in colorectal cancer is critical for diagnosis, treatment, and prognostic assessment. High levels of DNA methylation can often affect the expression of tumor suppressor genes by reducing the transcriptional activity of the genes, such as reducing expression or silencing expression, thereby indirectly inducing malignant tumors. High levels of DNA methylation of cancer-promoting genes usually inhibit the occurrence of malignant tumors by reducing the transcriptional activity of tumor-promoting genes. Therefore, by studying methylation patterns in colorectal cancer, potential therapeutic targets and prognostic markers can be discovered, providing an important basis for personalized treatment and clinical management.

In addition, with the advancement of technology, such as high-throughput sequencing and bioinformatics analysis, we can gain a more comprehensive understanding of the methylation profile in colorectal cancer and discover key methylation events related to tumor initiation and progression[11–13]. Thanks to this research, the molecular mechanisms behind colorectal cancer are better understood. Italso ooffers a crucial theoretical and practical foundation for the development of novel treatment approaches and techniques for estimating prognosis.

With the current development of computational biology, WGCNA is used to identify gene clusters with the same function, and machine learning methods are used to construct prognostic scores with significant prognostic effects[14]. Therefore, this study used three single-cell cohorts GSE166555[15], GSE146771[16], and EMTAB8107[17], Single-cell Hub 2 (TISCH2) database[18], and 5 transcriptome cohorts GSE17536[19], GSE39582[20], GSE72970[21] and TCGA-CRC (TCGA-COAD and TCGA-READ) for transcriptome and single-cell combination of colorectal cancer. Analyze and target the tumor heterogeneity of different methylation expression patterns in colon cancer and construct a prognostic model of colorectal cancer based on methylation-related regulatory genes. A new theoretical foundation and clinical practice guidelines for the customized treatment of colorectal cancer are anticipated to be provided by the development of this cutting-edge study, which is also anticipated to improve patients' health and quality of life.

CRC dataset sources and preprocessing

The databases of Genotype-Tissue Expression (GTEx), Gene Expression Omnibus (GEO), and The Cancer Genome Atlas (TCGA) were thoroughly searched for public gene-expression data and comprehensive clinical annotation. We removed patients without survival data from additional analysis. We also obtained information on somatic mutations from the TCGA database.

For CRC scRNA-seq datasets, we obtained three publicly free cohorts, including GSE166555^[15], GSE146771^[16], and EMTAB8107^[17], from the Tumor Immune Single-cell Hub 2 (TISCH2) database^[18]. The harmony approach was utilized to integrate three distinct scRNA-seq datasets that underwent thorough analysis in accordance with the Seurat process. Annotation of each cell cluster was performed by the SingleR R package^[22] and literature references.

Extraction of expression levels of m6A regulators

We summarized the intersect genes of from GSE17536^[19], GSE39582^[20], GSE72970^[21] and TCGA-CRC (TCGA-COAD and TCGA-READ), and merged the two TCGA datasets (termed as TCGA-CRC) using “sva” package. A total of 20 m6A regulators were extracted from one previous study, including 2 erasers (ALKBH5 and FTO), 7 writers (METTL3, METTL14, WTAP, VIRMA, RBM15, RBM15B, and ZC3H13) and 11 readers (YTHDC1, IGF2BP1, IGF2BP2, IGF2BP3, YTHDF1, YTHDF3, YTHDC2, HNRNPA2B1, YTHDF2, HNRNPC, and RBMX)^{[23, 24]}. Subsequently, for further investigation, unique m6A alteration patterns were found by unsupervised clustering analysis based on patient classification and the expression of 20 m6A regulators. The stability of the cluster numbers was determined using the consensus clustering algorithm. The "ConsensusClusterPlus"^[25] package was utilized for the above steps.

Gene set variation analysis (GSVA) and functional annotation

The evaluation of the differences in biological processes and pathways across patterns of m6A alteration was carried out by GSVA enrichment analysis using the "GSVA" R packages). From the MSigDB database, we downloaded the “c2.cp.kegg.v2022.1.Hs.symbols.gmt” and “c5.go.bp.v2022.1.Hs.symbols.gmt” gene sets and used them for GSVA analysis.

Differentially expressed genes (DEGs) identification between m6A distinct phenotypes

Based on the expression of 20 m6A regulators, the "ConsensusClusterPlus" program was utilized to identify m6A-related genes and divide patients into two different patterns of m6A alteration. The limma R package's empirical Bayesian method was used to determine DEGs between various modification patterns. Adjusted P value < 0.01 & |logFC| >1 were the criteria used to determine DEGs.

Next, the prognostic analysis was conducted using the univariate Cox regression model to identify genes associated with CRC prognosis. Based on these genes, ConsensusClusterPlus was repeatedly utilized to stratify the CRC cohorts into two groups (termed gene clusters) with distinct features.

Generation of m6A gene signature

In order to quantify the m6A modification patterns of individual tumors, a set of scoring criteria was created to evaluate the m6A modification patterns of individual CRC patients. The gene signature for the m6A was named m6Ascore. The m6A gene signature was established using the following protocols: The two gene Clusters-derived prognostic DEGs were identified using the limma and the univariate Cox regression model. Genes with significant prognoses were extracted for additional analysis. A prognostic signature was then developed by multivariate Cox regression analysis and the Least Absolute Shrinkage and Selector Operation (LASSO) ^[26] study. Each CRC's risk score was determined using the formula: Risk score = ExpGene1 × CoefGene1 + ExpGene2 × CoefGene2 + … ExpGene(n) × CoefGene(n). In this equation, "ExpGene" stood for gene expression and "CoefGene" for the regression coefficient.

Scoring of the m6A signature at the single-cell level

To calculate enrichment scores of m6A signature in single cells, the AddModuleScore() function of the Seurat package^[27] was performed. The visualization was completed using the SCP R package.

Drug sensitivity analysis

The "OncoPredict" R package was used to analyze the targeted and chemotherapeutic drugs ^[28]. OncoPredict simulates the relationship between drug IC50 values and tissue gene expression patterns taken from the Broad Institute Cancer Cell Line Encyclopedia and the Genomics of Drug Sensitivity in Cancer, respectively. From the Cancer Immunome Atlas (https://tcia.at/home), the immunophenotype scores (IPS) of CRC patients were obtained. The IPS was used to evaluate each patient's response to immune checkpoint inhibitor-based immunotherapy (anti-CTLA-4 and anti-PD-1).

Next, CRC sample groups with high and low risk were identified based on the median risk score value. The log-rank test of the KM survival curve and the time-dependent receiver operating characteristic (ROC) curve were used to assess the prognostic potential and effectiveness of the signature.

Statistical analysis

R 4.1.3 software was used for statistical analysis in the study. For comparisons between two groups of continuous and categorical data, the independent Student's t-test or the Mann-Whitney U test were used; for categorical variables, the χ2 test was used. For normally distributed variables, Pearson's correlation test was used in correlation analysis; for non-normally distributed variables, Spearman's correlation test was utilized. Statistical significance was defined as a two-tailed P-value < 0.05.

Landscape of m6A regulators’ genetic variation in CRC

We identified 20 m6A regulators, comprising 2 erasers, 7 writers, and 11 readers. In the beginning, we compiled copy number variations and somatic mutations in CRC regulators. 143 of the 583 CRC samples collected had m6A regulator mutations; the frequency of these mutations was determined to be 24.53%. We found that in CRC samples, the one that demethylases ALKBH5 had the low mutation, whereas ZC3H13 and VIRMA showed the highest frequency of mutation (Fig. 1a). Twenty regulators were discovered to have a common CNV modification when the frequency of CNV alteration was examined; copy number amplification accounted for the bulk of these modifications. On the other hand, the frequency of CNV deletion was distributed in ALKBH5, HNRNPC, METTL14, METTL3, RBM15, YTHDC1, YTHDC2, and YTHDF2 (Fig. 1b). Using the merged TCGA-CRC and GTEx dataset, Fig. 1c showed the expression patterns of the 20 m6A regulators, demonstrating the higher expression levels in CRC tumor tissues relative to CRC normal tissues. Similar findings were addressed by Fig. 1d.

20 regulators mediating patterns of m6A methylation modification

We utilized a m6A regulator network to depict the comprehensive landscape of patients with CRC, encompassing the connections, interactions, and intercorrelations of m6A regulators (Fig. 2a). The analysis focused on the difference in eraser gene expression between the wild and mutant forms of VIRMA since the writer gene showed a comparatively greater mutation frequency. The results indicated that mutant tumors exhibited increased expression of numerous m6A regulators, including IGF2BP3, METTL14, HNRNPA2B1, and ALKBH5, relative to wild-type tumors (Fig. 2b). Based on the expression of 20 m6A regulators, the R package of ConsensusClusterPlus was utilized to categorize patients with qualitatively distinct patterns of m6A modification. Ultimately, we were able to identify three distinct patterns of modification through unsupervised clustering, including Cluster 1 (referred to as C1) and Cluster 2 (referred to as C2) (Fig. 3a-c). Significant variations were also seen in the m6A transcriptional profile between the two different m6A modification patterns (Fig. 3d). The prognosis analysis conducted on the two primary subtypes of m6A modifications revealed a considerably notable survival benefit in the m6Acluster-C1 modification pattern (Fig. 4a).

TME cell’s infiltration characteristics in distinct modification patterns of m6A

Using GSVA enrichment analysis, the biological behaviors of these various m6A modification patterns were investigated. m6Acluster-C2 showed marked enrichment in KEGG pathways, including cell cycle, basal transcription factors, ubiquitin-mediated proteolysis, and mismatch repair (Fig. 4b). Additionally, m6Acluster-C2 had substantially elevated activation levels for GO-BP terms, including RNA splicing, mRNA processing, control of mRNA processing, and mRNA metabolic process (Fig. 4c).

Decoding the TME contexture by the CIBERSORT method unraveled that m6Acluster-C2 was characterized by highly infiltrated memory B cells, M0 macrophage, resting NK cells, and activated memory CD4 + T cells. Contrarily, m6Acluster-C1 was featured by enriched regulatory T cells and plasma cells (Fig. 4d).

Generation of m6A gene signatures

Each m6A modification pattern's potential biological function was examined by identifying m6A phenotype-associated intersect DEGs using the limma program. Subsequently, a univariate Cox regression prognostic analysis was performed for every DEG. The 20 genes with a significant prognosis are shown in Fig. 5a. We next used unsupervised clustering analysis to classify patients into distinct genomic subgroups based on the 20 m6A phenotype-associated genes as a guide. Results from the unsupervised clustering approach were in line with the m6A modification patterns clustering grouping. Additionally, two distinct m6A alteration genomic phenotypes, known as m6A gene clusters 1 and 2, were identified by these results (Fig. 5b-d). This indicates that m6A methylation modification occurs in CRC in two different patterns. In m6A gene cluster 2, the 20 m6A regulators were upregulated (Fig. 5e). Additionally, it was demonstrated that CRC patients in gene cluster 1 had a better prognosis than those in gene cluster 2, which had a worse prognosis.

Generation of m6A score

Using LASSO + multivariate Cox regression analysis, we generated a m6A-related signature and the m6A score based on the prognostic DEGs between the two m6A gene clusters for each CRC patient (Figs. 6a and b). The four distinct cohorts' prognostic prediction performance for the m6A scores was shown using KM survival analysis and time-dependent ROC analysis (Figs. 6c-f).

Using R 4.3.1, a scoring system known as m6Ascore was created to examine the complexity and individual heterogeneity of m6A alteration. The m6A alteration pattern in CRC patients is evaluated by this approach using principal component analysis (PCA).Based on their m6Ascore, there was a noticeable difference between the m6A clusters. Cluster 2 exhibited a lower score, but Cluster 1 displayed a higher score in m6A (Fig. 7a). The fact that compared to m6A gene cluster 1, m6A gene cluster 2 had significantly increased m6Ascore (Fig. 7b).

m6A modification’ characteristics in tumor somatic mutation

Subsequently, we examined the correlation between m6Ascore and microsatellite instability (MSI). The findings showed that, in comparison to the low m6Ascore group, the high m6Ascore group had a higher proportion of MSI-H, which is linked to a better prognosis (Fig. 7c). Next, the TCGA-CRC cohort's differences in somatic mutation distribution between the high and low m6Ascores were examined using the maftools package. The high m6Ascore group showed a greater tumor mutation burden than the low m6Ascore group (Fig. 7e, f). The substantial correlation between greater TMB and high m6Ascore tumors was validated by the TMB quantification analyses (Fig. 7d).

Characteristics of m6A modification in the immunotherapy

Subsequent findings revealed that the high and low m6Ascore groups had varied PD-L1 expression patterns, with the former showing high expression (Fig. 8a), suggesting a possible response to anti-PD-1/L1 therapy. In order to determine whether the patient's reaction to ICB could be predicted, the signature of the m6A alteration was investigated. We obtained CRC immunophenoscore (IPS) data from The Cancer Immunome Atlas (TCIA; https://tcia.at/home). When compared to patients from the high m6Ascore group who received anti-CTLA4 and/or anti-PD-1 medication, analysis revealed that patients from the low m6Ascore group had better results (Fig. 8b).

As previously reported, distinct sensitivities to ICB may be associated with obviously different TME features. Resting NK cells, M0 macrophages, and naïve B cells were more abundantly infiltrated in CRC patients with high m6Ascore, while resting mast cells and resting NK cells were more abundantly infiltrated in patients with low m6Ascore (Fig. 9a). Furthermore, investigation of the biomarkers of T cell exhaustion (Tex) and M2 macrophages showed that high m6Ascore was associated with elevated levels of these markers (Fig. 9b, c).

Characteristics of m6A modification in cancer hallmarks

To uncover the dysregulated pathways between m6A phenotypes, we identified significantly upregulated genes in m6A gene cluster 2 using the limma R package. By utilizing the GSEA method, we found that hypoxia, TGFβ signaling, hedgehog signaling, KRAS signaling, TNFA-NFκB signaling, and estrogen response pathways were activated in the high m6A score group. Fatty acid metabolism, oxidative phosphorylation, and MYC target V1 were downregulated in the group with a high m6A score (Fig. 10).

Characteristics of m6A modification in the chemotherapy

It was looked into whether the signature of the m6A modification may be used to predict a patient's sensitivity to chemotherapy. Analysis showed that patients with low m6A scores would be better candidates for treatment with vinblastine, paclitaxel, docetaxel, and cisplatin than those with high m6A scores (Fig. 11a). Further analysis of the mRNAsi index revealed that CRC patients with high m6A scores had significantly lower mRNAsi than those with low m6A scores (Fig. 11b).

m6A modification’ characteristics at single-cell levels

To elucidate the features of m6A signature at CRC single-cell levels, we obtained and integrated three datasets with 110802 cells (Fig. 12a). Under the resolution 0.6, 26 cell clusters were identified and annotated as 17 major cell sets (Fig. 12b, c). Markers of each cell set are listed in Fig. 12d. By AddMouduleScore() function, we qualified the m6A signature scores of each CRC cell and found that the score was elevated in the malignant cell set (Fig. 12e).

Colorectal cancer is the second most common cause of death connected to cancer ^[29]. Despite significant advances in treatment modalities over the past few decades, such as refinement in surgical techniques, radiotherapy, chemotherapy, targeted therapy, and immunotherapy, the efficacy of treatments remains constrained by the lack of early detection rates and personalized therapeutic strategies^[30–35]. Further research into the mechanisms underlying the formation and progression of colorectal cancer is required in order to provide a theoretical foundation for diagnosis, therapy, and prognostic evaluation.

We retrieved data on various types of CRC from several sources in order to examine the expression patterns of 20 m6A regulators. Regarding copy number variations, most regulators exhibited amplification. Because most of the regulators in the tumor group have high expression levels, we hypothesize that copy number amplification plays a role in the higher expression levels of these regulators in tumors. Additionally, various regulators showed different degrees of mutation. VIRMA exhibited a relatively high mutation frequency, and in its mutated group, expression of other m6A regulatory factors was upregulated. Studies have suggested that writers form complexes that act synergistically, with VIRMA serving as a scaffold protein that brings together WTAP, HAKAI, and ZC3H13, creating a binding site for METTL3 and METTL14, thereby facilitating their optimal catalytic function^[36]. We hypothesize that mutations and copy number variations in VIRMA promote its overexpression and enhance its ability to interact synergistically with other proteins in the complex, although further research is needed to confirm our conjecture.

We divided all patients into two groups according to the m6A regulator expression patterns using unsupervised clustering techniques. Cluster 1 (C1) exhibited significant survival advantages. Upon thoroughly comparing the similarities and differences between the two clusters, we observed significant enrichment and activation of pathways related to proliferation, protein metabolism, RNA processing, and others in Cluster 2 (C2). This suggests that cancer cells in C2 are more active and prone to progression. Immune infiltration analysis revealed a less active immune cell infiltration status in C2 than C1, indicating poor tumor immune levels. Integrating these results, we propose that active proliferation and immune-suppressive microenvironment may contribute to the inferior prognosis of C2 compared to C1.

Next, we removed two grouping groups of m6A modification patterns using unsupervised clustering techniques. In m6A gene cluster 2, 20 m6A regulatory factors were significantly upregulated, and patients in this cluster exhibited poorer prognosis.

We used LASSO + multivariate Cox regression analysis to create m6A-related characteristics and determine a m6A score for every CRC patient based on prognostic differentially expressed genes from the two m6A gene clusters. This score demonstrated promising prognostic prediction potential and diagnostic efficiency. Patients were split into groups based on their m6A score: high and low, with the median serving as the cutoff. Tumor mutation burden and MSI-H were increased in patients with elevated m6A scores. High levels of microsatellite instability result from the mismatch repair (MMR) system's inability to address replication errors in short repetitive DNA sequences (MSI-H). The MSI-H phenotype in colorectal cancer is linked to tumor proximity, a high density of infiltration by local lymphocytes, and a low rate of distant organ metastasis. Additionally, the prognosis for MSI-H colorectal cancer is better than that of MSS (microsatellite-stable, or MSS) cancer. However, other studies indicate that patients with MSI-H colorectal cancer respond less favorably to chemotherapy, potentially because of a higher rate of mutation that leads to the emergence of clones resistant to the drug. In the era of immunotherapy, MMR (MMRd) defects are increasingly recognized as hypermutators that contribute to enhanced antitumor immune responses, possibly augmented through checkpoint inhibition. Importantly, frameshift mutations generated in MMRd tumors result in alterations of the entire coding sequence downstream of the mutation site, thereby eliciting greater immunogenicity, whereas point mutations only potentially generate novel antigens at the mutation site. Consequently, rather than their total tumor mutation burden, it is hypothesized that the high frameshift mutation rate of MMRd tumors may be more responsible for the checkpoint inhibitors' notable pan-cancer efficacy ^[37]. Subsequent immune-related analyses corroborated this observation. Moreover, the types and abundance of immune cell infiltration varied across different m6A score levels, possibly contributing to differential responses to ICB. Therefore, further investigation into the differences in the immune microenvironment under different m6A modification levels is warranted to design personalized treatment strategies for optimal ICB efficacy.

Our m6A score also predicts patients' sensitivity to chemotherapy. Patients with lower m6A scores may be preferred candidates for treatments such as oxaliplatin, paclitaxel, docetaxel, and cisplatin. Our scoring system provides guidance for personalized chemotherapy regimen selection. More extensive investigations with larger sample sizes are necessary to further validate our findings. Additionally, patients with higher m6A scores exhibited significantly lower mRNA stemness index (mRNAsi), indicating higher differentiation levels and relatively lower malignancy. Thus, our m6A score holds great potential in predicting tumor malignancy and guiding treatment strategies.

Our study still has certain limitations. Although our sample size is relatively large, there is still a gap between the sample size needed for application to the entire population. Furthermore, we lack certain real-world relevant studies. In the future, we will focus on addressing the issues identified in this study and further demonstrate the value of our m6A score in CRC. Collaborating with multiple medical centers to conduct multicenter, large-sample studies will further enhance the reliability and generalizability of our conclusions.

Conflict of Interest

Not application.

Funding

The research was funded by the Youth Doctoral Incubation Plan Project of the Second Affiliated Hospital of Army Military Medical University 2022YQB092

Author Contribution

WD Sun and ZQ Zhang wrote the main manuscript text and YC Su prepared figures 1-12. All authors reviewed the manuscript.

Data Availability

The data used in the present study are openly accessible in the GTEx, GEO, and TCGA databases.

Sun C, Liu Y, Huang Y, Li B, Rang W. Colorectal Cancer Incidence and Mortality Trends and Analysis of Risk Factors in China from 2005 to 2015. Int J Gen Med. 2021. 14: 9965-9976.
Binefa G, Rodríguez-Moranta F, Teule A, Medina-Hayas M. Colorectal cancer: from prevention to personalized medicine. World J Gastroenterol. 2014. 20(22): 6786-808.
Xie YH, Chen YX, Fang JY. A comprehensive review of targeted therapy for colorectal cancer. Signal Transduct Target Ther. 2020. 5(1): 22.
Fearon ER. Molecular genetics of colorectal cancer. Annu Rev Pathol. 2011. 6: 479-507.
Yamagishi H, Kuroda H, Imai Y, Hiraishi H. Molecular pathogenesis of sporadic colorectal cancers. Chin J Cancer. 2016. 35: 4.
Hasbullah HH, Musa M. Gene Therapy Targeting p53 and KRAS for Colorectal Cancer Treatment: A Myth or the Way Forward. Int J Mol Sci. 2021. 22(21).
Liebl MC, Hofmann TG. The Role of p53 Signaling in Colorectal Cancer. Cancers (Basel). 2021. 13(9).
Malumbres M, Barbacid M. RAS oncogenes: the first 30 years. Nat Rev Cancer. 2003. 3(6): 459-65.
Testa U, Pelosi E, Castelli G. Colorectal cancer: genetic abnormalities, tumor progression, tumor heterogeneity, clonal evolution and tumor-initiating cells. Med Sci (Basel). 2018. 6(2).
Zhuang Y, Wang H, Jiang D, et al. Multi gene mutation signatures in colorectal cancer patients: predict for the diagnosis, pathological classification, staging and prognosis. BMC Cancer. 2021. 21(1): 380.
Kato I, Badsha KZ, Land S, et al. DNA/RNA markers for colorectal cancer risk in preserved stool specimens: a pilot study. Tumori. 2009. 95(6): 753-61.
Romero-Garmendia I, Garcia-Etxebarria K. From Omic Layers to Personalized Medicine in Colorectal Cancer: The Road Ahead. Genes (Basel). 2023. 14(7).
Loktionov A. Biomarkers for detecting colorectal cancer non-invasively: DNA, RNA or proteins. World J Gastrointest Oncol. 2020. 12(2): 124-148.
Wang Y, Yan K, Lin J, et al. CD8+ T Cell Co-Expressed Genes Correlate With Clinical Phenotype and Microenvironments of Urothelial Cancer. Front Oncol. 2020. 10: 553399.
Uhlitz F, Bischoff P, Peidli S, et al. Mitogen-activated protein kinase activity drives cell trajectories in colorectal cancer. EMBO Mol Med. 2021. 13(10): e14123.
Zhang L, Li Z, Skrzypczynska KM, et al. Single-Cell Analyses Inform Mechanisms of Myeloid-Targeted Therapies in Colon Cancer. Cell. 2020. 181(2): 442-459.e29.
Qian J, Olbrecht S, Boeckx B, et al. A pan-cancer blueprint of the heterogeneous tumor microenvironment revealed by single-cell profiling. Cell Res. 2020. 30(9): 745-762.
Han Y, Wang Y, Dong X, et al. TISCH2: expanded datasets and new tools for single-cell transcriptome analyses of the tumor microenvironment. Nucleic Acids Res. 2023. 51(D1): D1425-D1431.
Williams CS, Bernard JK, Demory Beckler M, et al. ERBB4 is over-expressed in human colon cancer and enhances cellular transformation. Carcinogenesis. 2015. 36(7): 710-8.
Marisa L, de Reyniès A, Duval A, et al. Gene expression classification of colon cancer into molecular subtypes: characterization, validation, and prognostic value. PLoS Med. 2013. 10(5): e1001453.
Del Rio M, Mollevi C, Bibeau F, et al. Molecular subtypes of metastatic colorectal cancer are associated with patient response to irinotecan-based therapies. Eur J Cancer. 2017. 76: 68-75.
Aran D, Looney AP, Liu L, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol. 2019. 20(2): 163-172.
Yu L, Wang L, Sun J, et al. N6-methyladenosine related gene expression signatures for predicting the overall survival and immune responses of patients with colorectal cancer. Front Genet. 2023. 14: 885930.
Su H, Wang Y, Li H. RNA m6A Methylation Regulators Multi-Omics Analysis in Prostate Cancer. Front Genet. 2021. 12: 768041.
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010. 26(12): 1572-3.
Tibshirani R. The lasso method for variable selection in the Cox model. Stat Med. 1997. 16(4): 385-95.
Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018. 36(5): 411-420.
Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform. 2021. 22(6).
Ferlay J, Colombet M, Soerjomataram I, et al. Estimating the global cancer incidence and mortality in 2018: GLOBOCAN sources and methods. Int J Cancer. 2019. 144(8): 1941-1953.
Zielińska A, Włodarczyk M, Makaro A, Sałaga M, Fichna J. Management of pain in colorectal cancer patients. Crit Rev Oncol Hematol. 2021. 157: 103122.
Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. 2021. 71(3): 209-249.
Morgan E, Arnold M, Gini A, et al. Global burden of colorectal cancer in 2020 and 2040: incidence and mortality estimates from GLOBOCAN. Gut. 2023. 72(2): 338-344.
Messersmith WA. NCCN Guidelines Updates: Management of Metastatic Colorectal Cancer. J Natl Compr Canc Netw. 2019. 17(5.5): 599-601.
Sánchez-Gundín J, Fernández-Carballido AM, Martínez-Valdivieso L, Barreda-Hernández D, Torres-Suárez AI. New Trends in the Therapeutic Approach to Metastatic Colorectal Cancer. Int J Med Sci. 2018. 15(7): 659-665.
van der Stok EP, Spaander M, Grünhagen DJ, Verhoef C, Kuipers EJ. Surveillance after curative treatment for colorectal cancer. Nat Rev Clin Oncol. 2017. 14(5): 297-315.
Uddin MB, Wang Z, Yang C. The m(6)A RNA methylation regulates oncogenic signaling pathways driving cell malignant transformation and carcinogenesis. Mol Cancer. 2021. 20(1): 61.
Danaher P, Warren S, Ong S, Elliott N, Cesano A, Ferree S. A gene expression assay for simultaneous measurement of microsatellite instability and antitumor immune activity. J Immunother Cancer. 2019. 7(1): 15.

No competing interests reported.

Download PDF

Editorial decision: Revision requested
02 Aug, 2024
Reviews received at journal
29 Jul, 2024
Reviews received at journal
25 Jul, 2024
Reviews received at journal
25 Jul, 2024
Reviewers agreed at journal
22 Jul, 2024
Reviewers agreed at journal
20 Jul, 2024
Reviewers agreed at journal
19 Jul, 2024
Reviewers invited by journal
19 Jul, 2024
Editor assigned by journal
19 Jul, 2024
Submission checks completed at journal
19 Jul, 2024
First submitted to journal
09 Jul, 2024

You are reading this latest preprint version

Characterizing m6A Modification Factors and Their Interactions in Colorectal Cancer: Implications for Tumor Subtypes and Clinical Outcomes

Status:

Version 1

Abstract

Background

Method

Results

Discussion

Figures

Introduction

Methods

CRC dataset sources and preprocessing

Extraction of expression levels of m6A regulators

Gene set variation analysis (GSVA) and functional annotation

Differentially expressed genes (DEGs) identification between m6A distinct phenotypes

Generation of m6A gene signature

Scoring of the m6A signature at the single-cell level

Drug sensitivity analysis

Statistical analysis

Results

Landscape of m6A regulators’ genetic variation in CRC

20 regulators mediating patterns of m6A methylation modification

TME cell’s infiltration characteristics in distinct modification patterns of m6A

Generation of m6A gene signatures

Generation of m6A score

m6A modification’ characteristics in tumor somatic mutation

Characteristics of m6A modification in the immunotherapy

Characteristics of m6A modification in cancer hallmarks

Characteristics of m6A modification in the chemotherapy

m6A modification’ characteristics at single-cell levels

Discussion

Declarations

Conflict of Interest

Funding

Author Contribution

Data Availability

References

Additional Declarations

Status:

Version 1