Decoding The Epitranscriptome: In Silico Insights Into m6A Regulatory Network In Breast Cancer

Summary

This protocol presents an approach for conducting in silico genetic, molecular, and prognostic analyses of m6A modification regulators by integrating mutation profiles, copy number alterations, gene expression, and clinical outcomes using publicly available datasets from the Cancer Genome Atlas (TCGA), the Genotype-Tissue Expression (GTEx) project, and microarray platforms.

Abstract

N6-methyladenosine (m6A) is the most abundant internal RNA modification in eukaryotic transcripts and plays a critical role in RNA metabolism, gene expression, and cellular homeostasis. Dysregulation of m6A regulators, including “writers,” “erasers,” and “readers”, has been increasingly implicated in cancer biology; however, their comprehensive roles in breast cancer remain to be understood. The primary objective of this methods article is to provide bioinformatics beginners with a step-by-step framework for utilizing publicly available cancer datasets to perform mutational analyses, assess gene expression alterations, and examine their associations with patient survival. As a case study, m6A regulators in breast cancer were analyzed using datasets from the Cancer Genome Atlas (TCGA), the Genotype-Tissue Expression (GTEx) project, and microarray platforms. Transcriptomic profiles were systematically analyzed to demonstrate workflows for evaluating the prognostic relevance of m6A regulatory components in breast cancer. Using this analytical framework, distinct patterns of genetic alterations and differential expression among key m6A regulators were identified. Several regulators, including METTL14, CBLL1, YTHDC1, HNRNPC, HNRNPA2B1, and RBMX, were associated with better patient survival, while YWHAG was associated with poor overall survival. This study provides a comprehensive systems genomics overview of m6A regulatory genes in breast cancer while demonstrating a practical and reproducible web-based bioinformatics workflow. These findings advance the understanding of epitranscriptomic regulation in breast cancer and offer a foundation for the development of novel m6A-based diagnostic and therapeutic strategies.

Introduction

Epitranscriptomic modifications represent an important layer of post-transcriptional gene regulation and contribute to diverse cellular processes and disease states. Among more than 170 RNA modifications identified to date, N6-methyladenosine (m6A) is the most prevalent and well-characterized in eukaryotic mRNAs¹. Installed by “writer” complexes including METTL3/METTL14, removed by “erasers” including FTO and ALKBH5, and interpreted by “reader” proteins including YTH and IGF2BP family members, m6A orchestrates RNA splicing, stability, transport, and translation, thereby influencing key biological processes including development, differentiation, and stress response^2,3.

Alterations in m6A regulatory components have been reported across a broad spectrum of malignancies⁴. In many cancers, aberrant m6A activity drives malignant phenotypes; e.g., elevated expression of METTL3 promotes prostate cancer initiation and progression by modulating the hedgehog pathway and MYC RNA methylation^5,6. Initially found to be implicated in exerting oncogenic effects in acute myeloid leukemia, FTO was shown to drive tumor progression in liver, lung, and colorectal cancers^7,8,9,10. However, context-dependent roles of FTO and ALKBH5 were identified that illustrate the dual nature of m6A-mediated regulation, which can promote both the oncogenic and tumor suppressive signaling^11,12,13,14. M6A readers, including YTHDF1/2/3,  heterogeneous nuclear ribonucleoproteins (hnRNPs), and insulin-like growth factor-2 mRNA-binding proteins (IGF2BP1-3), have also been found to be associated with carcinogenesis^15,16,17.

In breast cancer, increasing evidence suggests that m6A regulators are frequently dysregulated and may be associated with tumor subtypes, immune-related features, and clinical outcomes^18,19. Multiple mechanistic studies position METTL3 as a frequently upregulated pro-oncogenic factor in breast cancer. METTL3-mediated m6A installation can stabilize or enhance translation of transcripts that promote proliferation, epithelial-mesenchymal transition (EMT), metastasis, and chemoresistance²⁰. METTL3 has also been shown to promote breast cancer progression via targeting Bcl-2²¹. ALKBH5 has been implicated in regulating cancer stemness programs through NANOG and other stemness-related molecules, but its influence may vary by tumor context²².

As the list of m6A regulators continues to expand in recent years, an update on how the newly identified regulators might be dysregulated in breast cancer is needed. Table 1 provides a list of m6A regulators that include writers, readers, and erasers of m6A modification. Additionally, novel m6A regulators, including LRPPRC and YWHAG, have been identified with implications in cancer progression^23,24,25. Therefore, a comprehensive genetic and molecular characterization of all known m6A regulators was conducted in breast cancer using tools that can be employed by researchers with limited bioinformatic background.

