Exploring diagnostic m6A regulators in endometriosis

Data preprocessing

PCA analysis showed that EU samples could be easily distinguished from EC samples in the GSE7305, GSE7307, and E-MTAB-694 datasets. In the GSE51981 dataset, most NM samples could be distinguished from the EU samples (Figure 1A–1D). However, in the GSE6364 dataset, it was very difficult to differentiate EU and NM samples given that the central points of the PCA results almost overlapped (Figure 1E). To improve the quality of this analysis, we excluded the GSE6364 dataset from the subsequent combined analysis; this was due to non-satisfactory PCA performance.

Thus, we processed and merged raw CEL files for GSE7305, GSE7307, GSE51981, and E-MTAB-694. Then, we obtained a normalized matrix, that featured 34 normal (NM), 127 eutopic (EU), and 46 ectopic (EC) samples of the endometrium (Figure 1F). PCA analysis, carried out after adjustment for batch effects, indicated that the EC samples were obviously distinct from the NM samples while some EU samples were mixed with NM or EC samples. These findings suggested the innate heterogeneity of the EU samples (Figure 1G).

The mRNA expression landscape of m6A regulators in EMs

The expression patterns of 20 m6A regulators in 34 NM, 127 EU, and 46 EC samples are shown in Figure 2A. Exception for WATP, most m6A regulators were dysregulated among three groups. When compared with NM samples, most m6A regulators were significantly downregulated in EU samples, including KIAA1429, METTL14, METTL3, RBM15, ZC3H13, YTHDC1, YTHDC2, YTHDF2, YTHDF3, HNRNPC, HNRNPA2B1, RBMX, and FTO. In contrast, RBM15B, IGF2BP1, IGF2BP3, and ALKBH5 were significantly upregulated in EU samples. However, when compared with EU samples, only a small number of m6A regulators showed significant changes in EC samples: METTL3, YTHDF1, YTHDF2, YTHDF3, HNRNPC, and HNRNPA2B1 were downregulated, while IGF2BP2 and FTO were upregulated (Supplementary Table 2). Moreover, METTL3, YTHDF2, YTHDF3, HNRNPC, HNRNPA2B1, and FTO, were identified as differentially expressed m6A regulators that were shared between the EU vs. NM matrix and the EC vs. EU matrix (Figure 2B).

The generation of a heatmap for the expression of the 20 m6A regulators revealed significant heterogeneity in expression in both the EU and EC samples (Figure 2C). Spearman correlation analysis of the 20 m6A regulators indicated that the highest positive correlation coefficient was detected between HNRNPC and HNRNPA2B1 (r = 0.92, p < 0.05) (Figure 2D).

The relationship between selected m6A regulators and EMs

Differentially expressed m6A regulators that were shared between the EU vs. NM matrix and the EC vs. EU matrix (METTL3, YTHDF2, YTHDF3, HNRNPA2B1, HNRNPC, and FTO) were chosen for clinical correlation analysis. This analysis indicated that HNRNPC and HNRNPA2B1 showed differential expression in different r-AFS stages in EU samples (Figure 3A). However, no m6A regulators showed any significant changes across different r-AFS stages in EC samples (Figure 3B), in different phases of the menstrual cycle, in different age groups for both EU and EC samples (Figure 3C–3F), and in different races for the EU samples (Figure 3G). Moreover, we found that METTL3 and YTHDF2 were differentially expressed in different subtypes of EMs in EC samples: METTL3 was downregulated in ovarian EC samples when compared to peritoneal EC samples, while YTHDF2 was upregulated in the ovarian EC samples (Figure 3H).

The diagnostic value of HNRNPC and HNRNPA2B1 in EMS

Given that HNRNPA2B1 and HNRNPC were associated with the severity of EMs in EU samples, we selected these two molecules to evaluate their diagnostic value when compared to MKI67, CDH1 (E-cadherin), ACTA2 (α-SMA), PGR, ESR1, and ESR2, in EMs tissue samples. As shown in Figure 4A, HNRNPA2B1, HNRNPC, MKI67, CDH1 (E-cadherin), PGR, and ESR1, were all downregulated in the EU vs. NM matrix and the EC vs. EU matrix, while ACTA2 (α-SMA) and ESR2 were upregulated in the EU vs. NM and EC vs. EU matrices. Hence, 1/ACTA2 and 1/ESR2 were used to perform ROC analyses.

