Identification and Validation of m6A-GPI Signatures as a Novel prognostic Model for Colorectal Cancer

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

Download PDF

Research Article

Identification and Validation of m6A-GPI Signatures as a Novel prognostic Model for Colorectal Cancer

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

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

In order to develop an N6-methyladenosine-related gene prognostic index (m6A-GPI) that can predict the prognosis in colorectal cancer (CRC), we obtained m6A-related differentially expressed genes (DEGs) based on The Cancer Genome Atlas (TCGA) and m6Avar database, seven genes were screened by weighted gene co-expression network analysis (WGCNA) and least absolute shrinkage and selection operator (LASSO) analysis. Then, m6A-GPI was constructed based on the risk score. Survival analysis indicated that patients in the lower m6A-GPI group have more prolonged disease-free survival (DFS), and different clinical characteristic groups (tumor site and stage) also showed differential risk scores. In the analysis of the molecular characteristics, the risk score is positively associated with homologous recombination defects (HRD), copy number alterations (CNA), and the mRNA expression-based stemness index (mRNAsi). In addition, m6A-GPI also plays an essential role in tumor immune cell infiltration. The immune cell infiltration in the low m6A-GPI group is significantly higher in CRC. Moreover, we found that CIITA, one of the genes in m6A-GPI was up-regulated in CRC tissues based on real-time RT-PCR and Western blot. m6A-GPI is a promising prognostic biomarker that can be used to distinguish the prognosis of CRC patients in CRC.

Colorectal cancer

N6-methyladenosine

Prognostic signature

Immune infiltration

Colorectal cancer (CRC) has risen rapidly in recent years [1–3], and the treatment of CRC is mainly based on surgery, targeted therapy, neoadjuvant radiotherapy and adjuvant chemotherapy. Unfortunately, current treatments for colorectal cancer remain limited [4–6]. Therefore, we need to identify those high-risk CRC patients with poor prognosis, and further clarify the relevant mechanisms, so that individualized treatment can be implemented as soon as possible.

The tumor immune microenvironment (TIME) has an important impact on tumor prognosis and therapeutic effect [7, 8]. As the most abundant internal modification of eukaryotic mRNA and non-coding RNA, N6-methyladenosine (m6A) modification is not only associated with tumor growth, proliferation, and metastasis but also affects the process of immune cell recruitment and metabolic regulation of the TIME [9, 10], which will seriously affect the prognosis in CRC. Some studies have shown that blocking METTL3 can enhance the chemotherapeutic response and reduce stem cell frequency and tumor size both in vitro and in vivo [11]. Han et al. found that mice with blockade of YTHDF1 show an elevated antigen-specific CD8⁺ T cell antitumor response compared with wild-type mice, and the therapeutic efficacy of programmed cell death-ligand 1 (PD-L1) checkpoint blockade is enhanced in YTHDF1^−/− mice, which implies that YTHDF1 can serve as a potential therapeutic target in tumors [12]. In another study, the consumption of METTL3 or METTL14 in CT26 tumor mice with anti-PD-1 therapy significantly slowed tumor proliferation and prolonged the survival rate [13]. These results suggest that m6A modification may serve as a target affecting immune infiltration and survival time in CRC patients.

Despite this, there is no reliable tool for predicting the prognosis in CRC, and effective indicators are urgently needed. In this study, we developed a prognostic biomarker that can predict the prognosis and the immune infiltration in CRC patients. We focused on m6A-related genes, and seven genes were screened by weighted gene co-expression network analysis (WGCNA) and least absolute shrinkage and selection operator (LASSO) analysis. Then, we established an m6A-related gene prognostic index (m6A-GPI) based on the risk score. We conducted a series of stratification analysis and revealed the molecular and immune cell infiltration characteristics in the m6A-GPI subgroups. In addition, m6A-GPI was an independent predictor for CRC based on the univariate and multivariate Cox regression analysis results, and we constructed a nomogram including m6A-GPI to help clinicians accurately predict the prognosis of CRC patients. Moreover, we found that CIITA, one of the genes in m6A-GPI was up-regulated in CRC tissues based on real-time RT-PCR and Western blotting. These results showed that m6A-GPI is a reliable biomarker for predicting the prognosis of CRC.

1. Datasets Acquisition

We downloaded the RNA-seq data and clinical features of CRC from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov) database, which contained 616 tumor tissue samples and 51 paracancerous tissue samples. Mutation data were downloaded using the “TCGAbiolinks” packages in R language, and the independent validation datasets (GSE17538, 200 samples) were obtained from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/).

