m5C related-regulator-mediated methylation modification patterns and prognostic significance in breast cancer

Abstract

5-Methylcytosine (m⁵C) is closely associated with cancer. However, the role of m⁵C in breast cancer(BC)remains unclear. This study combined single-cell RNA sequencing (scRNA-Seq) and transcriptomics datasets to screen m⁵C regulators associated with BC progression and analyze their clinical values. Firstly, This study elucidates the mechanisms of the m⁵C landscape and the specific roles of m⁵C regulators in BC patients. we found that the dysregulation of m⁵C regulators with m⁵Cscore play the essential role of the carcinogenesis and progression in epithelial cells and myeloid cells of BC at single cell level. External validation was conducted using an independent scRNA-Seq datasets. Then, three distinct m⁵C modification patterns were identified by transcriptomics datasets. Based on the m⁵C differentially expressed regulators, the m⁵Cscore was constructed, and used to divide patients with BC into high and low m⁵Cscore groups. Patients with a high m⁵Cscore had more abundant immune cell infiltration, stronger antitumor immunity, and better prognoses. Finally, Quantitative real-time (PCR) and immunohistochemistry were used for the in vitro experimental validation, which had extensive prognostic value. In this study, we aimed to assess the expression of m⁵C regulators involved in BC and investigate their correlation with the tumor microenvironment, clinicopathological characteristics, and prognosis of BC. The m⁵C regulators could be used to effectively assess the cell specific regulation prognosis of patients with BC and develop more effective immunotherapy strategies.

Keywords: Breast cancer; Cell specific regulation; Epigenetics; Immunotherapy; m5C RNA modification.

Fig. 1

Single-cell transcriptomic landscape of m⁵C regulators regulating breast tissue key pathways in breast cancer (BC). (a) T-Distributed Stochastic Neighbor Embedding (TSNE) plot of normal cells and BC cells, colored by cell type. (b) Heatmap for differences in the expression of m⁵C regulators in different cell types between normal and BC samples. Red, up-regulation; blue, down-regulation. (c) Heatmap of each m⁵C regulators expression in each cell type. (d) Differential expression of genes (DEGs) in different cell types of BC patients compared with control samples. Red, up-regulation; blue, down-regulation. (e) Dot plot showing DEGs expression patterns of m⁵C regulators of each cell types. Each dot represents a regulator, of which the color saturation indicates the average expression level, and the size indicates the percentage of cells expressing the regulator. (f) Kaplan-Meier survival analysis based on m⁵C regulator expression. Red, high expression of m⁵C regulator; blue, low expression of m⁵C regulator. (g) Analysis of m⁵Cscore for 5 cell types. The two T-SNE plots of m⁵C regulators expression in BC samples. (h) ALYREF. (i) DNMT1.

Fig. 2

Regulation of tumor-related pathways by the m5C regulators. (a) Correlation analysis was used to analyze the association between the m⁵C regulators and tumor related pathways. Red, positive correlation; blue, negative correlation. (b–d) Differences in m⁵C signature between tumor and normal groups in different cell types, including endothelial cells (En), epithelial cells (Ep), and myeloid cells (mye). (e) Networks of WGCNA module which included ALYREF in myeloid cells. (f) Functional enrichment of module which included ALYREF in myeloid cells. (g) Networks of WGCNA module which included DNMT1 in myeloid cells. (h) Functional enrichment of module which included DNMT1 in myeloid cells. (i) UMAP plots showing the expression of different cell types in other BC tissues. (j) Proportion of m⁵C regulators expression in different cell types. (k) Analysis of m⁵C score for several cell types. (l) Networks of WGCNA module which included ALYREF in myeloid cells. (m) Functional enrichment of module which included ALYREF in myeloid cells. (n) Networks of WGCNA module which included DNMT1 in myeloid cells. (o) Functional enrichment of module which included DNMT1 in myeloid cells.

Fig. 3

Evaluation of m⁵C methylation modification patterns. (a) Differential expression of m⁵C regulators between breast cancer and normal breast tissues. Blue, normal tissue; red, tumor tissue. The lines in the boxes represent the median value, the bottoms and tops of the boxes represent the interquartile range, and the dots represent outliers. ***P < 0.001, **P < 0.01, *P < 0.05. Differences among the three modification patterns were tested by one-way ANOVA. (b) The interaction between m⁵C regulators in breast cancer. The lines connecting the m⁵C regulators represent the interaction between them. Blue, negative correlation; red, positive correlation. (c) Three different m⁵C modification subtypes were identified by unsupervised clustering based on m⁵C regulators (m⁵C cluster A, B, and C). (d) PCA derived from the m⁵C clusters showed a difference between the three clusters. Blue, m⁵C gene cluster A; yellow, m⁵C gene cluster B; and red, m⁵C gene cluster C. (e) Survival analysis based on the three m⁵C clusters in 1089 patients with breast cancer in the TCGA-BRCA cohort (P = 0.015, log-rank test). Blue, 567 patients in m⁵C cluster A; yellow, 334 patients in m⁵C cluster B; and red, 188 patients in m⁵C cluster C. (f) Unsupervised clustering of 17 m⁵C regulators in the TCGA-BRCA cohort identified a significant difference in the expression of regulators among the three modification patterns. The m⁵C clusters, TCGA project, age, sex, TNM classification, clinical stage, and survival status were used as patient annotations. Red, high expression of regulators; blue, low expression of regulators. (g–i) GSVA enrichment analysis showing the activation states of biological pathways in distinct m⁵C modification patterns. Red, activated pathways; blue, inhibited pathways. (g) m⁵C cluster A compared with m⁵C cluster B; (h) m⁵C cluster A compared with m⁵C cluster C; (i) m⁵C cluster B compared with m⁵C cluster C.