ROC analysis identified that HNRNPA2B1 (AUC = 0.892) and HNRNPC (AUC = 0.890) demonstrated higher diagnostic potential over MKI67 (AUC = 0.692), CDH1 (AUC = 0.667), 1/ACTA2 (AUC = 0.651), PGR (AUC = 0.801), ESR1 (AUC = 0.744) and 1/ESR2 (AUC = 0.744) in the EU vs. NM matrix (Figure 4C). HNRNPA2B1 (AUC = 0.716) and HNRNPC (AUC = 0.710) also exhibited relatively modest diagnostic value compared to MKI67 (AUC = 0.918), CDH1 (AUC = 0.906), 1/ACTA2 (AUC = 0.889), PGR (AUC = 0.908), ESR1 (AUC = 0.964) and 1/ESR2 (AUC = 0.693) in the EC vs. EU matrix (Figure 4B). Moreover, the scatter points showed a more concentrated distribution for HNRNPA2B1 and HNRNPC than for MKI67, CDH1, ACTA2, PGR and ESR1 (Figure 4A).

In addition, analysis of the validation dataset (GSE105764) also showed that HNRNPA2B1 and HNRNPC possessed diagnostic potential for discriminating between EC and EU samples, with an AUC of 0.969 and an AUC of 0.938 compared to MKI67 (AUC = 1.000), CDH1 (AUC = 1.000), 1/ACTA2 (AUC = 1.000), PGR (AUC = 0.969), ESR1 (AUC = 1.000) and 1/ESR2 (AUC = 0.828) (Figure 4D, 4E).

Functional annotation of HNRNPA2B1 and HNRNPC in EMs

To explore the putative function of HNRNPA2B1 and HNRNPC in EMs, we divided the EU and EC samples into high and low expression groups based on their median expression values. We then identified differentially expressed genes (DEGs) in 63 low-expression vs. 64 high-expression HNRNPA2B1 EU samples, 63 low-expression vs. 64 high-expression HNRNPC EU samples, 23 low-expression vs. 23 high-expression HNRNPA2B1 EC samples, and 23 low-expression vs. 23 high-expression HNRNPC EC samples. Volcano plots showed that the DEGs in EU samples were much more abundant than in the EC samples (Figure 5A–5D). GSEA analysis further indicated that HNRNPA2B1 and HNRNPC were involved in similar Biological Processes (BPs) and Reactome pathways, most of which were associated with immune and inflammatory responses (Figure 5E–5H). For example, in the EU samples, the enriched BP terms and Reactome pathways with the highest normalized enrichment score (NES) in the low-expression vs. high-expression HNRNPA2B1 samples, and the low-expression vs. high-expression HNRNPC samples, were both ‘antimicrobial humoral response’ and ‘initial triggering of complement’ (Figure 5E, 5F). The top 3 BP terms in each matrix were also presented as classical GSEA plots (Figure 5I–5L).

We also used the online Enrichr tool to carry out BP and KEGG analysis. This analysis also indicated that HNRNPA2B1 and HNRNPC were associated with immune and inflammatory pathways in both EU and EC samples. Moreover, these two molecules were related to RNA metabolic process, RNA transportation, RNA splicing, and cell cycle pathways (Supplementary Figure 1). Furthermore, 48 shared DEGs were obtained from screened DEGs (|log2FC| > 0.5 and p < 0.05) between low-expression vs. high-expression HNRNPA2B1 and low-expression vs. high-expression HNRNPC in respective EU and EC samples (Figure 5M). Most of these were enriched in Reactome pathways related to tissue remolding, such as ‘molecules associated with elastic fibers’, ‘elastic fiber formation’, ‘extracellular matrix organization’, and immune response such as ‘initial triggering of complement’ (Figure 5N).

