Classification and immune invasion analysis of breast cancer based on m6A genes
Introduction
Breast cancer (BRCA) is one of the most common malignant tumors in women. According to the latest global cancer data, there were 2.26 million new cases of BRCA and 680,000 deaths worldwide in 2020. The rapid growth of new cases of BRCA has replaced lung cancer as the most prolific cancer in the world (1), which seriously threatens the physical and mental health of women. Thanks to advances in local and systemic therapy, about 70–80% of patients with early metastasis-free BRCA can be cured, but advanced-stage patients remain unable to achieve remission with existing treatments (2). At the molecular level, BRCA is a highly heterogeneous disease (3), and multiple studies have shown that genetics, epigenetics, and transcriptome modifications influence the occurrence and progression of BRCA (4-6).
Methylation of N6 adenosine (m6A) is the most common and abundant methylation modification in eukaryotic messenger RNAs (mRNAs). Writers, erasers, and readers are the 3 important components to identify m6A modification sites and regulate the splicing, transport, translation, and stability of RNA. Modification of m6A can have positive or negative effects on the occurrence, differentiation, invasion, and metastasis of BRCA (6,7). In addition, there is growing evidence that m6A modification can control the immune response and plays an important role in both congenital and acquired immune responses in tumors (8-10). Study of the m6A protein and its inhibitors can reveal new opportunities for early diagnosis and treatment of BRCA, especially in combination with immunotherapy.
Based on the heterogeneity of BRCA, researchers have found significant differences in proliferation, apoptosis, invasion, and migration of tumor cells with different BRCA subtypes (11), as well as significant differences in response to treatment and prognosis. Therefore, continuing research on molecular typing of BRCA can assist in establishing the basis for precise and individualized treatment plans. To date, there have been few studies on molecular typing of the m6A regulatory gene and BRCA. Our study aimed to explore the typing of BRCA based on the m6A regulatory gene, evaluate the differences of immune cell infiltration, and construct a scoring system, to provide a theoretical basis for basic research and clinical treatment. We present the following article in accordance with the REMARK reporting checklist (available at https://dx.doi.org/10.21037/atm-21-3404).
Methods
Data set acquisition and data preprocessing
We obtained tumor mutation burden (TMB) data from 1,109 tumors and 113 normal samples in the breast cancer cohort of the Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov/). We downloaded copy number variation (CNV) data for TCGA queue from the University of California Santa Cruz (UCSC). In addition, we searched the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) to download GSE20685 breast cancer data sets (including tumor, 327 cases). All RNA-seq data was normalized using the ComBat function in the sva package (https://bioconductor.org/packages/sva/). The CNV data, expression differences, mutations, gene types, and correlations among 23 breast cancer-related m6A methylation genes were analyzed by the R software package (https://www.R-project.org/).
Survival analysis
Univariate Cox regression analysis and the Kaplan-Meier method were used to evaluate the prognostic role of 23 genes associated with m6A methylation in breast cancer. When the Kaplan-Meier method was used to plot survival curves, a log-rank P<0.05 was considered statistically significant.
Consistency cluster analysis and immunoinfiltration analysis
We performed Cluster typing using the Consensus Cluster Plus software package (https://bioconductor.org/packages/ConsensusClusterPlus/) based on the common expression of 23 m6A methylation regulatory genes. At the same time, the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm was used to calculate the immune cell content of tumor samples (28 immune cells) and compare the differences in the immune cell content between different m6A types. Then, we analyzed the gene differences between the groups according to the m6a typing (RNAseq expression matrix of TCGA and GEO datasets) and obtained differentially expressed genes (DEGs). We then performed univariate Cox analysis with a P<0.05 for filtering and obtained the prognostic DEGs. Based on the prognostic DEGs, consistency clustering was repeated to compare the m6A typing of the 2 different algorithms.
Construction of m6A score and its clinical significance
Principal component analysis (PCA) was performed on prognostic DEGs, major components 1 and 2 were extracted, and a new score (M6AScore = PCA1 + PCA2) was set to be obtained, which was named the m6AScore. The correlation between m6AScore and immune cells, TMB, different clinical subtypes, and PD-L1 was analyzed.
Statistical analysis
The software package of R language (https://www.r-project.org/) was used to analyze the data, the ComBat function in the sva package was used to normalize all RNA-seq data, and the Kaplan-Meier method was used to draw the survival curve. Logarithmic rank P<0.05 was considered statistically significant. The Consensus Cluster Plus software package was used for cluster typing and ssGSEA algorithm was used to evaluate tumor-infiltrating immune cells.
Ethical statement
The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).
Results
Genetic variation and prognosis of the m6A gene in BRCA
In this study, we first analyzed CNV mutations of 23 genes in BRCA. Among them, VIRMA and YTHDF1 mainly showed amplification, while WTAP, RBM15, ZC3H13, YTHDF2, and RBM15b had widespread deletion. In addition, the expression of these genes in tumors and normal tissues were analyzed. High expression of VIRMA, RBM15, YTHDF1, YTHDF2, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, and IGFBP2 was detected in tumors. The expressions of METL14, METL16, ZC3H13, RBM15b, YTHDC1, IGFBP1, IGFBP3, and FTO were significantly decreased in tumors (Figure 1). Patients with high expression of IGFBP1, IGFBP3, and RBM15b had a better prognosis, while those with low expression of ALKBH5, FMR, VIRMA, WTAP, YTHDF1, and YTHDF3 had a worse prognosis (Figure 2). We then analyzed the incidence of somatic mutations. A total of 57 (5.78%) of the 986 samples had genetic changes in the m6A gene, including missense mutations, nonsense mutations, and frameshift deletion mutations. The gene ZC3H13 had the highest mutation frequency, followed by LRPPRC and RBM15 (Figure 1). We analyzed the role of 23 m6A methylation genes in BRCA (“writers”: METTL3, METTL14, METTL16, WTAP, ZC3H13, RBM15, and RBM15b; “readers”: YTHDC1, YTHDC2, YTHDF1, YTHDF3, YTHDF2, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP2, IGFBP3, and RBMX; and “erasers”: ALKBH5 and FTO), which suggested that the interaction between the genes HNRNPC and hnRNPA2B1 might be risk factors, while IGFBP3 and YTHDF3 might be potential protective factors, and were associated with tumor progression and prognosis (Figure 1). Based on the above analysis, there are significant differences and associations between the m6A gene and BRCA genome and transcription, and the expression changes and genetic variation of the m6A gene play an important role in regulating the occurrence and development of BRCA and are closely related to prognosis.
Genotyping based on m6A methylation regulatory genes
Based on 23 m6A methylation regulatory genes, TCGA queue and BRCA samples from the GSE20685 dataset were grouped into K subtypes using the partitioning around medoids (PAM) algorithm in Consensus Cluster Plus. According to the cumulative distribution function (CDF) curve of the consensus score and the area under the curve (AUC), when k=3, 3 m6A subtypes (Figure 3A-3D) were distinguished. To further understand whether these 3 m6A subtypes can reflect different clinical outcomes of patients. Based on overall survival (OS) analysis, there was no significant difference in prognosis among the 3 subtypes (P=0.439) (Figure 3E). Surprisingly, when we tried to use the ssGSEA algorithm to calculate the content of immune cells in tumor samples (23 types of immune cells) and compared the infiltration of immune cells among different m6A subtypes, we found that 3 subtypes were filled with a large number of immune cells, and the infiltration level of most subtypes was high (Figure 3F). Then, we conducted paired genetic difference analysis for the 3 subtypes and obtained 3,407 DEGs (Figure 3G). Univariate Cox analysis and P value were used for screening, 276 prognostic DEGs were obtained, and then classified again. This revealed a further division of m6A into 3 types, and there was a significant difference in prognosis between each type (P<0.001) (Figure 3H-3K). Therefore, based on the prognostic DEGs of 23 m6A methylated genes, BRCA can be divided into 3 molecular subtypes (A, n=452; B, n=588; C, n=377), all of which were closely related to immune cell infiltration, and differences in expression were detected.
Construction of m6A score
Cluster analysis was performed on 276 genes to further evaluate the correlation between prognostic DEGs genes and clinical significance. We found that older patients with late clinical-stage and poor clinical prognosis were mainly molecular subtype C (Figure 4A). Although our study found that m6A plays an important role in immune infiltration and prognosis of BRCA, these analyses were based on patient populations only and were not able to predict the pattern of action of m6A in individual patients. Therefore, PCA analysis was performed on 276 differential genes to develop an m6A score (m6AScore = PCA1 + PCA2) to further identify the relevant clinical features of m6A. The results showed that the m6A score was closely related to A, B, and C 3 subtypes, and each gene cluster was closely related. Patients with a high m6A score had a better prognosis (P=0.002), and patients with poor prognosis were mainly concentrated in type C (Figure 4B-4E).
Clinical significance of m6A score
We performed validation the m6A score to clarify its potential prognostic value. Correlation analysis of immune cells showed that central memory CD8 T cells and neutrophils were positively correlated with high m6A score, while most of the remaining immune cells were negatively correlated with low m6A score (Figure 5A). In terms of TMB, a high m6A score correlated with higher TMB, and was concentrated in the B and C types. Patients with high TMB had worse OS (P=0.001), and high m6A score with high TMB had a poor prognosis, consistent with the expected results (Figure 5B-5E). In terms of age, a higher m6A score was associated with a worse prognosis regardless of age ≥45 or age <45 years (Figure 5F-5G). In terms of tumor, node, metastasis (TNM) staging, in addition to patients with early metastases, patients with high m6A scores had a better prognosis (Figure 5H-5M). A Low m6A score was associated with low PD-L1 expression, and was concentrated in patients with TNBC (Figure 5N-5P), which may be the reason for the poor prognosis in the group with a low m6A score. In conclusion, our study found that the m6A score was related to immune cell infiltration and response to immunotherapy, and had a certain prognostic predictive ability.
Discussion
Based on TCGA cohort and GEO dataset, the genetic mutation of the m6A regulatory gene in BRCA was firstly analyzed. The results showed that BRCA exhibited unique somatic mutations and genetic disorders due to its heterogeneity. There were deletion and amplification of m6A regulatory genes in CNV to varying degrees, but the mutation frequency was not high, and the mutation rate of somatic cells was only 5.78%, which mainly presented as meaningless mutation. A study of the m6A regulatory gene spanning 33 cancers showed that the frequency of mutations in the m6A regulatory gene was very low in tumors, ranging from 0.02% to 8.07% (12), which is consistent with our findings. The IGFBP3 and YTHDF3 genes may be protective factors, while HNRNPC and HNRNPA2B1 may be risk factors, further affecting the gene expression between BRCA and normal tissues. Surprisingly, several studies have shown that HNRNPC is overexpressed in BRCA and cell lines and that knocking down HNRNPA2B can reduce the proliferation of breast cancer cells, induce apoptosis, and prolong the S phase of the cell cycle in vitro (13). Subsequently, Wu et al. further demonstrated that inhibition of HNRNPC in breast cancer cells (MCF7 and T47D) inhibited cell proliferation and tumor growth (14). In addition, YTHDF3 overexpression predicted a poor prognosis for OS, and YTHDF3 upregulation was an independent prognostic factor for OS in BRCA patients. Externally, YTHDF3 overexpression was demonstrated to be present in BRCA patients with brain metastasis (15,16). The loss of glucose-regulated protein 78 can reduce the entry of IGFBP-3 into cells, thereby altering its tumor-promoting effect and predicting a poor prognosis in BRCA patients (17), which is quite different from our results. In general, mutations may have little influence on m6A modification, but they play an indispensable role in the genetic or epigenetic factors of BRCA, and are closely related to prognosis.
There is growing evidence that the m6A regulatory gene plays a role in inflammation, immunity, and tumors through regulators (18-20); however, no studies have been conducted to explore the relationship between m6A regulatory genes and BRCA molecular typing. Therefore, identification of the corresponding molecular subtypes in BRCA based on m6A regulatory genes may facilitate understanding of the mechanism of m6A modification of BRCA and its potential clinical applications. In this study, we identified 3 different BRCA subtypes, which are characterized by a large number of immune cell infiltrates and complex relationships between gene clusters within each subtype. Previous studies have shown that interactions between different m6A regulatory genes influence tumor differentiation and pathogenesis (21), which explains the different gene expression patterns among molecular subtypes. It was found that IGF2BP1, IGF2BP2, and IGF2BP3 were upregulated in TNBC and downregulated in other types of BRCA and that CBLL1 was upregulated in estrogen receptor (ER)-positive breast cancer (22,23). In addition, there is growing evidence that the tumor microenvironment plays an important role in tumor progression and immunotherapy (24). Han et al. demonstrated that binding of the m6A reader YTHDF1 promotes the translation of encoded lysosomal proteases, thereby reducing the cross-presentation of dendritic cells (DCs) to tumor antigens, and that inhibition of YTHDF1 reduces immune escape and thus improves the therapeutic effect of PD-1 (25). Based on these findings, we hypothesize that the m6A regulatory gene plays an important role in immune activation and tumor cell infiltration and that some patients may benefit from immunotherapy.
To better quantify the molecular subtype of m6A regulatory genes and define different m6A modification methods, we established the m6A score. Neutrophils and central memory CD8 T cells showed a high m6Ascore for infiltrative characteristics, while the rest of the immune cells showed a low m6A score. Neutrophils play a key role in tumors. Neutrophils have a direct role in killing tumors. Neutrophils combined with lymphocytes can predict the prognosis of BRCA patients, and targeting neutrophils is a potential treatment approach for BRCA (26,27). Further analysis showed that a higher m6A score was positively associated with mutations in TMB and had a poorer prognosis. Immunotherapy based on immune checkpoint inhibitors has some efficacy in TNBC, and these patients often show higher TMB and immune cell infiltration (28-30). In addition, we observed a strong correlation between the m6A score and PD-L1, a predictor of the immune response, suggesting that the m6A gene may influence the efficacy of immunotherapy and that this score could be used in clinical practice to assist in the formulation of treatment plans. There were some limitations in our study, and more clinical samples, especially those from eastern populations, are needed for verification of our findings. Conversely, more clinicopathological features should be included to improve the accuracy of the predictive model.
Conclusions
In this study, we comprehensively evaluated BRCA samples based on 23 m6A regulatory genes and performed subtype analysis to construct a scoring system. At the same time, we evaluated the immune infiltration, which helped us to improve the understanding of the role of m6A in BRCA and may provide a theoretical basis for immunotherapy.
Acknowledgments
Funding: None.
Footnote
Reporting Checklist: The authors have completed the REMARK reporting checklist. Available at https://dx.doi.org/10.21037/atm-21-3404
Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://dx.doi.org/10.21037/atm-21-3404). The authors have no conflicts of interest to declare.
Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).
Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.
References
- 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:209-49. [Crossref] [PubMed]
- Harbeck N, Penault-Llorca F, Cortes J, et al. Breast cancer. Nat Rev Dis Primers 2019;5:66. [Crossref] [PubMed]
- Larsen MJ, Thomassen M, Gerdes AM, et al. Hereditary breast cancer: clinical, pathological and molecular characteristics. Breast Cancer (Auckl) 2014;8:145-55. [Crossref] [PubMed]
- Bianchini G, Balko JM, Mayer IA, et al. Triple-negative breast cancer: challenges and opportunities of a heterogeneous disease. Nat Rev Clin Oncol 2016;13:674-90. [Crossref] [PubMed]
- Wang FW, Ao X, Fu SM. Expression of SOX11 and HER2 and their association with recurrent breast cancer. Transl Cancer Res 2019;8:248-54. [Crossref]
- Chen Y, Lin Y, Shu Y, et al. Interaction between N6-methyladenosine (m6A) modification and noncoding RNAs in cancer. Mol Cancer 2020;19:94. [Crossref] [PubMed]
- Chen B, Li Y, Song R, et al. Functions of RNA N6-methyladenosine modification in cancer progression. Mol Biol Rep 2019;46:1383-91. [Crossref] [PubMed]
- Geng Y, Guan R, Hong W, et al. Identification of m6A-related genes and m6A RNA methylation regulators in pancreatic cancer and their association with survival. Ann Transl Med 2020;8:387. [Crossref] [PubMed]
- Zhang C, Fu J, Zhou Y. A Review in Research Progress Concerning m6A Methylation and Immunoregulation. Front Immunol 2019;10:922. [Crossref] [PubMed]
- Paramasivam A, Vijayashree Priyadharsini J. Novel insights into m6A modification in circular RNA and implications for immunity. Cell Mol Immunol 2020;17:668-9. [Crossref] [PubMed]
- Perou CM, Sørlie T, Eisen MB, et al. Molecular portraits of human breast tumours. Nature 2000;406:747-52. [Crossref] [PubMed]
- Li Y, Xiao J, Bai J, et al. Molecular characterization and clinical relevance of m6A regulators across 33 cancer types. Mol Cancer 2019;18:137. [Crossref] [PubMed]
- Hu Y, Sun Z, Deng J, et al. Splicing factor hnRNPA2B1 contributes to tumorigenic potential of breast cancer cells through STAT3 and ERK1/2 signaling pathway. Tumour Biol 2017;39:1010428317694318 [Crossref] [PubMed]
- Wu Y, Zhao W, Liu Y, et al. Function of HNRNPC in breast cancer cells by controlling the dsRNA-induced interferon response. EMBO J 2018;37:e99017 [Crossref] [PubMed]
- Liu L, Liu X, Dong Z, et al. N6-methyladenosine-related Genomic Targets are Altered in Breast Cancer Tissue and Associated with Poor Survival. J Cancer 2019;10:5447-59. [Crossref] [PubMed]
- Chang G, Shi L, Ye Y, et al. YTHDF3 Induces the Translation of m6A-Enriched Gene Transcripts to Promote Breast Cancer Brain Metastasis. Cancer Cell 2020;38:857-871.e7. [Crossref] [PubMed]
- Zielinska HA, Daly CS, Alghamdi A, et al. Interaction between GRP78 and IGFBP-3 Affects Tumourigenesis and Prognosis in Breast Cancer Patients. Cancers (Basel) 2020;12:3821. [Crossref] [PubMed]
- Chen XY, Zhang J, Zhu JS. The role of m6A RNA methylation in human cancer. Mol Cancer 2019;18:103. [Crossref] [PubMed]
- Jian D, Wang Y, Jian L, et al. METTL14 aggravates endothelial inflammation and atherosclerosis by increasing FOXO1 N6-methyladeosine modifications. Theranostics 2020;10:8939-56. [Crossref] [PubMed]
- Shulman Z, Stern-Ginossar N. The RNA modification N6-methyladenosine as a novel regulator of the immune system. Nat Immunol 2020;21:501-12. [Crossref] [PubMed]
- He L, Li H, Wu A, et al. Functions of N6-methyladenosine and its role in cancer. Mol Cancer 2019;18:176. [Crossref] [PubMed]
- Makdissi FB, Machado LV, Oliveira AG, et al. Expression of E-cadherin, Snail and Hakai in epithelial cells isolated from the primary tumor and from peritumoral tissue of invasive ductal breast carcinomas. Braz J Med Biol Res 2009;42:1128-37. [Crossref] [PubMed]
- Zheng F, Du F, Qian H, et al. Expression and clinical prognostic value of m6A RNA methylation modification in breast cancer. Biomark Res 2021;9:28. [Crossref] [PubMed]
- Galon J, Bruni D. Approaches to treat immune hot, altered and cold tumours with combination immunotherapies. Nat Rev Drug Discov 2019;18:197-218. [Crossref] [PubMed]
- Han D, Liu J, Chen C, et al. Anti-tumour immunity controlled through mRNA m6A methylation and YTHDF1 in dendritic cells. Nature 2019;566:270-4. [Crossref] [PubMed]
- Gershkovitz M, Caspi Y, Fainsod-Levi T, et al. TRPM2 Mediates Neutrophil Killing of Disseminated Tumor Cells. Cancer Res 2018;78:2680-90. [Crossref] [PubMed]
- Zhang W, Shen Y, Huang H, et al. A Rosetta Stone for Breast Cancer: Prognostic Value and Dynamic Regulation of Neutrophil in Tumor Microenvironment. Front Immunol 2020;11:1779. [Crossref] [PubMed]
- Karn T, Denkert C, Weber KE, et al. Tumor mutational burden and immune infiltration as independent predictors of response to neoadjuvant immune checkpoint inhibition in early TNBC in GeparNuevo. Ann Oncol 2020;31:1216-22. [Crossref] [PubMed]
- Keup C, Kimmig R, Kasimir-Bauer S. Liquid Biopsies to Evaluate Immunogenicity of Gynecological/Breast Tumors: On the Way to Blood-Based Biomarkers for Immunotherapies. Breast Care (Basel) 2020;15:470-80. [Crossref] [PubMed]
- Gao C, Li H, Liu C, et al. Tumor Mutation Burden and Immune Invasion Characteristics in Triple Negative Breast Cancer: Genome High-Throughput Data Analysis. Front Immunol 2021;12:650491 [Crossref] [PubMed]
(English Language Editor: J. Jones)