Fig. 4

Construction of m⁵C gene signatures and functional annotation. (a) Overlapping m⁵C phenotype-related DEGs in the three m⁵C clusters. (b) Three different genomic subtypes identified by unsupervised clustering based on the overlapping m⁵C phenotype-related DEGs. (c) PCA derived from the three m⁵C gene clusters showed a difference between gene clusters. Blue, m⁵C gene cluster A; yellow, m⁵C gene cluster B; and red, m⁵C gene cluster C. (d) Survival analysis based on the three m⁵C gene clusters in 1089 patients from the TCGA-BRCA cohort (P = 0.025, log-rank test). Blue, 326 patients with m⁵C gene cluster A; yellow, 621 patients with m⁵C gene cluster B; and red, 142 patients with m⁵C gene cluster C. (e) GO functional enrichment analysis of 2312 overlapping m⁵C phenotype-related DEGs. (f) KEGG pathway enrichment analysis of 2312 overlapping m⁵C phenotype-related DEGs. (g) Differences in the m⁵Cscore among the three m⁵C clusters in the TCGA-BRCA cohort (P < 0.001, Kruskal–Wallis test). Blue, m⁵C cluster A; yellow, m⁵C cluster B; and red, m⁵C cluster C. (h) Differences in the m⁵Cscore among the three m⁵C gene clusters in the TCGA-BRCA cohort (P < 0.001, Kruskal–Wallis test). Blue, m⁵C gene cluster A; yellow, m⁵C gene cluster B; and red, m⁵C gene cluster C.

Fig. 5

Association between the m⁵Cscore and response to pharmacotherapy. (a) The relative distribution of TIDE scores was compared between the low and high m⁵Cscore groups. The lines in the boxes represent the median value, the bottoms and tops of the boxes represent the interquartile range, and the dots represent outliers. Blue, low m⁵Cscore group; red, high m⁵Cscore group. (b) Survival analysis for 858 patients with a high TIDE score and 231 with a low TIDE score (P < 0.05, log-rank test). (c) Survival curves of TIDE scores combined with m⁵C scores (P < 0.05, log-rank test). (d) Comparisons of the proportions of non-responders and responders to ICB between the low and high m⁵Cscore groups; the ROC curves of the TIDE score model in patients with BC (AUC: 0.846, 95% CI: 0.801–0.888). Blue, non-responder groups; red, responder groups. (e) Relative distribution of immune dysfunction scores between the low and high m⁵Cscore groups (P < 0.001). Blue, low m⁵Cscore group; red, high m⁵Cscore group. (f) Survival analysis stratified by both m⁵Cscore and immune dysfunction scores (P = 0.001, log-rank test). (g) survival analysis for 977 patients with a high immune exclusion and 112 with a low immune exclusion (P < 0.05, log-rank test). (h) Survival analyses stratified by both the m⁵Cscore and immune exclusion (P < 0.001, log-rank test). (i–t) Predicted response of patients to six chemotherapeutic drugs based on the m⁵Cscore. (I, j) bortezomib; (k, l) erlotinib; (m, n) roscovitine; (o, p) salubrinal; (q, r) sorafenib; (s, t) vinorelbine. Blue, low m⁵Cscore group; red, high m⁵Cscore group.

Fig. 6

Correlation between clinicopathological characteristics and the m⁵Cscore. (a) m⁵Cscore based on survival status (P < 0.001). The lines in the boxes represent the median value, the bottoms and tops of the boxes represent the interquartile range, and the dots represent outliers. Blue, living patients; red, deceased patients. (b) The proportions of living and dead patients with BC in the low and high m⁵Cscore groups. In the low m⁵Cscore group, 82% of patients were alive and 18% were dead, and in the high m⁵Cscore group, 90% of patients were alive and 10% were dead. Blue, living patients; red, deceased patients. (c) The m⁵Cscore based on N stage. The Kruskal–Wallis test was used to compare the statistical difference between five N stage groups. The lines in the boxes represent the median value, the bottoms and tops of the boxes represent the interquartile range, and the dots represent outliers. Blue, N₀ stage group; red, N₁ stage group; yellow, N₂ stage group; purple, N₃ stage group; green, N_x stage group. (d) HER2 expression status in the low and high m⁵Cscore groups. In the low m⁵Cscore group, 81% and 19% of patients had HER2 + and HER2- disease, respectively; in the high m⁵Cscore group, 87% and 13% of patients had HER2+ and HER2− disease, respectively. Blue, patients with HER2− disease; red, patients with HER2+ disease. (e) Sankey diagram showing the flow of m⁵C cluster, m⁵C gene cluster, m⁵Cscore, and survival status. (f) Sankey diagram showing the flow of m⁵C cluster, m⁵C gene cluster, m⁵Cscore, and age. (g–l) Kaplan–Meier survival analysis based on the m⁵Cscore in subgroups with different clinical characteristics. Red, high m⁵Cscore group; blue, low m⁵Cscore group. (g) patients ≤ 65 years; (h) patients with T₂ stage; (i) patients with stage III disease; (j) patients with N₁ disease; (k) patients with M₀ disease; (l) Women.

Fig. 7

The expression of m5C-related genes was verified using RT-qPCR and IHC. (a–d) Differential expression of m⁵C-related genes in normal and tumor groups. (e–m) IHC of the ALYREF, DNMT1 and DNMT3a in tumor tissue. (n–p) Percentage of positive staining for m⁵C-related genes between the paracancerous and tumor groups.