The objective of this Methods article is to present a step-by-step platform-based bioinformatics protocol for analyzing m6A regulators in breast cancer using publicly available cancer genomics resources. Using datasets from The Cancer Genome Atlas (TCGA) (www.cancer.gov/tcga), the Genotype Tissue Expression (GTEx) project²⁶, and web-based analytical platforms such as cBioPortal and UCSC Xena, this protocol demonstrates reproducible workflows for assessing mutational profiles, gene expression alterations, and association with patient survival. This visualized and accessible approach is intended to facilitate the adoption of epitranscriptomic data analysis by researchers new to cancer bioinformatics.

Protocol

NOTE: The list of genes encoding m6A methylation regulators, categorized as writers, readers, and erasers, is presented in Table 1. All listed genes were included in the subsequent analyses of mutations, expression patterns, and overall survival. All software and tools used in this study are listed in the Table of Materials.

1. Identification of genetic alterations in m6A regulators

1. Access the cBioPortal for cancer genomics. Navigate to the cBioportal website (www.cbioportal.org)^27,28. From the homepage, select the “Query” tab to begin a new analysis.
2. Select the appropriate cancer study and cohort.
3. In the “Select Studies for Visualization and Analysis” search bar, type “Breast Invasive Carcinoma“ and select “Breast Invasive Carcinoma (TCGA, Pan-Cancer Atlas)”.
4. At the bottom, select “Query by Gene”.
   CRITICAL: Ensure that the selected cohort (996 samples) includes both mutation and copy-number alteration (CNA) data.
5. Define the Genetic Query. In the “Enter Genes”, enter the HUGO gene symbols for the full list of m6A regulators under investigation.
   NOTE: Genes can be entered as a list separated by spaces. Under “Select Genomic Profiles”, ensure the following two data types are checked: Mutations and Copy-number Alterations.
6. Under “Select Patient/Case Set,” choose the default sample set that corresponds to the cohort with all profile cases.
7. Click the blue “Submit Query” button.
8. Retrieve and interpret the genetic alteration data. Upon submission, the result will load on the “summary” tab. The central “OncoPrint” visualization provides an immediate overview of genetic alterations across all queried genes in the cohort, which can be downloaded.
9. Next to the OncoPrint, locate the “Cancer Type Summary” plot. This provides a quantitative breakdown of alterations across breast cancer subtypes.
10. Perform a Pan-Cancer Analysis. Return to the cBioPortal homepage and select “TCGA PanCancer Atlas Studies” under the “Query” tab.
11. In the gene input box, enter the same list of m6A regulator genes.
12. Click “Submit Query” and on the results page, navigate to the “Cancer Types Summary” tab. This provides a pan-cancer view.

2. Comparative transcriptomic analysis of m6A regulators using UCSC Xena.