Notably, PGR was one of the 48 shared DEGs and was downregulated in the low-expression HNRNPA2B1 and low-expression HNRNPC samples compared to high-expression HNRNPA2B1 and high-expression HNRNPC samples (Figure 5N). PGR was also positively correlated with the expression of HNRNPA2B1 (r = 0.77, p < 0.05) and HNRNPC (r = 0.72, p < 0.05) in both the training and validation datasets (Figure 5O, 5P).

The association of HNRNPA2B1 and HNRNPC with infiltrating immune cells in EMs

ssGSEA analysis indicated that compared with the high-expression HNRNPA2B1 or HNRNPC EU samples, multiple types of immune cells showed significant elevation in the low-expression HNRNPA2B1 or HNRNPC EU samples (Figure 6A, 6B). Similarly, the infiltration scores for immune cells were generally higher in the low-expression HNRNPA2B1 or HNRNPC EC samples compared to the high-expression HNRNPA2B1 or HNRNPC EC samples, although only B cells showed a significant elevation in the low-expression vs. high-expression HNRNPA2B1 EC samples (Figure 6C, 6D). Moreover, the enrichment of immune cells was negatively correlated with the expression of HNRNPA2B1 and HNRNPC in both EU and EC samples (Figure 6E–6H, Supplementary Figure 2).

Furthermore, the MCP-counter method also showed that immune cells were more abundant in the low-expression HNRNPA2B1 or HNRNPC EU and EC samples compared to the high-expression HNRNPA2B1 or HNRNPC EU and EC samples (Supplementary Figure 3). The negative correlation was also observed between MCP-counter scores and the expression of HNRNPA2B1 and HNRNPC in both EU and EC samples (Supplementary Figure 4).

The potential regulatory mechanisms underlying the role of HNRNPA2B1 and HNRNPC in EMs

NetworkAnalyst version 3.0 predicted 9 transcription factors (TFs) for HNRNPA2B1 and 16 for HNRNPC. These predicted TFs were then used to construct a gene-TFs network (Figure 7A). Several of the predicted TFs were differentially expressed in the training dataset, among which SRF, ELK1, USF2, FOXC1, HNF4A were upregulated, and BRCA1, ESR1, YY1, NFYA were downregulated, both in the EU vs. NM matrix and the EC vs. EU matrix. The analysis also showed that TP53, E2F1, GATA2, NRF1, and MEF2A, presented with contrasting trends for these two matrices; for instance, TP53 was upregulated in the EU vs. NM matrix but was downregulated in the EC vs. EU matrix (Figure 7A, 7B and Supplementary Table 3).

Both HNRNPA2B1 and HNRNPC were significantly and positively correlated with the expression levels of NFYA (one of the common predicted TFs for HNRNPA2B1 and HNRNPC), and ESR1 (one of the putative TFs for HNRNPC), both of which were downregulated in the EU vs. NM and EC vs. EU matrices (Figure 7B, 7C). Similar results were observed in the GSE105764 validation dataset (Figure 7D, 7E and Supplementary Table 4).

Data resources

