Identification of Endometriosis Pathophysiologic-related Genes Based on Meta-analysis and Bayesian Approach

1. Introduction

Endometriosis is an enigmatic disease characterized by the presence of endometrial glands and stroma outside the uterine endometrium, inhabiting diverse locations, including the ovaries, peritoneum, and extraperitoneal organs [1]. Endometriosis, characterized by a prevalence of as high as 10% among women of reproductive age, gradually progresses and evolves into a chronic inflammatory condition. Medical treatments for endometriosis, such as gonadotropin-releasing hormone agonists, dienogest, and oral contraceptives, along with surgical treatments often cause a reduction in functional ovarian tissue and ovarian reserve, making pregnancy challenging for both natural conception and assisted reproductive techniques [1,2].

Endometriosis manifests as a complex condition with an intricate etiology, weaving together a tapestry of hormonal, immunological, and environmental influences [1]. The preeminent conceptualization of diverse pathophysiological frameworks emerged from the work of John A. Sampson in 1927 [3]. He postulated that the endometrial tissue enters the peritoneal cavity through retrograde menstruation, a process in which expelled endometrial fragments traverse the fallopian tubes, affix, and infiltrate the peritoneal epithelium. Inadequate immune responses may permit the persistence and proliferation of these implants, thereby preventing effective clearance from the implantation sites [1,4]. Another theory of the pathophysiology of endometriosis involves metaplasia of the coelomic epithelium, which is believed to contribute substantially to the development of endometriosis, primarily induced by environmental factors. Although retrograde menstruation in women is a universal physiological phenomenon during the reproductive period, coupled with the observation that only approximately 10% develop endometriosis, this suggests an interplay of diverse factors, including altered cellular immunity, metastasis, genetic underpinnings, environmental influences, and a complex mode of inheritance marked by interactions between specific genes and the environment [1,4,54].

Given the exceedingly intricate pathophysiology of endometriosis, researchers are yet to establish definitive biomarkers for diagnosis [5]. Although laparoscopic surgery has traditionally been considered the gold standard for diagnosing endometriosis [11,12], it is no longer routinely performed for diagnostic purposes. With medical treatment, independent of surgical intervention, constituting a substantial portion of therapeutic approaches, the need non-invasive diagnostic methods distinct from historical approaches have increased; however, they require high specificity [2,11]. Similar to other complex diseases of multifactorial origin, the identification of genes associated with high-penetrance endometriosis remains elusive. Genetic exploration holds promise for the identification of critical pathophysiological pathways [2,11]. In this study, we analyzed gene expression data related to endometriosis from open databases and identified genes associated with the presence and severity of the disease using Bayesian analysis (Figure 1).

2. Results

2.1. Data Exploration and Selection

Obstetrics and gynecology specialists selected data from the Gene Expression Omnibus (GEO) for the meta-analysis. GSE6364 was used to analyze gene expression in endometrial tissue from patients with moderate-to-severe endometriosis and a control group. Similarly, GSE73622 examined gene expression in endometrial tissue samples from a patient-control study in which researchers reprogrammed collected cells into endometrial mesenchymal stem cells (MSC) and subsequently differentiated them into endometrial stromal fibroblasts (SF) using progesterone. GSE141549 explored gene expression in endometrial and peritoneal tissues from patients with endometriosis and healthy women as well as in various lesions within the peritoneal cavities of patients. Uniquely, data from 115 patients were analyzed on a separate platform (GPL10558), while data including from additional 53 controls were analyzed on the GPL13376 platform. Analysts excluded data from the control group in GSE141549-GPL13376 if medication usage or menstrual cycle details were not recorded.

2.2. Differential Expression Analysis

A laboratory medicine specialist conducted gene expression analyses across all five GSE datasets. Given that GSE6364 contains data with gene expression values exceeding 8,000, the analysis utilized log-transformed values. Each GSE sample was trimmed using the NormalizeQuantile. Trimmed data were subjected to principal component analysis (PCA) to identify any batch effects, and any identified batches were designated as confounders during the limma analysis [13].