1. Access the UCSC Xena Platform. Navigate to the UCSC Xena website (https://xena.ucsc.edu)²⁹.
2. From the homepage, click on the “Launch Xena” button to enter the main analysis browser.
3. In the Xena browser, click “DATA SETS”.
4. Among the datasets, select “TCGA TARGET GTEx”. This contains uniformly processed RNA-Seq data from the TCGA and GTEx projects' normal tissues.
5. On the next page, click “VISUALIZE”.
6. Define the Phenotype (sample group) variable. In the “Select Your First Variable”, select “Main Category” in the Phenotypic data type.
7. Click "TO SECOND VARIABLE”. Then, in the Genomic data type, tick “Gene Expression” in the dataset. Add the gene list in the “Add Gene or Position” box. Click “Done”.
8. Visualize expression patterns with a heatmap.
9. To separate the breast (TCGA+GTEx) samples from the TCGA TARGET GTEx, type ”Breast” and use the filter option to keep samples.
10. The heatmap is now visible and can be downloaded as a PDF.
11. Generate comparative box plots for individual genes. To quantify and visualize expression differences for a specific gene, use the “View as chart”. Using this option, data can be viewed as a box plot, a dot plot, and a violin plot, comparing the expression distribution between the two sample groups.
12. Use “Download as PDF” option to download the charts.
13. Statistical significance (p-value) can be obtained by clicking on “STATISTICS”.

3. Assessing prognostic significance of m6A regulators using Kaplan-Meier Plotter.

1. Access the Kaplan-Meier Plotter tool. Navigate to the Kaplan-Meier Plotter website (https://kmplot.com/analysis)³⁰.
2. From the homepage, select the “breast cancer” tab to initiate an analysis specific to breast cancer datasets.
3. Configure the gene query for a single gene.
4. In the primary input section, locate the “Gene symbol” box.
5. Enter the official symbol of the m6A regulator gene to be analyzed (e.g., METTL3).
   CRITICAL: Directly below the gene input box, locate and enable the checkbox for “Only JetSet best probe set”. This ensures the most reliable and specific microarray probe is automatically selected for your gene, optimizing data quality and reproducibility.
6. Define the survival analysis parameters. In the “Survival” section, select “Overall survival (OS)” as the primary endpoint for this analysis. The tool will automatically utilize data from 1880 breast cancer patients when this setting is selected.
7. Ensure the “Split patients by” option is set to “median”. This will stratify patients into two equal groups; high-expression and low-expression, based on the median expression value of the queried gene across all samples.
8. “Follow-up threshold” can be used to select the follow-up period. For this study, 180 months was selected.
9. Generate and interpret the Kaplan-Meier Plot.
10. Click the “Draw Kaplan-Meier Plot” button.
11. A new window will load, displaying the survival curve.
12. Interpret the key plot elements; The X-axis indicates time in months, the Y-axis shows the probability of overall survival, the two colored lines represent the survival curves for the high-expression (red) and low-expression (black) patient groups. The log-rank P-value is displayed, indicating the statistical significance of the difference between the two survival curves.

Results

Mutational landscape of m6a methylation regulators in breast cancer

In an earlier study on the genomic analysis of TCGA datasets, recurrent mutations in several genes encoding regulators of DNA methylation were reported³¹. In the present study, cBioPortal was utilized to analyze the “Breast Invasive Carcinoma (TCGA, PanCancer Atlas)” dataset in order to examine mutational profiles of genes encoding the writers, readers, and erasers of m6A RNA methylation. This analysis revealed diverse genetic alterations among breast cancer patients, with alteration frequencies varying substantially across genes-from 0.4% in CNBP and RBM15B to 12% in VIRMA (Figure 1A). Gene amplification represented the most common alteration, while additional events included deep deletions, base substitutions, and multiple concurrent alterations. Notably, alterations in genes regulating m6A-related functions were detected in 476 patients (48% of the cohort) (Figure 1B), underscoring the significance of m6A modification dynamics in breast cancer. Although the frequency of different alteration types varied, such mutations were observed across all molecular subtypes of breast cancer (Figure 1C). For validation, PIK3CA, TP53, CDH1, and GATA3 were included as reference control genes (Figure 1A). Strikingly, alterations in the m6A regulatory machinery were not restricted to breast cancer. Analysis of 10,967 samples from 10,953 patients from 32 studies in the TCGA Pan-Cancer Atlas revealed conserved mutational patterns in a wide range of cancer types. It has been recently shown that the m6A pathway is frequently altered in prostate cancer (PCa) and overall exerts a pro-oncogenic role³². These findings indicate that mutations affecting genes encoding the writers, readers, and erasers of m6A RNA modification are a common feature across multiple cancers (Figure 2).

Aberrant gene expression profiles in breast cancer

Emerging evidence highlights transcriptomic disruptions as key contributors to tumorigenesis, with aberrant gene expression offering potential as biomarkers in breast cancer. To investigate this, transcript levels of genes regulating m6A modification were analyzed using data from the TCGA and the Genotype Tissue Expression (GTEx) project representing normal breast tissue. As illustrated in Figure 3A, various m6A-associated genes exhibited significant dysregulation in breast cancer samples. Both upregulation and downregulation were observed in tumor tissues compared to normal controls. METTL3 and WTAP, both components of the writer complex, were downregulated among other genes, while several other genes, including VIRMA, YTHDF1, and YTHDF3, were upregulated. Figure 3B further delineates the differential expression profiles of individual genes across the TCGA and GTEx cohorts. Collectively, these findings indicate that genes encoding m6A methylation writers, readers, and erasers undergo extensive transcriptional deregulation in breast cancer, underscoring their potential relevance in disease progression.

m6A machinery genes and their role in patient prognosis

Following the observation that genetic alterations and gene expression changes are highly prevalent among cancer patients, the prognostic relevance of these expression changes in breast cancer was investigated. Utilizing the Kaplan-Meier (KM) Plotter tool³⁰, which integrates microarray datasets, overall survival (OS) was assessed in a cohort of 1880 breast cancer patients according to the expression of m6A regulator genes. This analysis revealed that elevated expression of METTL14, CBLL1, YTHDC1, HNRNPC, HNRNPA2B1, and RBMX was significantly associated with improved overall survival. In contrast, YWHAG overexpression correlated with poor survival outcomes (Figure 4). As controls, CCND2 and TOP2A, known markers of better and poor prognosis, respectively, were included. Other genes encoding m6A regulators did not show statistically significant correlations with patient survival (Supplementary figure). These findings highlight a subset of m6A methylation regulatory genes with potential utility in breast cancer prognostication.

Figure 1: Genetic alterations in m6A writers, readers, and erasers genes in breast cancer. (A) The distribution of alterations across 996 breast cancer patients is shown, with each grey line representing an individual case. The color-coded bars denote different alteration types, including missense mutations, deep deletions, amplifications, in-frame mutations, and truncating mutations. Well-characterized genes, PIK3CA, TP53, CDH1 , and GATA3, are included as positive controls due to their established mutation frequencies. (B) Overall alteration frequency for m6A regulatory genes across the patient cohort. (C) Genetic alteration patterns in m6A regulator genes by breast cancer subtypes. Please click here to view a larger version of this figure.

Figure 2: Frequency of genetic alterations in genes encoding m6A writers, readers, and erasers across diverse cancer types. The analysis is based on data from the TCGA pan-cancer atlas, comprising  10,967 samples from 10,953 patients across 32 cancer studies. Please click here to view a larger version of this figure.

Figure 3: Expression anomalies in genes encoding m6A writers, readers, and erasers. (A) The overexpression (red bars) and underexpression (blue bars) of all the genes are displayed. Data from GTEx and TCGA were used to compare normal vs. breast cancer samples. (B) This figure presents a comparison of individual gene expression in normal vs breast cancer patients. Xena employs Welch's t-test to determine the p-values for each gene. Please click here to view a larger version of this figure.

Figure 4: Expression profiles of m6A writers, readers, and erasers and their association with prognosis in breast cancer. Kaplan- Meier survival curves depict overall patient survival, with the X axis indicating time (months), and the Y-axis showing overall survival probability. Red lines represent the high-expression group, while black lines represent the low-expression group. The patients were stratified based on the median gene expression levels. p-values were determined using the Log-Rank test. Please click here to view a larger version of this figure.

Supplementary Figure: Members of m6A regulators do not show significant correlation with patient overall survival, as shown by Kaplan- Meier survival curves. Red lines represent the high-expression group, while black lines represent the low-expression group.Please click here to download this file.

Type	Gene Symbol
Writers	METTL3
	METTL14
	ZC3H13
	WTAP
	RBM15
	RBM15B
	METTL16
	CBLL1
	KIAA1429/VIRMA
Readers	YTHDF1
	YTHDF2
	YTHDF3
	YTHDC1
	YTHDC2
	HNRNPA2B1
	HNRNPC
	HNRNPG/RBMX
	IGF2BP1
	IGF2BP2
	IGF2BP3
	CNBP
	ELAVL1
	SND1
	PRRC2A
	PRRC2B
	PRRC2C
	EIF3A
	FMR1
	FXR1
	FXR2
	LRPPRC
	MSI2
Erasers	ALKBH5
	FTO

Table 1: Genes encoding writers, readers, and erasers of m6A. Table 1 provides an overview of the major gene families responsible for installing, recognizing, and removing m6A modification in eukaryotic RNA.

Discussion

This Method's article provides a comprehensive, accessible, and integrated workflow for the systematic multi-omics profiling and clinical translation of any gene signature in cancer research, demonstrated here through the analysis of m6A RNA methylation regulators in breast cancer. By combining these major public bioinformatics platforms, this approach enables researchers to efficiently progress from genomic discovery to clinically relevant hypotheses without requiring advanced computational expertise.

The primary strength of this methodology is its modular, hypothesis-generating pipeline. The protocol guides the user through a logical sequence; first identifying which genes are genetically altered (using cBioPortal), then assessing expression dysregulation in a batch-corrected environment (using UCSC Xena), and finally evaluating the clinical impact of that dysregulation on patient survival (using Kaplan-Meier Plotter). This stepwise analysis, from DNA to RNA to clinical outcome, effectively prioritizes candidate genes for further study. For instance, applying this workflow to m6A regulators efficiently pinpointed genes like YWHAG (frequent alteration, prognostic for poor survival) as high-priority targets for functional validation.

The protocol’s design for pan-cancer analysis further enhances its utility, allowing researchers to rapidly determine if a molecular signature is specific to one cancer or is a shared feature of tumorigenesis, as was observed with widespread alterations in m6A machinery in this study. Pan-cancer analysis revealed that mutations affecting m6A writers, readers, and erasers were not confined to breast cancer but are shared across multiple malignancies. This aligns with accumulating evidence that aberrant m6A regulation is a distinctive feature of oncogenesis across diverse tumor types, as it influences multiple hallmarks of cancer and physiological processes, including RNA splicing, stability, translation, and non-coding RNA activity³³.

This methodological approach is highly adaptable. While demonstrated with m6A regulators, the identical workflow can be immediately applied to characterize immune checkpoint genes, metabolic enzymes, or novel gene signatures from RNA-seq experiments in any cancer type available within these databases. This step-by-step format lowers the barrier for wet-lab scientists to perform sophisticated in silico analyses, accelerating the transition from genomic data to biological insight.

In conclusion, this protocol provides a robust framework for the contextualization of cancer-related genes. In addition to delineating the mutational landscape, the findings revealed bidirectional expression changes in the regulators of m6A. This finding underscores the complexity of the epitranscriptome and reinforces the established paradigm of context-dependent function, as exemplified by the dual roles reported for METTL3 and YTHDF family proteins in different cancers^34,35. m6A axis has also been shown to play a role in regulating proliferation, metastasis, and immune evasion in triple-negative breast cancer^36,37. Interestingly, the survival analysis identified a subset of m6A regulators with prognostic significance. Elevated expression of METTL14, CBLL1, YTHDC1, HNRNPC, HNRNPA2B1, and RBMX was associated with favorable outcomes, whereas YWHAG expression correlated with poor overall survival. These findings support the potential clinical utility of m6A regulators as prognostic biomarkers. CBLL1 was also identified as one of the factors with a favorable prognosis in an earlier study³⁸. However, the current analysis incorporating updated members of m6A regulators identified additional members, like RBMX and YWHAG, with better and worse overall survival, respectively. The observation that distinct regulators can predict either adverse or favorable prognosis underscores the dual and context-specific functions of m6A modifications in cancer biology. Although YTHDF1 and YTHDF3 are significantly upregulated in tumors, their lack of correlation with overall survival may reflect functional redundancy among the readers, context-dependent roles across cancer subtypes, or the need to consider their ratio or the net m6A regulatory network rather than individual expression levels. Additionally, although VIRMA exhibited the highest alteration frequency (12%, primarily amplification), its expression did not significantly correlate with overall survival in breast cancer. One possible explanation is that hIgh VIRMA expression indicates only an elevated potential for m6A deposition, not whether the requisite downstream readers or target mRNAs are present to translate this into aggressive tumor behavior. Notably, while YTHDF1 and YTHDF3 were overexpressed in the cohort, YTHDC1 and YTHDC2 were significantly downregulated. This mismatched expression pattern suggests that the functional combination of specific writers and readers required for oncogenic output may not be operative in breast cancer. Thus, despite its high expression, VIRMA may not serve as a dominant driver in this context (Supplementary figure).

The authors acknowledge a primary limitation of this in silico pipeline. The analyses are inherently correlative; they identify strong associations but do not establish mechanistic causality. Such causality can be established through functional genomics or by using small-molecule inhibitors targeting components of the m6A pathway³⁹. Notably, the first METTL3-targeting peptide inhibitor, RSM3, was recently developed and demonstrated anticancer potential in prostate cancer models in vivo⁴⁰. This methodological workflow therefore represents a valuable tool for identifying candidate targets and stratifying patient populations most likely to benefit from these therapeutic interventions.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
cBioPortal	Memorial Sloan-Kettering Cancer Center	https://www.cbioportal.org
Genotype-Tissue Expression (GTEx)	GTEx Consortium	https://gtexportal.org
Kaplan-Meier Plotter	Gyorffy lab/A5 Genetics Ltd	https://kmplot.com
The Cancer Genome Atlas (TCGA)	National Cancer Institute (NCI)	https://www.cancer.gov/tcga
UCSC Xena Browser	University of California Santa Cruz	https://xenabrowser.net