m5C regulator-mediated modification patterns and tumor microenvironment infiltration characterization in colorectal cancer: One step closer to precision medicine

Abstract

Background: The RNA modification 5-methylcytosine (m5C) is one of the most prevalent post-transcriptional modifications, with increasing evidence demonstrating its extensive involvement in the tumorigenesis and progression of various cancers. Colorectal cancer (CRC) is the third most common cancer and second leading cause of cancer-related deaths worldwide. However, the role of m5C modulators in shaping tumor microenvironment (TME) heterogeneity and regulating immune cell infiltration in CRC requires further clarification.

Results: The transcriptomic sequencing data of 18 m5C regulators and clinical data of patients with CRC were obtained from The Cancer Genome Atlas (TCGA) and systematically evaluated. We found that 16 m5C regulators were differentially expressed between CRC and normal tissues. Unsupervised cluster analysis was then performed and revealed two distinct m5C modification patterns that yielded different clinical prognoses and biological functions in CRC. We demonstrated that the m5C score constructed from eight m5C-related genes showed excellent prognostic performance, with a subsequent independent analysis confirming its predictive ability in the CRC cohort. Then we developed a nomogram containing five clinical risk factors and the m5C risk score and found that the m5C score exhibited high prognostic prediction accuracy and favorable clinical applicability. Moreover, the CRC patients with low m5C score were characterized by "hot" TME exhibiting increased immune cell infiltration and higher immune checkpoint expression. These characteristics were highlighted as potential identifiers of suitable candidates for anticancer immunotherapy. Although the high m5C score represented the non-inflammatory phenotype, the CRC patients in this group exhibited high level of sensitivity to molecular-targeted therapy.

Conclusion: Our comprehensive analysis indicated that the novel m5C clusters and scoring system accurately reflected the distinct prognostic signature, clinicopathological characteristics, immunological phenotypes, and stratifying therapeutic opportunities of CRC. Our findings, therefore, offer valuable insights into factors that may be targeted in the development of precision medicine-based therapeutic strategies for CRC.

Keywords: 5-methylcytosine; RNA methylation; colorectal cancer; immune infiltrates; precision medicine; tumor microenvironment.

Figure 1

Landscape of m5C methylation regulators. (A) Graphical summary of the study protocol. (B) Overview of the m5C genes and their functions in different cancer types. (C) Pie charts showing the types of m5C regulators (top). The heatmap depicts the mRNA expression levels of 18 m5C regulators between normal mucosae and CRC tissues according to transcriptome data from TCGA and GTEx (bottom).

Figure 2

Landscape of the genetic alterations and transcriptional variations in the m5C genes in CRC. (A) Boxplot comparison of the differential expression levels of m5C genes between tumor and normal tissues from the TCGA-CRC dataset. (B) Circos plot showing the chromosomal distribution of 18 m5C genes. (C) PCA of the expression profiles of 18 m5C regulators. (D) CNV frequencies of the 18 m5C regulators. Column height represents the change in frequency. (E) Spearman’s correlation analysis of the 18 m5C genes from the TCGA-CRC dataset. (F) Mutation frequencies of the 18 m5C genes from the TCGA-CRC cohort. *P < 0.05, ***P < 0.001, ns, no significant.

Figure 3

Different m5C modification patterns showing distinct biological characteristics. (A) Correlations and correlation coefficients between the 18 m5C regulators in CRC. Each circle represents an individual gene, and the size of the circle represents the associated prognosis. Data were generated using the log-rank test (range: 0.1–0.0001). The green or purple dots represent favorable factors or risk factors for OS, respectively, and red or blue lines indicate positive or negative correlations between the regulators, respectively. (B) Consensus clustering matrix (k = 2). (C) Survival analysis of the patients in the clusters generated according to m5C scores from the TCGA dataset. (D) Heatmap generated using DEGs between m5C clusters A and B. (E) Heatmap showing the GSVA analysis, which showed the activation or inhibition of biological pathways according to the m5C clusters. (F) GO and (G) KEGG analyses of the DEGs between m5C clusters A and B. (H) Heatmap showing the immunotherapy-predicted pathways between m5C clusters A and B.

Figure 4