GSE141549-GPL10558 and GSE141549-GPL13376 include the disease grade results for each patient. These two datasets can be used to analyze disease progression. To analyze the genes associated with endometriosis, we adjusted for GSE141549-GPL13376 and generated a binomial distribution for the patient-control groups instead of using continuous variables for disease grades. Consequently, the GSE6364, GSE73622-MSC, GSE73622-SF, and GSE141549-GPL13376 datasets acquired binomial distributions, allowing us to calculate fold-change (FC) values for the presence of endometriosis. Enrichment analyses were performed using Gene Ontology [14] and the Kyoto Encyclopedia of Genes and Genomes [15]. When common genes between differentially expressed genes (DEGs) and genes within the enrichment pathway were revealed, the false discovery rate (FDR) was calculated. A variety of pathways, including lipid metabolism, protein degradation, cell metabolism, and mRNA signaling (Figure S1–6), were identified in each GSE dataset.

Figure 2 compares the FC distribution across the five GSE datasets and the FC correlation for the same genes identified in the two GSE datasets. A positive FC correlation was observed only between the GSE73622-MSC and GSE73622-SF datasets, which used the same platform as the experimental group. However, GSE141549-GPL13376 and GSE141549-GPL10558, although both were sourced from endometriosis tissues of the same patients and used different platforms, showed a negative correlation in FC for the same observed genes. This difference, attributed to GSE141549-GPL13376 focusing on presence and GSE141549-GPL10558 focusing on severity, led us to separate the meta-analyses for presence and severity.

2.3. Meta-Analysis

The results of the meta-analysis for presence, which are based on binomial groups, cannot be directly used for the severity meta-analysis. However, by transforming the controls and patients into continuous groups with values of 0 and 1, respectively, they were included in the severity analysis. Additionally, the GSEs used for presence analysis were indirectly incorporated into the severity analysis. Each gene had log fold-change (logFC) and standard error (SE) data for each GSE, which were used to perform the meta-analysis using METAL based on the inverse variance-weighted average method (IVW) [16]. Analysis using METAL provided the z-scores and P for each gene. DEGs which p less than 0.05 were related to cellular responses and protein degradation through pathway enrichment analysis (Figure S7–8). Genes that were significant (P < 0.05) had a z-score absolute value greater than 1.96. Approximately, 160 genes had a z-score of less than 1.96 in both meta-analyses (Figure 3).

2.4. Bayesian Analysis of Endometriosis Severity Related Genes

For the 160 genes filtered by the P obtained from the meta-analyses results, Bayesian analysis was applied to identify endometriosis-related upstream genes. Based on a scoring matrix consisting of five types of prior knowledges: Genome wide association study (GWAS) for SNP-associated with the development of endometriosis [6]; human transcription factors (TFs) catalog that affect transcription of pathological genes [7]; SNP-associated with gene expression (eSNP) as known as expression quantitative trait loci (eQTL) that affect the amount of pathological protein expression [8]; DigSee from disease-gene database [9]; and interactome database containing protein-protein interactions (PPI) of proteins made from pathological genes with other proteins [10], the genes were scored based on the number of datasets in which they were included, and the confirmed genes with higher scores were selected as priority genes.

The final list contained 24 genes, with the highest priority PPARA and HLA-DQB1 genes colored in purple. Magenta was selected for three Bayesian datasets, including EP300, MAP2K6, ZC3HAV1, EIF2S1, VIM, ZNF436, and MRRF, and the green color was set for the next ranking including SETBP1, CEP152, TSPAN14, GNG5, BRD4, RPS11, GDI2, MAPK7, TXN, UTP15, ZNF134, ZNF304, ZNF786, ZNF24, and ZNF550. (Figure 4). Higher levels indicated a stronger association with the pathogenesis of endometriosis.

2.5. Correlation and Network Analysis

Pearson correlation analysis of the 24 genes was conducted based on the t-values identified from the six GSE datasets, and 21 gene pairs with P below 0.05 were identified (Figure 5b). This analysis revealed a major cluster consisting of 19 genes and several individual genes surrounding these clusters (Figure 5c). The major cluster was centered on HLA-DQB1 and ZNF24. In particular, HLA-DQB1 scored high in the Bayesian analysis (Figure 4). PPARA, which received the highest score in Bayesian analysis, was depicted as the surrounding gene paired with ZNF134 and UTP15. The ZNF family, which includes highly upstream genes (Figure 5a), lies together, whereas ZNF24 is mainly centered on HLA-DQB1.

3. Discussion