2. Identification Of M6a-related Genes

Then, pertinent references were searched and 21 m6A-related genes were screened, we identified 6,797 m6A-related genes in colorectal cancer from the m6Avar database (http://rmvar.renlab.org/). The differential expression analysis was performed in the "limma" R package, and differentially expressed genes (DEGs) were obtained in this process (adj. P < 0.05, log2FC > 0.585 or < 0.67). Consensus clustering, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were then performed in Metascape (https://metascape.org/) after consideration in the context of the m6A-related genes obtained from TCGA and m6Avar.

3. Wgcna Analysis

We conducted WGCNA analysis to identify hub genes. First, the similarity matrix was transformed into an adjacency matrix and then into a topological overlap matrix (TOM); TOM distances were used to cluster genes into WGCNA modules, and modules were determined by the dynamic pruning tree with a minimum of 30 genes per module.

4. Construction And Validation Of M6a-gpi

In order to identify prognostic genes, using the R package “survival” to perform univariate Cox regression analysis. Next, we implemented LASSO Cox regression to construct m6A-GPI that can predict the disease-free survival (DFS) of CRC patients, m6A-GPI was calculated based on the coefficient of genes, and the formula is:

$$ {r} {i} {s} {k} {s} {c} {o} {r} {e}=\sum _{ {i}= 1}^{ {n}}{ {C} {o} {e} {f}}_{ {i}} {*}{ {x}}_{ {i}}$$

where ${ {C} {o} {e} {f}}_{ {i}}$ is the coefficient, and the ${ {x}}_{ {i}}$ is the FPKM value of the m6A-related genes. We constructed Kaplan–Meier (K-M) survival curves of two subgroups and analyzed the 1-, 3-, and 5-year DFS rates of the cases. Validation datasets were downloaded from the GEO database, and we combined TCGA clinical information and explored the stability of m6A-GPI with different clinical characteristics.

5. Comprehensive Analysis Of Molecular And Immune Characteristics In Different M6a-gpi Subgroups

DNA changes are the basic factor in the development of cancer and play an important role in promoting the progress of cancer [14], and tumor immune escape mechanisms indicate that malignant tumors are capable of evading the immune response. To explore the immunogenicity of CRC, we analyzed the effect of m6A-GPI on mutation load, homologous recombination defects (HRD), neoantigen loads, copy number alterations (CNA) and the mRNA expression-based stemness index (mRNAsi). We obtained immune characteristics from the Genomic Data Commons (GDC) data portal (https://gdc.cancer.gov/about-data/publications/panimmune), and the HRD score comes from PMID: 29617664. Furthermore, we analyzed the somatic mutation difference between the low- and high-risk groups by the R package “maftools.”

The CIBERSORT (https://cibersort.stanford.edu/) is a novel algorithm and it can evaluate gene expression data from RNA sequences and assess the immune cell compositions of complex tissues [15, 16]. CIBERSORT can be used to calculate the content of 22 kinds of human immune cell phenotypes and the sum of all estimates of immune cell type fractions yields one. We compared the relative proportions of 22 immune cells between the two subgroups and presented the results in a landscape map.

In the tumor microenvironment, non-tumor components are divided into two types that are valuable for tumor diagnosis and prognostic evaluation, immune cells and stromal cells [17]. To determine the impact of immune cell infiltration (such as T cells, Tregs, NK cells and macrophages) on the treatment of ICIs, we calculated the immune score, matrix score, and tumor purity in each CRC sample based on the ESTIMATE algorithm.

6. Independent Prognostic Factor And Nomogram

To verify whether m6A-GPI can serve as an independent prognostic factor in CRC, we conducted univariate and multivariate Cox regression analysis. In order to provide doctors with a quantitative method for predicting the prognosis of patients with CRC, we constructed a nomogram using the risk status, age, cancer type, sex, cancer stage, and cancer site, then established calibration plots of DFS at 1-, 3-, and 5-year in the TCGA cohorts.

7. Crc Tissue Samples

Colorectal samples and their adjacent normal tissues were collected from Liaoning Cancer Hospital and Institute (Shenyang, China), and the colorectal samples showed a confirmed histological diagnosis of CRC. The study was approved by the institutional ethics committee, and individual consent forms were signed by each patient.

8. Quantitative Real-time Rt-pcr

Total RNA was extracted from colorectal tumors and their adjacent normal tissues of 16 patients by using Trizol reagent (Life Technologies, Carlsbad, CA) following the manufacturer recommendations. Concentration of RNA was quantified by Nanodrop 2000 (Thermo, Wilmington, DE). Reverse-transcription to cDNA (50 ng per sample) was using with iScript cDNA Supermix (TaKaRa, Dalian, China). Quantitative RT-PCR was performed using a reaction mixture containing SYBR mix (TaKaRa,Dalian,China), and real-time fluorescence was detected by QuantStudio 6 Flex (ABI, Foster City, CA). The primers were designed and synthesized by Life Technologies. The sequences of the primer pairs were as follows, GBP2: forward 5′-CTATCTGCAATTACGCAGCCT-3′, reverse 5′-TGTTCTGGCTTCTTGGGATGA-3′, CXCL10: forward 5′-GTGGCATTCAAGGAGTACCTC-3′, reverse 5′-TGATGGCCTTCGATTCTGGATT-3′, CXCL13: forward 5′-GCTTGAGGTGTAGATGTGTCC-3′, reverse 5′-CCCACGGGGCAAGATTTGAA-3′, FASLG: forward 5′- TGCCTTGGTAGGATTGGGC-3′, reverse 5′-GCTGGTAGACTCTCGGAGTTC-3′, CIITA: forward 5′-CCTGGAGCTTCTTAACAGCGA-3′, reverse 5′-TGTGTCGGGTTCTGAGTAGAG-3′, IL12RB1: forward 5′-TAGGGACCTGAGATGCTATCG-3′, reverse 5′-CCCGGAGCTAAGGCAACAC-3′, CXCR6: forward 5′-GACTATGGGTTCAGCAGTTTCA-3′, reverse 5′-GGCTCTGCAACTTATGGTAGAAG-3′, GAPGH: forward 5′-TCCCATCACCATCTTCCA-3′, reverse 5′-ACTCACGCCACAGTTTCC-3′.

9. Western Blot Analysis

CRC samples and their adjacent normal tissues of 7 patients were lysed according to the kit (PC101-PC104, Epizyme Biomedical Technology, Shanghai, China). Total protein concentration was quantified using a BCA protein assay kit (DQ111-01, TransGen Biotech Co., Ltd., Beijing, China). 50 µg of protein was separated by SDS–PAGE and was electrophoretically transferred to polyvinylidene difluoride (PVDF) membranes. The membranes were blocked with 3% BSA in Tris-buffered saline (TBS) containing 0.05% Tween-20 for 1 h. After blocking, the membranes were incubated with corresponding primary antibodies overnight at 4°C. The primary antibody for CIITA (#55099-1-AP; 1:1000) was purchased from Proteintech Group, Inc. (Wuhan, China). Membranes were washed by 1× TBST, followed by incubation with anti-Rabbit IgG-HRP (ZB-2301, Zhong Shan-Golden Bridge Biological Technology Co., Beijing, China) for 1 h. Immunoreactive bands were visualized by using Tanon 5500 (Tanon, Shanghai, China). Equal loading of proteins was verified by GAPDH (#60004-1-Ig; 1:3000, Proteintech Group, Inc., Wuhan, China).

10. Statistical analysis

SPSS 24.0 (IBM Corporation, Armonk, United States) and R programming language (version 4.0.2) were used to perform the statistical analysis. Kaplan–Meier curves and the log-rank test were used to compare the DFS between various subgroups. The prognostic ability of the predictors for 1-, 3-, and 5-year DFS was evaluated by receiver operating characteristic (ROC) curves and the area under the curve (AUC) values. Univariate and multivariate Cox regression analysis were utilized to evaluate the independent prognostic value of the model. A two-sided P < 0.05 was considered significant.

1. Identification of m6A-related genes

Based on our workflow (Fig. 1), we obtained 6,803 m6A-related and referred to as GeneSet 1. These genes intersected with 5,220 DEGs from the TCGA database, and 1,291 differentially expressed m6A-related genes were screened for GO and KEGG analysis. In the KEGG and GO analysis, we extracted a total of 5,642 genes from the significantly enriched pathways, which were referred to as GeneSet 2. In addition, we also clustered 1,291 differentially expressed m6A-related genes, and finally obtained two sets of samples, and screened out 644 DEGs between them, which were referred to as GeneSet 3. The three gene sets contained a total of 10,893 genes (Fig. 2A). The screening process is shown in Supplementary Fig. 1.

2. 33 Key Genes Were Identified By Wgcna Analysis

To obtain key genes related to m6A modification, we performed WGCNA on 10,893 genes and finally obtained 18 modules (Fig. 2B-C). Then, the GS value of each module was calculated, with a larger GS indicating that the module was more related to the phenotypic characteristics of the sample (Fig. 2D). We calculated the Pearson correlation coefficient between each module and the phenotypic characteristics of the sample (Fig. 2E).

According to the results, we identified two key modules, greenyellow and brown, and selected the genes in these two modules for subsequent analysis. We constructed a protein-protein interaction (PPI) network based on two modules, and 105 hub genes were screened in the network (Supplementary Fig. 2A). At the same time, 185 hub genes were screened in these two modules according to the thresholds of MM > 0.6 and GS > 0.2 (Supplementary Fig. 2B-C). The intersection contains 33 genes, which are considered to be the key genes related to m6A (Fig. 2F).

3. Construction Of The M6a-gpi In The Tcga Dataset

Univariate Cox analysis was performed on 33 key genes, and the results showed that 10 genes (IL12RB1, IL2RB, IFNG, FASLG, CXCL9, CXCL13, GBP2, CXCL10, CXCR6, and CIITA) had a significant relationship with the prognosis of CRC (P < 0.05). The forest plot is presented (Fig. 3A). To further determine the genes used to construct m6A-GPI, LASSO analysis was performed to identify the 7 most important genes and their coefficients (IL12RB1, FASLG, CXCL13, GBP2, CXCL10, CXCR6, and CIITA) (Fig. 3B-C), we utilized m6A-GPI to determine the patient’s risk score and divided all patients into a high-risk group and a low-risk group based on the median risk score. (Fig. 3D).

Compared with patients with high-risk scores, lower risk scores represent better DFS and a relatively longer survival time in the K-M curves (P < 0.05) (Fig. 3E). At the same time, a ROC curve was used to test the accuracy of m6A-GPI in predicting patient survival. The AUC of the 1-, 3-, and 5-year DFS rates reach 0.66, 0.67 and 0.65, respectively (Fig. 3F), which indicated that m6A-GPI has the potential to predict the DFS of patients in the TCGA cohort.

4. Validation Of The M6a-gpi

Generally, the pathological stage is of great significance to the prognosis of CRC [18], but other factors such as age and gender can also affect the prognosis. Therefore, we tested the stability of m6A-GPI in different clinical characteristics. In the stratified samples, the results showed that the high- and low-risk groups still had significant survival differences after distinguishing age, sex, and stage (P < 0.05) (Fig. 4A-F), which indicates that the m6A-GPI has good stability in stratified samples.

The external dataset was obtained from the GSE17538 cohort, and we used the same formula to calculate the risk score of the patients in this cohort. Similarly, patients were divided into high- and low-risk groups (Fig. 4G). Patients with higher risk scores had poor DFS in the GEO cohort (P < 0.05), which is consistent with the previous analysis of the TCGA cohort (Fig. 4H).

Furthermore, we compared the risk score in the clinical characteristics of the TCGA cohort (age, cancer type, sex, stage and cancer site), and the results showed that the risk score was significantly different in stages and cancer sites (P < 0.05) (Fig. 5A-E). As the stage increases, the risk score has an upward trend, and the risk score of left colon cancer is higher than that of right colon cancer.

5. The Molecular And Mutation Characteristics Of Different M6a-gpi Groups

We compared some potential factors that determine tumor immunogenicity in two subgroups, and the results indicated that the risk score was positively correlated with HRD, CNA, and mRNAsi (Fig. 6). HRD mainly include loss of heterozygosity (LOH), telomere allele imbalance (TAI), and large-scale transition (LST). These three indicators can be used to determine the genomic instability score (GIS) and then evaluate the HRD status. CNA, LOH, TAI, and LST represent the level of chromosome instability [19]. The mRNAsi is an index that can assess the similarity between tumor cells and stem cells and is related to the active biological processes in stem cells and the high degree of tumor dedifferentiation [20].

In addition, we used the “maftools” R package to analyze the distribution of somatic mutations between two subgroups in the TCGA cohort. Then we sorted the genes according to the mutation rate and identified the genes with the highest mutation rate in two groups. The mutation rates of APC, TP53, KRAS, TTN, MUC16, PIK3CA, SYNE1, FAT4, OBSCN, and MUC4 were higher than 20% in both groups (Fig. 7A-B). After we grouped samples according to the risk score, the high-risk samples showed significant amplifications on chromosomes 8, 11, 12, 17, and 20, while deletions were found on chromosomes 1, 3 to 6, 8, 10, 12, 16 to 20. However, the low-risk samples showed significant amplifications on chromosomes 5, 6, 8, 10, 12, 13, 16, 17, 19, and 20, while deletions were found on chromosomes 1, 3 to 8, 10, 15 to 22 (Fig. 7C-D).

6. Immune Characteristics Of Different M6a-gpi Groups

The effect of tumor treatment depends not only on the tumor immunogenicity of the tumor but also on the TIME. TIME is formed by various cells, including immune cells (such as T cells, Tregs, NK cells, and macrophages), endothelial cells, and inflammatory mediators [21]. The role of immune cells is particularly important, and it may affect the patient's response to treatment. To compare the distribution of immune cells in m6A-GPI subgroups, we analyzed the relative proportions of immune cells between the two m6A-GPI subgroups. Compared with low-risk patients, high-risk patients showed more infiltration of NK cells, M0 macrophages, and mast cells. However, in this group, there were fewer CD4 + T cells, CD8 + T cells, M1 macrophages, and M2 macrophages (Fig. 8A).

Then we explored the level of the stromal score, immune score and tumor purity among the two groups, and the results showed that the high-risk group had higher tumor purity and that the low-risk group had a higher stromal score and immune score (Fig. 8B-D). It suggests that the high-risk patients had a higher proportion of cancer cells in the tissue, and the TIME of the low-risk group contains contained abundant immune or matrix components.

7. Independent Prognostic Factor And Nomogram

By using univariate Cox regression analysis and multivariate Cox regression analysis, we sought to determine whether m6A-GPI was an independent prognostic factor for patients with CRC. Univariate Cox analysis showed that m6A-GPI was closely related to the prognosis of CRC [hazard ratio (HR) = 3.041, 90% CI: 2.06–4.5, P < 0.001]. Multivariate Cox analysis further showed that m6A-GPI can be used as an independent predictor in CRC [hazard ratio (HR) = 2.4, 90% CI: 1.6–3.59, P < 0.001] (Fig. 9A). At the same time, we constructed a nomogram and calibration plots of DFS at 1-, 3-, and 5-year in the TCGA dataset (Fig. 9B-E), which provides doctors with a method to quantitatively predict the prognosis of CRC. The accuracy of prediction at 3 years can be increased to 0.75 after combining the risk score and clinical characteristics (Fig. 9F).

8. Validation Of M6a-related Genes Expression Levels In Crc Tissues

To further investigate the expression levels of m6A-related genes in CRC clinical tissues. We first detected the mRNA expression of m6A-related genes (IL12RB1, FASLG, CXCL13, GBP2, CXCL10, CXCR6, and CIITA). Notably, the expression of CIITA was substantially increased in CRC tissues (Fig. 10A). However, there was no significant difference of others genes between CRC and adjacent normal tissues (Fig. 10A). Subsequently, protein level of CIITA was detected in CRC and adjacent normal tissues. Consistent with mRNA levels, the protein level of CIITA was also obviously increased in CRC tissues (Fig. 10B).

Several lines of evidence indicate that m6A modification has become an important target in tumor immunity [22, 23]. We believe that a novel prognostic marker based on m6A-related genes will help predict the prognosis and immune infiltration of colorectal cancer. For this reason, we screened 7 m6A-related genes (IL12RB1, FASLG, CXCL13, GBP2, CXCL10, CXCR6, and CIITA) through WGCNA and LASSO regression in this study and established an m6A-GPI in colorectal cancer. Survival analysis based on m6A-GPI showed that a lower risk score means longer DFS for patients. At the same time, we found that groups with different clinical characteristics (such as tumor site and stage) also showed differences in risk scores. We explored the molecular factors that affect tumor immunogenicity. The risk score was positively correlated with the HRD, CNA, and mRNAsi, and people with higher risk scores may have chromosomal instability. It is worth noting that m6A-GPI not only has good prognostic predictive ability but is also related to tumor immunogenicity, immune cell infiltration in CRC patients, which indicates that m6A-GPI may become a predictive indicator of tumor treatment. We also confirmed that m6A-GPI is an independent prognostic factor for CRC patients, which will provide useful guidance for clinical treatment strategies. Finally, we verified the expression levels of m6A-related genes (IL12RB1, FASLG, CXCL13, GBP2, CXCL10, CXCR6, and CIITA) in CRC tissues. Our results indicated that CIITA might play a crucial role in the prognosis of colorectal cancer.

At present, many studies have shown that m6A-related genes play an indispensable role in cancer progression and metastasis [24], but how they regulate immune cell infiltration in CRC is still unclear. In our study, the tumor immune microenvironment of patients with higher risk scores had increased infiltration of immune cells, such as resting NK cells, and M0 macrophages. Previous studies have shown that macrophages can be recruited to tumor tissues and contribute to tumor angiogenesis [25], which may cause poor DFS in high-risk groups. In addition, a significantly higher proportion of CD4⁺ T cells, CD8⁺ T cells, memory B cells, and M1 macrophages were found in patients in the low-risk group, indicating that there is a greater proportion of T cells and B cells in low-risk CRC tumors. CD8⁺ T cells play a major role in tumor immunity. CD8⁺ T cells differentiate into cytotoxic T cells in the body, and cytotoxic T cells can enter the tumor microenvironment and inhibit the growth of the tumor [26]. In the TIME, there is a tendency for CD8⁺ T cells to increase in METTL3- or METTL14-null tumors, accompanied by increased secretion of IFN-γ, CXCL9 and CXCL10 [13]. Meanwhile, studies have confirmed that the expression of ALKBH5 is specifically upregulated when T cells are activated, and ALKBH5 increases m6A modification on IFN-γ and CXCL2 in CD4⁺ T cells, thereby affecting mRNA stability and protein expression. These modifications lead to changes in the response of CD4⁺ T cells [27]. The infiltration of M1 macrophages can promote inflammation and inhibit tumor cells in the TIME, M1 cells can be activated by IFN-γ and destroy tumors by producing nitric oxide, type 1 cytokines, and chemokines [28], which is consistent with the trend we observed in the low-risk group. Finally, this study found that compared with the high-risk group, CRC patients in the low-risk group had higher immune scores and stromal scores, and had lower tumor purity.

m6A modification changes the TIME, which largely affects the therapeutic response to antitumor immunotherapy. Approximately 85% of CRC patients have mismatch-repair‐proficient or microsatellite instability‐low (pMMR-MSI-L) tumors. This type of patient failed to benefit from any single immunotherapy, but microsatellite instability-high (pMMR-MSI-H) CRC responds well to immunotherapy because it can recruit a large number of immune cells such as CD8⁺ / CD4⁺ T cells and macrophages into the microenvironment [29–32]. Wang et al. proposed that the destruction of m6A methyltransferase enhances the immunotherapy response of pMMR-MSI-L colorectal cancer by regulating the tumor microenvironment and tumor-infiltrating cells [13]. In fact, the loss of METTL3 or METTL14 enhances the interaction between the tumor and the immune system through the IFN-γ-STAT1-IRF1. In another study on “eraser” ALKBH5, researchers found that the knockout of ALKBH5 in mice with CT26 colorectal cancer or B16 melanoma significantly reduced tumor growth and prolonged the survival rate of mice during immunotherapy. This may be related to ALKBH5 inhibiting immune cells in the tumor microenvironment and regulating lactic acid. These processes increase the response to anti-PD-L1 therapy and the loss of ALKBH5 changes the composition of immune cells and metabolite tumor microenvironment [10]. In the subgroups based on m6A-GPI, the TIME and immune cells have changed. We found that tumors in the low-risk group recruited more CD4⁺ / CD8⁺ T cells, our research provides important guidance for predicting the proportion of immune cells in CRC.

Our research provides ideal predictors for the prognosis and immune cell infiltration, but it is undeniable that there are several limitations in this study. First, we used retrospective data from public databases to construct and verify the m6A-GPI, and it would be more rigorous to use a larger-scale prospective data to evaluate its reliability. Second, the population in our study was mainly Americans, and different countries may have deviations in the results due to ethnic differences. In fact, this manuscript is the first part of our research, and our next research project will stem from these results.

In conclusion, m6A-GPI is a promising m6A-related prognostic biomarker. Our study divides patients into different risk subgroups based on m6A-GPI, which will help doctors identify the molecular and immune characteristics and predict the progression, prognosis of CRC. Moreover, m6A-GPI may be a potential indicator in the adjustment of tumor treatment strategies.

Author Contributions

Bin Ma conceived and designed this study, Bin Ma performed the bioinformatic analysis and visualization. Simeng Bao performed the experimental validation. Yongmin Li collected data and performed the statistical analysis. Bin Ma and Simeng Bao wrote the original draft and Bin Ma revised the manuscripts. All authors revised and approved the final manuscript.

Funding

This study was supported by the National Natural Science Foundation of China (81902383); the Revitalizing Liaoning Talents Program (XLYC1907004); the Young and Middle-aged Scientific & Technological Innovation Talent Support Plan of Shenyang City (RC200223); Cultivation Program of National Science Foundation of Liaoning Cancer Hospital (2021-ZLLH-03); Beijing Xisike Clinical Oncology Research Foundation (Y-QL202101-0039).

Acknowledgments

We obtained data from and we sincerely thank The Cancer Genome Atlas (TCGA) database for the establishment and management of the database.

Ethical approval statement and patient consent

Consent for publication

Not applicable.

Competing Interests

The authors have declared that no competing interest exists.

Data Availability Statement

The datasets can be found in https://portal.gdc.cancer.gov/ and https://www.ncbi.nlm.nih.gov/geo/

Thrumurthy, S.G., S.S. Thrumurthy, C.E. Gilbert, P. Ross, and A. Haji, “Colorectal adenocarcinoma: risks, prevention and diagnosis,” Bmj, 2016. 354: p. i3590.
Arnold, M., M.S. Sierra, M. Laversanne, et al., “Global patterns and trends in colorectal cancer incidence and mortality,” Gut, 2017. 66(4): p. 683-691.
Mauri, G., A. Sartore-Bianchi, A.G. Russo, et al., “Early-onset colorectal cancer in young individuals,” Mol Oncol, 2019. 13(2): p. 109-131.
Mao, R., F. Yang, Z. Wang, et al., “Clinical Significance of a Novel Tumor Progression-Associated Immune Signature in Colorectal Adenocarcinoma,” Front Cell Dev Biol, 2021. 9: p. 625212.
Lefèvre, J.H., L. Mineur, M. Cachanado, et al., “Does A Longer Waiting Period After Neoadjuvant Radio-chemotherapy Improve the Oncological Prognosis of Rectal Cancer?: Three Years' Follow-up Results of the Greccar-6 Randomized Multicenter Trial,” Ann Surg, 2019. 270(5): p. 747-754.
Van der Jeught, K., H.C. Xu, Y.J. Li, X.B. Lu, and G. Ji, “Drug resistance and new therapies in colorectal cancer,” World J Gastroenterol, 2018. 24(34): p. 3834-3848.
Lim, H., D. He, Y. Qiu, et al., “Rational discovery of dual-indication multi-target PDE/Kinase inhibitor for precision anti-cancer therapy using structural systems pharmacology,” PLoS Comput Biol, 2019. 15(6): p. e1006619.
Li, X., S. Warren, V. Pelekanou, et al., “Immune profiling of pre- and post-treatment breast cancer tissues from the SWOG S0800 neoadjuvant trial,” J Immunother Cancer, 2019. 7(1): p. 88.
Vassilakopoulou, M., M. Avgeris, V. Velcheti, et al., “Evaluation of PD-L1 Expression and Associated Tumor-Infiltrating Lymphocytes in Laryngeal Squamous Cell Carcinoma,” Clin Cancer Res, 2016. 22(3): p. 704-13.
Li, N., Y. Kang, L. Wang, et al., “ALKBH5 regulates anti-PD-1 therapy response by modulating lactate and suppressive immune cell accumulation in tumor microenvironment,” Proc Natl Acad Sci U S A, 2020. 117(33): p. 20159-20170.
Li, T., P.S. Hu, Z. Zuo, et al., “METTL3 facilitates tumor progression via an m(6)A-IGF2BP2-dependent mechanism in colorectal carcinoma,” Mol Cancer, 2019. 18(1): p. 112.
Han, D., J. Liu, C. Chen, et al., “Anti-tumour immunity controlled through mRNA m(6)A methylation and YTHDF1 in dendritic cells,” Nature, 2019. 566(7743): p. 270-274.
Wang, L., H. Hui, K. Agrawal, et al., “m(6) A RNA methyltransferases METTL3/14 regulate immune responses to anti-PD-1 therapy,” Embo j, 2020. 39(20): p. e104514.
Huang, D., F. Zhang, H. Tao, et al., “Tumor Mutation Burden as a Potential Biomarker for PD-1/PD-L1 Inhibition in Advanced Non-small Cell Lung Cancer,” Target Oncol, 2020. 15(1): p. 93-100.
Newman, A.M., C.L. Liu, M.R. Green, et al., “Robust enumeration of cell subsets from tissue expression profiles,” Nat Methods, 2015. 12(5): p. 453-7.
Jiang, Y.Z., D. Ma, C. Suo, et al., “Genomic and Transcriptomic Landscape of Triple-Negative Breast Cancers: Subtypes and Treatment Strategies,” Cancer Cell, 2019. 35(3): p. 428-440.e5.
Duan, J., Y. Pan, X. Yang, et al., “Screening of T Cell-Related Long Noncoding RNA-MicroRNA-mRNA Regulatory Networks in Non-Small-Cell Lung Cancer,” Biomed Res Int, 2020. 2020: p. 5816763.
Qiang, J.K., R. Sutradhar, V. Giannakeas, et al., “Impact of diabetes on colorectal cancer stage and mortality risk: a population-based cohort study,” Diabetologia, 2020. 63(5): p. 944-953.
Chao, A., C.H. Lai, T.H. Wang, et al., “Genomic scar signatures associated with homologous recombination deficiency predict adverse clinical outcomes in patients with ovarian clear cell carcinoma,” J Mol Med (Berl), 2018. 96(6): p. 527-536.
Zhang, Y., J.T. Tseng, I.C. Lien, et al., “mRNAsi Index: Machine Learning in Mining Lung Adenocarcinoma Stem Cell Biomarkers,” Genes (Basel), 2020. 11(3).
Wang, M., J. Zhao, L. Zhang, et al., “Role of tumor microenvironment in tumorigenesis,” J Cancer, 2017. 8(5): p. 761-773.
Li, M., X. Zha, and S. Wang, “The role of N6-methyladenosine mRNA in the tumor microenvironment,” Biochim Biophys Acta Rev Cancer, 2021. 1875(2): p. 188522.
Gu, Y., X. Wu, J. Zhang, et al., “The evolving landscape of N(6)-methyladenosine modification in the tumor microenvironment,” Mol Ther, 2021. 29(5): p. 1703-1715.
He, L., H. Li, A. Wu, et al., “Functions of N6-methyladenosine and its role in cancer,” Mol Cancer, 2019. 18(1): p. 176.
Fu, L.Q., W.L. Du, M.H. Cai, et al., “The roles of tumor-associated macrophages in tumor angiogenesis and metastasis,” Cell Immunol, 2020. 353: p. 104119.
Iwahori, K., “Cytotoxic CD8(+) Lymphocytes in the Tumor Microenvironment,” Adv Exp Med Biol, 2020. 1224: p. 53-62.
Zhou, J., X. Zhang, J. Hu, et al., “m(6)A demethylase ALKBH5 controls CD4(+) T cell pathogenicity and promotes autoimmunity,” Sci Adv, 2021. 7(25).
Gui, T., A. Shimokado, Y. Sun, T. Akasaka, and Y. Muragaki, “Diverse roles of macrophages in atherosclerosis: from inflammatory biology to biomarker discovery,” Mediators Inflamm, 2012. 2012: p. 693083.
Ganesh, K., Z.K. Stadler, A. Cercek, et al., “Immunotherapy in colorectal cancer: rationale, challenges and potential,” Nat Rev Gastroenterol Hepatol, 2019. 16(6): p. 361-375.
Llosa, N.J., M. Cruise, A. Tam, et al., “The vigorous immune microenvironment of microsatellite instable colon cancer is balanced by multiple counter-inhibitory checkpoints,” Cancer Discov, 2015. 5(1): p. 43-51.
Deschoolmeester, V., M. Baay, F. Lardon, P. Pauwels, and M. Peeters, “Immune Cells in Colorectal Cancer: Prognostic Relevance and Role of MSI,” Cancer Microenviron, 2011. 4(3): p. 377-92.
Le, D.T., J.N. Uram, H. Wang, et al., “PD-1 Blockade in Tumors with Mismatch-Repair Deficiency,” N Engl J Med, 2015. 372(26): p. 2509-20.

No competing interests reported.

Supplementalinformation.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Identification and Validation of m6A-GPI Signatures as a Novel prognostic Model for Colorectal Cancer

Status:

Version 1

Abstract

Figures

Introduction

Methods

1. Datasets Acquisition

2. Identification Of M6a-related Genes

3. Wgcna Analysis

4. Construction And Validation Of M6a-gpi

5. Comprehensive Analysis Of Molecular And Immune Characteristics In Different M6a-gpi Subgroups

6. Independent Prognostic Factor And Nomogram

7. Crc Tissue Samples

8. Quantitative Real-time Rt-pcr

9. Western Blot Analysis

10. Statistical analysis

Results

1. Identification of m6A-related genes

2. 33 Key Genes Were Identified By Wgcna Analysis

3. Construction Of The M6a-gpi In The Tcga Dataset

4. Validation Of The M6a-gpi

5. The Molecular And Mutation Characteristics Of Different M6a-gpi Groups

6. Immune Characteristics Of Different M6a-gpi Groups

7. Independent Prognostic Factor And Nomogram

8. Validation Of M6a-related Genes Expression Levels In Crc Tissues

Discussion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1