Construction of a prognostic signature using m5C-related genes. (A) Clustering dendrogram of the CRC samples and associated clinical traits. (B) Clustering dendrogram of the genes with dissimilarity on the basis of topological overlap with the corresponding module colors. (C) Heatmap of the association between module eigengenes and clinical phenotypes of CRC. (D) LASSO coefficient profiles of m5C-related genes. (E) Partial likelihood deviance for the LASSO coefficient profiles. (F) Survival analysis of the CRC patients stratified by the m5C risk score. (G) ROC curves for predicting the sensitivity and specificity of 1-, 3-, and 5-year OS based on the m5C score. (H) Alluvial diagram of subtype distributions in the groups with different m5C scores and survival outcomes. (I, J) Ranked dot and scatter plots showing the m5C score distribution and patient survival statuses. (K) Heatmap of the expression of eight m5C-related genes in the different m5C risk groups. (L, M) Univariate and multivariate Cox analyses of the m5C risk scores and clinical variables.

Figure 5

Relationship between the m5C score, clinical features, and immunological characteristics. (A) Heatmap of the distribution of clinical characteristics and corresponding m5C risk score in each CRC sample. (B–E) Heatmap and table indicating the distribution of the clinicopathological features between the high- and low-m5C-score groups. (F) Relationships between the m5C risk score and CSC index. (G) Spearman’s correlation analysis of the m5C scores and TMB. (H) TMB in the different m5C risk groups. (I, J) Relationships between the m5C risk score and MSI. *P < 0.05, ***P < 0.001.

Figure 6

Correlation between the m5C risk score and immune phenotypes. (A) Heatmap showing the associations between the m5C score and the enrichment scores of several therapeutic signatures. (B) Differences in the expression levels of the immune effector genes between the two m5C score groups. (C) Spearman’s correlation analysis of the m5C score with the activities of cancer immunity cycles (left) and immune-related pathways analyzed by the ssGSEA (right). (D) Differences in the expression levels of immune checkpoint genes between the two m5C score groups. (E) Heatmap showing the significant differential expression of immunomodulators between the two risk groups. *P < 0.05, **P < 0.01, ***P < 0.001.

Figure 7

Mutation profiles and drug-susceptibility analysis. (A) Correlations between the TIDE scores and clinicopathological features (survival status, TNM stage, gender, age, and m5C score subtypes). (B) Comparison of the TIDE score between the two m5C score subgroups. (C) Survival analysis of the CRC patients in the high- and low-TIDE-score groups. (D, E) Waterfall chart depicting the somatic mutation landscapes in the low- and high-m5C-score groups. (F) Relationship between the m5C score and chemotherapeutic sensitivity. (G) Association between the m5C score and targeted treatment sensitivity. ***P < 0.001.

Figure 8

Construction and validation of the m5C score-based nomogram. (A) Development of the nomogram for predicting the 1-, 3‐ or 5‐year OS of CRC patients. (B) ROC curves for the nomogram for predicting the 1-, 3-, and 5-year OS. (C–E) Calibration plots of the nomogram for predicting the 1-, 3- and 5-year OS. (F) DCA for the nomogram assessing clinical utility. (G) Kaplan–Meier survival curves on the basis of the m5C score calculated using the nomogram. *P < 0.05, **P < 0.01, ***P < 0.001.

Figure 9

External validation of the m5C score using a GEO CRC dataset (GSE38832) and an independent CRC cohort. (A) Kaplan–Meier survival curve for patients with low and high m5C scores. (B) Heatmap of the associations between the m5C score and the enrichment scores of several therapeutic signatures in the GSE38832 dataset. (C) Spearman’s correlation analysis of the m5C score with activities of cancer immunity cycles (left) and immune-related pathways using the GSE38832 dataset. (D) Heatmap of the eight m5C-related risk gene profiles in 40 pairs of CRC tumor tissues and adjacent normal tissues. (E) Histogram showing the ratio of the AJCC, T, and N stages between the low- and high-risk groups. (F) Relationships between the m5C risk scores and clinicopathological characteristics of 40 CRC patients in our cohort. (G) The relative mRNA expression levels of several immune checkpoint genes were examined using RT-qPCR. (H, I) PD-L1 and CTLA4 expression was detected using immunofluorescence between the CRC patient samples in the low- (left) and high- (right) m5C-score groups. *P< 0.05, **P< 0.01, ns, no significant.