The treatment of endometriosis remains challenging even with the use of specific combinations of medical, surgical, and psychological approaches. Infiltration into various organs, including the uterus, ovaries, and pelvic reproductive organs as well as the rectosigmoid colon, bowel, peritoneum, and diaphragm, complicates management and contributes to high recurrence rates. Meta-analyses and systematic reviews have shown that pain recurs in 20.5% of cases at three years and 43.5% at five years, whereas the recurrence of lesions larger than 10 mm is reported in 9% and 28% of cases, respectively [12]. Severe forms of endometriosis, such as deep infiltrative endometriosis, cause chronic pelvic pain, dysmenorrhea, and dyspareunia, significantly reducing quality of life. It also leads to infertility, which affects women across multiple life stages [5]. As a chronic inflammatory disease with high recurrence rates, severe forms of endometriosis, including deep and extra-abdominal endometriosis, require new therapeutic approaches beyond the traditional medicosurgical treatments. Our analysis provides a novel perspective on the severe forms of endometriosis that continue to recur despite existing treatments, and offers an alternative understanding of this condition.

The comparative analysis of gene expression profiles between pathological and healthy tissues is a robust methodology that provides valuable insights into the underlying cellular processes implicated in the etiology of various diseases. This approach serves as the cornerstone for comprehending the molecular intricacies that drive pathological conditions [17]. Previous studies have analyzed the changes in gene expression associated with endometriosis using various microarray platforms. Some studies have compared gene expression between endometriotic lesions and normal endometrial tissues [18,19]. In certain studies, comparisons were made between gene expression in patients with endometriosis and healthy controls [20,21].

HLA-DQB1, also known as major histocompatibility complex class II DQ beta 1, is located on the short arm of chromosome 6 at position 21.31. Fagerberg et al. demonstrated that HLA-DQB1 was expressed at low levels in all tissues, particularly in the lungs, lymph nodes, and spleen [22]. Typically, specific genotypes of HLA-DQB1 been studied with respect to the development of autoimmune diseases [23,24]. However, a recent study by Xu et al. revealed that increased methylation of HLA-DQB1 is associated with rheumatoid arthritis [25]. In our study, Bayesian analysis was conducted following differential expression analysis, allowing for the identification of significant genes without determining whether they were upstream or downstream. This approach has yielded important insights, suggesting that endometriosis may also be classified as an autoimmune disease, indicating a significant achievement in understanding its pathogenesis.

PPARA, referred to as peroxisome proliferator-activated receptor alpha, is located on the long arm of chromosome 22 at position 13.31. Fagerberg et al. revealed that PPARA is weakly expressed in all organs, but high levels are observed in the kidney, heart, and small intestine [22]. The activation of PPARA is well-known to be associated with lipid metabolism [26], and lipid agonists are being explored for their therapeutic effects in various diseases based on this mechanism [27]. PPARA scored highly in the Bayesian analysis as a gene linked to obesity. In this study, PPARA was identified to be downstream of endometriosis. However, given that lipid metabolism has not been strongly associated with endometriosis, the connection between PPARA and endometriosis is thought to be more related to angiogenesis in endometriotic lesions, as recently reported by Pergialiotis et al. [28], rather than lipid metabolism.

The zinc finger protein family has been linked, either individually or collectively, to a variety of diseases, including adenocarcinoma, squamous cell carcinoma, and Parkinson's disease [29,30,31,32], highlighting the lack of a clearly defined pathogenic mechanism. However, in this study, several ZNF family members were identified upstream and formed part of the main cluster in the network analysis. Well-known genes such as HLA-DQB1 and PPARA, which seem unrelated, were found to be downstream, whereas genes from the ZNF family, for which the specific pathogenic mechanisms remain unknown, were upstream. This intricately woven network, in which seemingly disparate genes come together in a synergistic relationship, was revealed in this study. This suggests that endometriosis should be viewed from a completely new and different perspective than that of previous studies.

4. Materials and Methods

4.1. Pre-Processing

The expression values of the downloaded GSE datasets were adjusted to ensure that the maximum and minimum values were within a consistent range. When the maximum values were excessively large, log transformation was applied. Because a log transformation cannot be performed for negative values, an appropriate constant was added to all values to ensure that the minimum value was greater than zero before the transformation. The data were normalized to ensure that each sample exhibited a consistent distribution of expression values [13]. PCA was conducted to examine the clustering patterns of the samples and assess the relationship between various factors included in the sample information. Factors showing strong correlations were identified as potential confounders.