The analysis flowchart is shown in Figure 8. The GEO (https://www.ncbi.nlm.nih.gov/geo/) and Array-Express (https://www.ebi.ac.uk/arrayexpress/) databases were used to mine datasets associated with EMs. The search strategy was previously described by Poli et al. [43] and used the following keywords for the initial screen: ‘endometriosis’ or ‘endometrium’ and ‘tissue’ with ‘homo sapiens’ or ‘human’. Subsequently, we used the ‘platform GPL570’ (Affymetrix Human Genome U133 Plus 2.0 Array) as an additional filter condition to reduce the ‘platform effect’.

Our searches identified five candidate microarray datasets: GSE7305 [44], GSE7307, GSE51981 [45], GSE6364 [46], and E-MTAB-694 [47]. Further screening led to the exclusion of 37 endometrium samples in GSE51981 and 20 endometrium samples in GSE6364 due to the fact that these samples were associated with other forms of uterine or ovarian pathology (e.g., leiomyomata, adenomyosis, and ovarian cyst). Thus, 10 EU and 10 EC samples from GSE7305; 23 EU and 18 EC samples from GSE7307; 77 EU and 34 NM samples from GSE51981, 21 EU and 7 NM samples from GSE6364, and 17 EU and 18 EC samples from E-MTAB-694 were reserved. Most of the RNA-seq datasets associated with EMs usually involved small sample sizes and did not usually include NM, EU, and EC samples simultaneously. Therefore, we selected an independent RNA-seq dataset, GSE105764 [48], from the GEO database as a validation cohort; this dataset included 8 paired EU and EC tissue samples and was based on the GPL20301 platform (Illumina HiSeq 4000) (Table 1).

Clinical information associated with these datasets was retrieved from published articles (Supplementary Table 1). Clinical information was not available for GSE7307 as this data analysis had not been published. Ethical approval was not necessary as our analysis was based on public data resources, and did not involve direct experiments on humans or animals.

Data preprocessing

First, raw CEL files were downloaded from the GEO and Array-Express databases. Second, we adjusted data for background correction and normalization using the robust multi-array average (RMA) method in the R package ‘Oligo’ [49]. Finally, probes were annotated according to the GPL570 annotation file; for this, we used the official Affymetrix website (http://www.affymetrix.com). Before combining data, we then performed principal component analysis (PCA) on each microarray dataset to examine dataset quality. Datasets that showed good performance in the PCA were selected for combing as the training dataset; datasets that showed poor performance in the PCA were excluded. The raw CEL files of the qualifying datasets were then combined in the R package to create a merged and normalized matrix containing NM, EU, and EC samples. We also used the ‘removeBatchEffect’ function in the R package ‘LIMMA’ [50] to adjust batch effects, and then carried out PCA analysis to evaluate the performance of batch effect adjustment.

The mRNA expression landscape of m6A regulators in EMs

Following pretreatment, the expression levels of 20 m6A regulators were evaluated in samples of NM, EU, and EC tissues by the ‘Kruskal.test’ in the R Package. We also used the R package ‘LIMMA’ to identify differentially expressed m6A regulators in two contrasting matrices: EU vs. NM and EC vs. EU [50] with a threshold of p < 0.05. The differentially expressed m6A regulators that were shared between the EU vs. NM matrix and the EC vs. EU matrix were then screened out by intersection analysis. The expression of these 20 m6A regulators was then demonstrated in a heatmap using ‘pheatmap’ R package; this was based on Euclidean distance and hierarchical clustering [51]. We also evaluated the correlations between these m6A regulators in EMs using Spearman’s correlation analysis with a statistical threshold of p < 0.05.

The relationships between the selected m6A regulators and EMs

The differentially expressed m6A regulators that were shared between the EU vs. NM matrix and the EC vs. EU matrix were then used for clinical correlation analysis. Clinical information is provided in Supplementary Table 1. Using this information, we investigated the associations between the selected m6A regulators and the key clinical features of patients with EMs, including the revised American Fertility Society (r-AFS) stages, menstrual cycle phases, age, races, and disease subtypes (ovarian EMs or peritoneal EMs). This analysis was carried out for both EU and EC samples and involved the ‘Wilcox.test’ in the R package with a statistical threshold of p < 0.05.

We did not analyze samples if clinical information was missing. We found that information relating to the race was missing from the EC samples; consequently, we did not perform race analysis on the EC samples. Furthermore, subtype correlation analysis was not performed for the EU samples as many of these EU samples were associated with multiple diagnoses simultaneously; for instance, some patients with ovarian EMs also had peritoneal EMs.

Diagnostic value of selected m6A regulators in EMs

Next, we selected the m6A regulators that were significantly associated with the severity of EMs to perform receiver operating characteristic (ROC) analysis. This allowed us to compare the diagnostic potential of these m6A regulators with certain histological mRNA markers of EMs, including MKI67, CDH1 (E-cadherin), ACTA2 (α-SMA), PGR, ESR1, and ESR2. MKI67 is a classic marker of cell proliferation; proliferative activity is known to be limited in ectopic lesions due to the reduced ratio of epithelial/stromal cells [1, 52]. The process of epithelial-mesenchymal transition (EMT) is often marked by the loss of certain epithelial markers, including CDH1 (E-cadherin), and the acquisition of certain mesenchymal markers, such as ACTA2 (α-SMA). These processes, along with fibroblast-to-myofibroblast trans-differentiation (FMT), are thought to contribute to the increased production of collagen and smooth muscle metaplasia (SMM), thus resulting in fibrosis in EMs [53, 54]. Moreover, previous research has shown that the downregulation of ESR1 and PGR, and the upregulation of ESR2, play a central role in the estrogen-driven inflammation and progesterone resistance in EMs [55]. The results derived from our analyses were validated in the GSE105764 dataset. A p-value < 0.05 was regarded as statistically significant.

Functional annotation of diagnostic m6A regulators in EMs

To explore the putative function of diagnostic m6A regulators in EMs, we divided the EU and EC samples into two groups according to their median expression values: a low expression group and a high expression group. We then used the ‘LIMMA’ R package to identify differentially expressed genes (DEGs) between these two groups. Next, we used the ‘Clusterprofiler’ R package to perform gene set enrichment analysis (GSEA) analysis [56]. And we used two gene sets to act as reference gene sets: ‘c5.bp.v7.0.symbols.gmt’ for Biological Process (BP) terms for Gene Ontology (GO) analysis, and ‘c2.cp.reactome.v7.0.symbols.gmt’ for Reactome Pathway analysis. These gene sets were downloaded from the Molecular Signature Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). P-values < 0.01 and adjusted p-values < 0.05 were considered to be statistically significant.

Then, the DEGs that were identified between the low expression and high expression groups, which showed a |log2FC| > 1 and p < 0.05 in the EU samples, and DEGs with a |log2FC| > 0.5 and a p < 0.05 in the EC samples, were used for BP and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. For this analysis, we used Enrichr (https://amp.pharm.mssm.edu/Enrichr/), a useful online tool for gene annotation. Terms retrieved with a p-value < 0.05 were considered to be statistically significant.

Next, to investigate the common functions of candidate m6A regulators in EU and EC samples, the DEGs identified between the low and high expression groups (with a |log2FC| > 0.5 and a p < 0.05) in the EU samples were intersected with those from the EC samples. These intersecting DEGs were further evaluated by Reactome pathway analysis; this analysis was carried out on the Reactome website (https://reactome.org).

The association between diagnostic m6A regulators and infiltrating immune cells in EMs

To investigate the correlation between diagnostic m6A regulators and infiltrating immune cells in EMs, we performed single-sample gene set enrichment analysis (ssGSEA); this is a method used to evaluate gene set enrichment in a single sample and was carried out using the ‘GSVA’ R package [57]. We used a variety of markers for each type of immune cell in EMs, as described by Vallvé et al. [38]: macrophages type 1 (Mφ1), macrophages type 2 (Mφ2), immature dendritic cells (iDCs), mature dendritic cells (mDCs), Natural killer cells (NK), mast cells (MC), Eosinophils (EN), Neutrophils (NT), B cells, CD8+T cells, CD4+Th17 cells, CD4+Th1 cells, CD4+Th2 cells, Treg cells, NKT cells, and γδ T cells. We made these markers into a ‘.gmt’ file as the reference gene set. The ssGSEA enrichment scores obtained for each type of immune cell were then scaled and compared between the low expression and high expression groups of m6A regulators. Furthermore, we also used the Microenvironment Cell Populations-counter (MCP-counter) method [58] to evaluate the infiltration of immune cells in the low expression and high expression groups of m6A regulators for both EU and EC samples. We also analyzed the relationship between m6A regulators and immune infiltration cells by Spearman correlation; p < 0.05 was considered to be statistically significant.

Potential regulatory mechanisms of diagnostic m6A regulators in EMs

Transcription factors (TFs) play a critical role in the regulation of gene expression. The dysregulation of TFs was frequently described in EMs [1]. To explore the potential regulatory mechanisms of diagnostic m6A regulators in EMs, we predicted TFs and constructed a gene-TF network using Network Analyst version 3.0 (https://www.networkanalyst.ca/) using the JASPAR database (http://jaspar.genereg.net/) as a source. The Gene-TF network was then visualized by software Cytoscape [59]. The expression levels of TFs and their correlation with diagnostic m6A regulators were then examined in both training (‘Kruskal.test’) and validation datasets (‘DEseq2’ package); adjusted p-value < 0.05 was considered to be statistically significant.