4.2. Differenctial Expression Analysis

After reviewing the sample information, patient and control groups were designated. For differential expression analysis, factors identified as confounders during preprocessing were incorporated into the model design, and the limma package was used for analysis [13]. The results provided FC between the patient and control groups for each gene, presented as log values. Additionally, P and t-values were analyzed. As the product of the t-value and SE equals logFC, the SE was manually calculated. P and logFC were used to generate volcano plots.

4.3. Meta-Analysis

The GSE6364, GSE73622-MSC, and GSE73622-SF datasets provide gene expression values for both healthy individuals and patients, allowing for binomial differential expression analysis. GSE141549-GPL13376 also included both healthy individuals and patients with the disease stage recorded, enabling both binomial and continuous differential expression analysis. GSE141549-GPL10558 included only patients, allowing for continuous differential expression analysis based on the disease stage. Meta-analysis of the presence of endometriosis was performed using four GSE datasets suitable for binomial analysis. For endometriosis severity, a meta-analysis was conducted by integrating two GSE datasets suitable for continuous differential expression analysis with the results of the endometriosis-presence meta-analysis. Because the present meta-analysis results pertained to the control and patient groups, they were reassigned as 0 and 1, respectively, and then transformed into a continuous scale for meta-analysis. A meta-analysis was conducted using previously obtained logFC and SE values from the METAL program [16]. p and z-scores were calculated for each gene.

4.4. Bayesian Analysis

Bayesian analysis was used to identify the upstream genes associated with endometriosis through evidence collection. The five databases used as evidence included GWAS, TFs, eQTL, DigSee, and PPIs.

The statistical summary related to endometriosis obtained from the GWAS catalog contained information regarding genetic variation by calculating the relationship between SNPs and disease occurrence using logistic regression [6]. This database identifies SNP-trait associations through literature searches, and then extracts reported traits, significant SNP-trait associations, and sample metadata.

The TF catalog was derived from the study by Lambert et al. [7], who identified human transcription factors and their functional characterization based on the important role of TFs as master regulators of chromatin and transcription. This database contains 2,765 genes for which proteins were examined manually by assembling lists of inferred TFs from multiple sources, and 1,639 genes were identified as actual TFs based on validated experiments.

Genomic variants in uterine tissues extracted using the gene expression quantification method at the transcriptomic level were retrieved from the eQTL catalog [8]. This database provides statistical results for cis-eQTLs, which indicate the association between gene expression and SNPs, calculated by linear regression.

DigSee is a search engine that applies text-mining methods to MEDLINE abstracts for cancer-related research on gene-disease interactions [9]. By entering the keyword endometriosis in the search engine, 847 disease-related genes were identified.

The PPI network information was acquired from the STRING database [10]. STRING is a large database and search engine that collects materials on functional associations from publicly accessible sources, such as individual high-throughput studies. Seven lines of evidence support STRING's use of public data as an information base: three prediction streams based on genomic contextual information, and one each for co-expression, text mining methods, experiment-based data, and previously curated pathway and protein complex knowledge.

4.5. Correlation and Network Analysis

Pearson's correlation analysis was conducted on the 24 upstream genes identified in the Bayesian analysis, using t-values from the six differential expression analyses (Figure 2) and z-scores from the two meta-analyses (Figure 3a and 3b). In the correlation analysis, only the top 5% and bottom 5% of the derived degrees were used to form the edges of each node, which were then applied to the network analysis.

5. Conclusions

Selecting specific genes and validating their pathophysiological roles in analyzing gene-disease associations requires extensive time and labor [33]. This study integrated the results from multiple gene expression studies and conducted differential analyses of gene presence in patient-control comparisons and severity associations in stage analysis. This approach systematically narrowed down 160 genes from over 20,000 based on evidence. Furthermore, Bayesian analysis was employed to identify upstream genes. The regulatory genes identified through this process demonstrated substantial statistical significance from the perspective of data analytics. Although experimental validation of the roles of these genes in disease development is still required, this study established a strong foundation by identifying key genes for future research, making a significant contribution to the field.