Identification of Endometriosis Pathophysiologic-Related Genes Based on Meta-Analysis and Bayesian Approach

Abstract

Endometriosis is a complex disease with diverse etiologies, including hormonal, immunological, and environmental factors; however, its exact pathogenesis remains unknown. While surgical approaches are the diagnostic and therapeutic gold standard, identifying endometriosis-associated genes is a crucial first step. Five endometriosis-related gene expression studies were selected from the available datasets. Approximately, 14,167 genes common to these 5 datasets were analyzed for differential expression. Meta-analyses utilized fold-change values and standard errors obtained from each analysis, with the binomial and continuous datasets contributing to endometriosis presence and endometriosis severity meta-analysis, respectively. Approximately 160 genes showed significant results in both meta-analyses. For Bayesian analysis, endometriosis-related single nucleotide polymorphisms (SNPs), the human transcription factor catalog, uterine SNP-related gene expression, disease–gene databases, and interactome databases were utilized. Twenty-four genes, present in at least three or more databases, were identified. Network analysis based on Pearson’s correlation coefficients revealed the HLA-DQB1 gene with both a high score in the Bayesian analysis and a central position in the network. Although ZNF24 had a lower score, it occupied a central position in the network, followed by other ZNF family members. Bayesian analysis identified genes with high confidence that could support discovering key diagnostic biomarkers and therapeutic targets for endometriosis.

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 is recognized as a complex condition influenced by hormonal, immunological, and environmental factors [1]. In 1927, John A. Sampson proposed that retrograde menstruation allowed endometrial tissue to enter the peritoneal cavity, adhere, and invade the epithelium, with inadequate immune responses enabling its persistence and growth [1,3,4]. Another theory suggested that coelomic epithelium metaplasia, driven by environmental factors, played a key role in its development. Although retrograde menstruation occurred universally in reproductive aged women, only about 10% developed endometriosis, indicating the involvement of additional factors such as altered immunity, metastasis, genetic predispositions, environmental triggers, and a complex interplay of gene–environment interactions [1,4,5].

Given the exceedingly intricate pathophysiology of endometriosis, researchers have yet to establish definitive biomarkers for diagnosis [5]. Although laparoscopic surgery was traditionally considered the gold standard for diagnosing endometriosis [6,7], 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 for non-invasive diagnostic methods distinct from historical approaches has increased; however, they require high specificity [2,6]. Similarly 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,6].

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 by combining multiple genome-wide endometriosis-related transcriptomes (Figure 1). However, there are crucial challenges of gene-level studies, such as heterogeneity among datasets, high cost, and relatively small dataset size, which could provide false positive or negative findings [8,9]. We implemented meta-analysis [10], a traditional method to successfully deal with the genetic study’s challenges, to identify reliable endometriosis-related biomarkers (Figure 1). Recently, a Bayesian approach, which integrates computational results with known knowledge (database and published papers), was used to screen key genes involving disease pathophysiology [11,12,13,14]. We utilized the Bayes method, also called convergence functional genomics, to narrow down the hub genes involving the pathophysiology of endometriosis.

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 (MSCs) and subsequently differentiated them into endometrial stromal fibroblasts (SFs) 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 from 53 additional 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. Indentification of Differentially Expressed Genes

A laboratory medicine specialist conducted gene expression analyses across all five GSE datasets. For handling gene expression data, R (R Foundation for Statistical Computing, Vienna, Austria) was used. Given that GSE6364 contains data with gene expression values exceeding 8000, the analysis utilized log-transformed values. Each GSE sample was trimmed using normalizeQuantile function. 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 [20].

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 [21] and the Kyoto Encyclopedia of Genes and Genomes [22]. 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 (Figures S1–S6), 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 on DEGs

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 the 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) [10]. Analysis using METAL provided the z-scores and p for each gene. DEGs with p less than 0.05 were related to cellular responses and protein degradation through pathway enrichment analysis (Figures S7 and S8). 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. The scoring matrix consisted of the following five types of prior knowledges: genome-wide association study (GWAS) for SNPs associated with the development of endometriosis [15]; human transcription factors (TFs) catalog the affect transcription of pathological genes [16]; SNPs associated with gene expression (eSNP), known as expression quantitative trait loci (eQTL), that affect the amount of pathological protein expression [17]; DigSee from disease–gene database [18]; and interactome database containing protein–protein interaction (PPI) of proteins made from pathological genes with other proteins [19]. 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 being the 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. Network Construction by Correlation Analysis

Pearson correlation analysis of the 24 genes was conducted based on the t-values identified from the 6 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 the 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 Bayesian method includes two main parts: the optimization of likelihood and verification by prior knowledge. We arranged the results of the meta-analysis as the likelihood part and multiple lines of evidence, including GWAS [15], TF [16], eQTL [17], and the disease–gene database [18], as the prior part. Different datasets derived DEG results using samples collected from entirely distinct individuals, leading to varying coefficient values for the same gene. However, this heterogeneity was utilized as a hyperparameter to structure each dataset into models of identical format. By integrating these models, the meta-analysis addressed uncertainties observed across raw datasets by optimizing the likelihood [23]. Furthermore, leveraging datasets obtained from actual experiments, rather than entirely virtual foundations, enabled the rapid construction of models based on prior knowledge [24]. This approach provided robust evidence for building models aimed at validating research hypotheses. Bayesian’s flexible nature allows it to be incorporated into meta-analysis-based studies and helps reduce the time and cost of biomarker discovery.

The Bayesian approach could reduce false positive results by using prior pieces of knowledge that are sufficiently reliable without special statistical calculations [13,14]. Additionally, integrating DEG analysis results from various datasets through meta-analysis involved iterative modeling, which, as a combination of prior and likelihood, delivered more reliable outcomes compared to simple statistical analyses [24]. Furthermore, validating the selected genes using credible genetic network databases enhanced the quality of pathophysiological identification. By adopting sufficiently well-founded signals, the Bayesian methods can show another way of looking at results beyond “p-value hysteria and fantasy” [25].

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 [7]. 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 [26]. 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 [27,28]. In certain studies, comparisons were made between gene expression in patients with endometriosis and healthy controls [29,30].

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 [31]. Typically, specific genotypes of HLA-DQB1 were studied with respect to the development of autoimmune diseases [32,33]. However, a recent study by Xu et al. revealed that increased methylation of HLA-DQB1 is associated with rheumatoid arthritis [34]. In our study, the 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 [31]. The activation of PPARA is well known to be associated with lipid metabolism [35], and lipid agonists are being explored for their therapeutic effects in various diseases based on this mechanism [36]. 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 was not 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. [37], rather than lipid metabolism.

PPARα activates genes required for fatty acid β-oxidation and represses IL-6 gene expression by regulating the AP-1 and NF-kB pathway [38]. HLA-DQB1 is expressed on antigen-presenting cells, which present exogenously derived peptides to CD4+ T cells to elicit an immune response. The etiological link between HLA-DQB1 and endometriosis was found at the DNA level [39]. From the results of this study, it remained unclear exactly how and to what extent PPARA and HLA-DQB1 contributed to the presence and severity of endometriosis. However, based on their association with pathways involved in immune response regulation, it was suggested that endometriosis was linked to an abnormal immune state. Given the intricately intertwined signaling pathways of the immune response, it was likely that the precise pathophysiology of endometriosis had not yet been fully elucidated. The findings of this study may have served as a guiding light in navigating this uncertain landscape.

The zinc finger protein family was linked, either individually or collectively, to a variety of diseases, including adenocarcinoma, squamous cell carcinoma, and Parkinson’s disease [40,41,42,43], 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. ZNF436 was selected in meta-analysis and ranked high in the Bayesian approach (Figure 4). ZNF436 was suggested to function as a negative regulator of mitogen-activated protein kinases (MAPKs), in which the signaling pathway regulated cell proliferation [44,45]. The ectopic expression of ZNF304 in primary CD4+ cells silences human immunodeficiency viral transcription [46]. The upstream of the ZNF family was hypothesized to contribute to immune dysregulation and an inadequate immune response, potentially leading to the onset and progression of endometriosis.

4. Materials and Methods

4.1. Dataset Selection Criteria

The selection of appropriate gene expression datasets from the GEO involved narrowing down results through several steps. The initial search term was “endometriosis,” and studies focusing on Homo sapiens were retained. To analyze the entire human genome, the study type was specified as “Expression profiling by array”. Datasets containing sufficient sample sizes with both patient and normal control groups were reviewed. Comparative studies evaluating responses to drugs were excluded to focus on genes associated with the pathophysiological characteristics of the disease. GSE6364 was identified as a dataset containing 21 patient samples and 16 control samples, providing gene profiling data of the endometrium. GSE73622 included 22 patient samples and 28 controls, analyzing the gene expression of dedifferentiated cells from the endometrium in endometriosis. GSE141549, with 335 patient samples and 73 controls, featured the largest sample size and provided transcriptome analysis across various stages of endometriosis and multiple tissue types.

4.2. 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 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 [20]. The log transformation used for normalization employed natural logarithms. Since the distribution of gene expression values varied across samples, errors could arise during analysis. To address this, normalizeQuantiles function from the limma package was applied to complete the normalization process [20]. 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.3. Differential Expression Analysis

After reviewing the sample information, patient and control groups were designated. For differential expression analysis, factors identified as confounders during pre-processing were incorporated into the model design, and the limma package was used for analysis [20]. Variables identified as confounders through PCA were handled as factors using the as.factor function within R’s model.matrix, effectively mitigating batch effects. This design enabled the use of the lmFit function to construct linear models for continuous variables in the samples. The resulting linear models were analyzed statistically using the eBayes function. The topTable function organized the statistical outcomes into a data frame, facilitating the selection of genes. The results provided FC values 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.4. 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 [10]. p and z-scores were calculated for each gene.

4.5. 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 [15]. 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. [16], 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 2765 genes for which proteins were examined manually by assembling lists of inferred TFs from multiple sources, and 1639 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 [17]. 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 [18]. 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 [19]. 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.6. 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 6 differential expression analyses (Figure 2) and z-scores from the 2 meta-analyses (Figure 3a,b). 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 [47]. 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.

Author Contributions

Conceptualization, J.K. and T.L.; methodology, T.L., S.H. and K.A.; validation, K.A. and J.O.; formal analysis, T.L.; investigation, K.A.; resources, J.K. and S.J.C.; data curation, J.K. and T.L.; writing—original draft preparation, J.K. and T.L.; writing—review and editing, K.A., J.K. and J.O.; visualization, K.A.; supervision, S.J.C.; project administration, Y.U. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available in Gene Expression Omnibus (GEO) at https://www.ncbi.nlm.nih.gov/geo/ (accessed on 15 March 2024), reference numbers GSE6364, GSE73622, GSE141549. These data were derived from the following resources available in the public domain: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse6364 (accessed on 15 March 2024); https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse73622 (accessed on 15 March 2024); https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse141549 (accessed on 15 March 2024).

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. A scheme for the identification of genes associated with endometriosis. (STEP 1) A multi-disciplinary team including obstetrics and gynecology doctors, bioinformatics experts, and laboratory medicine doctors conducted literature and database reviews to select candidate transcriptome datasets related to endometriosis. The GSE6364 dataset analyzed gene expression in endometrial tissues from control groups and patients with moderate-to-severe endometriosis (indicated by red arrows). Similarly, GSE73622 differentiated endometrial mesenchymal stem cells and endometrial stromal fibroblasts from the endometrial tissue using progesterone were utilized to analyze gene expression (indicated by green arrows). GSE141549 examined endometrial and peritoneal tissues as well as various lesions in the peritoneal cavity of patients with endometriosis and healthy women (indicated by yellow arrows). (STEP 2) The algorithm conducts a differential expression analysis on datasets from the three GSE studies. A meta-analysis using the inverse variance-weighted method identified genes that were ultimately associated with endometriosis. (STEP 3) The Bayes method is implemented to identify the endometriosis upstream genes using five lines of external evidences: genome-wide association studies (GWAS) on single nucleotide polymorphism (SNP) associated with disease occurrence [15]; transcription factors (TF) affecting the transcription (T) of pathological genes [16], which identifies genes related to pathophysiology; expression quantitative trait loci (eQTL) reflecting the expression of SNP related to pathological proteins [17]; disease–gene database, known as DigSee [18]; and protein–protein interaction (PPI) data showing interactions between pathological proteins (P) and other proteins (Pr) [19].

Figure 2. A comparison of global endometriosis-related signatures among five transcriptomic datasets. The volcano plots placed in the diagonal panels indicate cohort-specific endometriosis-related signatures (FC between two conditions or among three statuses). The specified order of GSE datasets, GSE6364, GSE73622-MSC, GSE73622-SF, GSE141549-GPL13376, and GSE141549-GPL10558 was followed, with volcano plots from differential analysis for each dataset arranged diagonally in the center. In the volcano plots, dots represent genes, with those having a p of 0.05 or higher shown in pale gray. Among the significant genes, red indicates upregulated genes, while blue means downregulated genes. Values placed in the upper triangle panels indicate the degree of correlation matched with the lower triangle plots evaluated by PCC and SCC. The p for each correlation coefficient appears in parentheses. A gray dot of correlational plots located in the lower triangle panels denotes a single identical gene contained in both gene expression datasets, and the matched red dashed line describes the correlation between two dataset obtained from linear regression. Abbreviations: FC, fold-change; MSC, mesenchymal stem cell; SF, stromal fibroblast; PCC, Pearson’s correlation coefficient; SCC, Spearman’s correlation coefficient; Bi, binomial; Or, ordinal.

Figure 3. The identification of the endometriosis-related genes based on meta-analysis. (a) The IVW-based meta-analysis combined the four lists of global FC values between the presence and absence of endometriosis of four gene expression datasets (GSE6364, GSE73622-MSC, GSE73622-SF, and GSE141549-GPL13376). (b) The IVW combined the two lists of endometriosis severity-related signatures (GSE141549-GPL13376 and GSE141549-GPL10558) with the meta-analyzed FC values related to endometriosis presence (a). (c) The scatter plot indicates the comparison of meta-analyzed z-scores between presence and severity. A point represents an individual gene, and red points (about 160 genes) denote cases having 1.96 or more z-scores in both binomial (x-axis) and step-wise (y-axis) statuses. Abbreviations: IVW, inverse variance-weighted average method; FC, fold-change.

Figure 4. The identification of regulators related to endometriosis based on the Bayesian approach. Cells placed on the upper panel indicate candidate regulators. The colors of each cell are determined based on the strength of external evidence. Cells on the lower panel indicate the prior knowledge section of the Bayesian approach. The prior knowledge includes multiple lines of external evidence, including endometriosis-related SNPs (GWAS) [15], the human TF catalog [16], uterus eSNP (eQTL) [17], the disease–gene database (DigSee) [18], and the interactome database [19]. Abbreviations: SNP, single nucleotide polymorphism; GWAS, genome-wide association study; TF, transcription factor; eSNP, SNP-associated gene expression; eQTL, expression quantitative trait loci.

Figure 5. The identification of network patterns among the 24 upstream endometriosis meta genes. (a) The first and second rows represent the z-score calculated by IVW in METAL [10]. The z-scores in the first and second rows are the meta-analyzed combined endometriosis presence and severity signatures, respectively. The order of the genes was determined by calculating the average of two z-scores. The third to sixth rows reveal notable binomial distribution (Bi) of endometriosis presence. The seventh to eighth rows show significant ordinal distribution (Or) of endometriosis severity. (b) Colors, direction, density, and shape were colored differently according to Pearson’s correlation coefficient (PCC) values between the upstream meta genes. *, **, and *** indicate p less than 0.05, 0.01, 0.001, respectively. (c) The gene network is constructed based on the PCC matrix (b) and cut-off (p [less than] 0.05). Red, orange, yellow, and khaki points indicate genes having 4 or more, 3, 2, and 1 or less of the number of interactions, respectively. The size of gene labels varies according to their significance. The color of the edges is determined by PCC, with positive values shown in pink and negative values in sky blue, while smaller p results in thicker edges.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ijms26010424/s1.

References

1. Giudice, L.C.; Kao, L.C. Endometriosis. Lancet; 2004; 364, pp. 1789-1799. [DOI: https://dx.doi.org/10.1016/S0140-6736(04)17403-5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15541453]

2. Raja, M.H.R.; Farooqui, N.; Zuberi, N.; Ashraf, M.; Azhar, A.; Baig, R.; Badar, B.; Rehman, R. Endometriosis, infertility and MicroRNA’s: A review. J. Gynecol. Obstet. Hum. Reprod.; 2021; 50, 102157. [DOI: https://dx.doi.org/10.1016/j.jogoh.2021.102157]

3. Sampson, J.A. Peritoneal endometriosis due to the menstrual dissemination of endometrial tissue into the peritoneal cavity. Am. J. Obstet. Gynecol.; 1927; 14, pp. 422-469. [DOI: https://dx.doi.org/10.1016/S0002-9378(15)30003-X]

4. Burney, R.O.; Giudice, L.C. Pathogenesis and pathophysiology of endometriosis. Fertil. Steril.; 2012; 98, pp. 511-519. [DOI: https://dx.doi.org/10.1016/j.fertnstert.2012.06.029] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22819144]

5. Horne, A.W.; Missmer, S.A. Pathophysiology, diagnosis, and management of endometriosis. BMJ; 2022; 379, e070750. [DOI: https://dx.doi.org/10.1136/bmj-2022-070750]

6. Szaflik, T.; Romanowicz, H.; Szyłło, K.; Kołaciński, R.; Michalska, M.M.; Samulak, D.; Smolarz, B. Analysis of long non-coding RNA (lncRNA) UCA1, MALAT1, TC0101441, and H19 expression in endometriosis. Int. J. Mol. Sci.; 2022; 23, 11583. [DOI: https://dx.doi.org/10.3390/ijms231911583] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36232884]

7. Bozdag, G. Recurrence of endometriosis: Risk factors, mechanisms and biomarkers. Women’s Health; 2015; 11, pp. 693-699. [DOI: https://dx.doi.org/10.2217/whe.15.56]

8. Neupane, B.; Loeb, M.; Anand, S.S.; Beyene, J. Meta-analysis of genetic association studies under heterogeneity. Eur. J. Hum. Genet.; 2012; 20, pp. 1174-1181. [DOI: https://dx.doi.org/10.1038/ejhg.2012.75]

9. Wookward, A.A.; Urbanowicz, R.J.; Naj, A.C.; Moore, J.H. Genetic heterogeneity: Challenges, impacts, and methods through an associative lens. Genet. Epidemiol.; 2022; 46, pp. 555-571. [DOI: https://dx.doi.org/10.1002/gepi.22497]

10. Willer, C.J.; Li, Y.; Abecasis, G.R. METAL: Fast and efficient meta-analysis of genomewide association scans. Bioinformatics; 2010; 26, pp. 2190-2191. [DOI: https://dx.doi.org/10.1093/bioinformatics/btq340]

11. Hägg, S.; Skogsberg, J.; Lundström, J.; Noori, P.; Nilsson, R.; Zhong, H.; Maleki, S.; Shang, M.M.; Brinne, B.; Bradshaw, M. et al. Multi-organ expression profiling uncovers a gene module in coronary artery disease involving transendothelial migration of leukocytes and LIM domain binding 2: The Stockholm Atherosclerosis Gene Expression (STAGE) study. PLoS Genet.; 2009; 5, e1000754. [DOI: https://dx.doi.org/10.1371/journal.pgen.1000754] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19997623]

12. Koplev, S.; Seldin, M.; Sukhavasi, K.; Ermel, R.; Pang, S.; Zeng, L.; Bankier, S.; Di Narzo, A.; Cheng, H.; Meda, V. et al. A mechanistic framework for cardiometabolic and coronary artery diseases. Nat. Cariovasc. Res.; 2022; 1, pp. 85-100. [DOI: https://dx.doi.org/10.1038/s44161-021-00009-1] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36276926]

13. Xu, M.; Zhang, D.F.; Luo, R.; Wu, Y.; Zhou, H.; Kong, L.L.; Bi, R.; Yao, Y.G. A systematic integrated analysis of brain expression profiles reveals YAP1 and other prioritized hub genes as important upstream regulators in Alzheimer’s disease. Alzheimers Dement.; 2018; 14, pp. 215-229. [DOI: https://dx.doi.org/10.1016/j.jalz.2017.08.012]

14. Bertsch, B.; Ogden, C.A.; Le-Niculescu, H.; Kuczenski, R.; Niculescu, A.B. Convergent functional genomics: A Bayesian candidate gene identification approach for complex disorders. Methods; 2005; 37, pp. 274-279. [DOI: https://dx.doi.org/10.1016/j.ymeth.2005.03.012] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16308156]

15. Backman, J.D.; Li, A.H.; Marcketta, A.; Sun, D.; Mbatchou, J.; Kessler, M.D.; Benner, C.; Liu, D.; Locke, A.E.; Balasubramanian, S. et al. Exome sequencing and analysis of 454,787 UK Biobank participants. Nature; 2021; 599, pp. 628-634. [DOI: https://dx.doi.org/10.1038/s41586-021-04103-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34662886]

16. Lambert, S.A.; Jolma, A.; Campitelli, L.F.; Das, P.K.; Yin, Y.; Albu, M.; Chen, X.; Taipale, J.; Hughes, T.R.; Weirauch, M.T. The human transcription factors. Cell; 2018; 172, pp. 650-665. [DOI: https://dx.doi.org/10.1016/j.cell.2018.01.029]

17. Lonsdale, J.; Thomas, J.; Salvatore, M.; Phillips, R.; Lo, E.; Shad, S.; Hasz, R.; Walters, G.; Garcia, F.; Young, N. et al. The Genotype-Tissue Expression (GTEx) project. Nat. Genet.; 2013; 45, pp. 580-585. [DOI: https://dx.doi.org/10.1038/ng.2653] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23715323]

18. Kim, J.; Kim, J.-J.; Lee, H. An analysis of diease-gene relationship from Medline abstracts by DigSee. Sci. Rep.; 2017; 7, 40154. [DOI: https://dx.doi.org/10.1038/srep40154]

19. Szklarczyk, D.; Kirsch, R.; Koutrouli, M.; Nastou, K.; Mehryary, F.; Hachilif, R.; Gable, A.L.; Fang, T.; Doncheva, N.T.; Pyysalo, S. et al. The STRING database in 2023: Protein-protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res.; 2023; 51, pp. D638-D646. [DOI: https://dx.doi.org/10.1093/nar/gkac1000] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36370105]

20. Lee, T.; Hwang, S.; Seo, D.M.; Shin, H.C.; Kim, H.S.; Kim, J.-Y.; Uh, Y. Identification of cardiovascular disease-related genes based on the co-expression network analysis of genome-wide blood transcriptome. Cells; 2022; 11, 2867. [DOI: https://dx.doi.org/10.3390/cells11182867]

21. The Gene Ontology Consortium. The gene ontology resource: 20 years and still GOing strong. Nucleic Acids Res.; 2019; 47, pp. D330-D338. [DOI: https://dx.doi.org/10.1093/nar/gky1055]

22. Kanehisa, M.; Goto, S. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res.; 2000; 28, pp. 27-30. [DOI: https://dx.doi.org/10.1093/nar/28.1.27] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10592173]

23. Lotfi, S.; Izmailov, P.; Benton, G.; Goldblum, M.; Wilson, A.G. Bayesian model selection, the marginal likelihood, and generalization. Proceedings of the 39th International Conference on Machine Learning; Baltimore, MD, USA, 17–23 July 2022; Volume 162, pp. 14223-14247. Available online: https://proceedings.mlr.press/v162/lotfi22a.html (accessed on 31 December 2024).

24. Mikkola, P.; Martin, O.A.; Chandramouli, S.; Hartmann, M.; Pla, O.A.; Thomas, O.; Pesonen, H.; Corander, J.; Vehtari, A.; Kaski, S. et al. Prior knowledge elicitation: The past, present, and future. Bayesian Anal.; 2024; 19, pp. 1129-1161. [DOI: https://dx.doi.org/10.1214/23-BA1381]

25. Spanagel, R. Convergent functional genomics in addiction research—A translational approach to study candidate genes and gene networks. Silico Pharmacol.; 2013; 1, 18. [DOI: https://dx.doi.org/10.1186/2193-9616-1-18] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25505662]

26. Trevino, V.; Falciani, F.; Barrera-Saldaña, H.A. DNA microarrays: A powerful genomic tool for biomedical and clinical research. Mol. Med.; 2007; 13, pp. 527-541. [DOI: https://dx.doi.org/10.2119/2006-00107.Trevino] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17660860]

27. Borghese, B.; Mondon, F.; Noël, J.-C.; Fayt, I.; Mignot, T.-M.; Vaiman, D.; Chapron, C. Gene expression profile for ectopic versus eutopic endometrium provides new insights into endometriosis oncogenic potential. Mol. Endocrinol.; 2008; 22, pp. 2557-2562. [DOI: https://dx.doi.org/10.1210/me.2008-0322]

28. Hever, A.; Roth, R.B.; Hevezi, P.; Marin, M.E.; Acosta, J.A.; Acosta, H.; Rojas, J.; Herrera, R.; Grigoriadis, D.; White, E. et al. Human endometriosis is associated with plasma cells and overexpression of B lymphocyte stimulator. Proc. Natl. Acad. Sci. USA; 2007; 104, pp. 12451-12456. [DOI: https://dx.doi.org/10.1073/pnas.0703451104]

29. Tamaresis, J.S.; Irwin, J.C.; Goldfien, G.A.; Rabban, J.T.; Burney, R.O.; Nezhat, C.; DePaolo, L.V.; Giudice, L.C. Molecular classification of endometriosis and disease stage using high-dimensional genomic data. Endocrinology; 2014; 155, pp. 4986-4999. [DOI: https://dx.doi.org/10.1210/en.2014-1490]

30. Zhao, L.; Gu, C.; Ye, M.; Zhang, Z.; Han, W.; Fan, W.; Meng, Y. Identification of global transcriptome abnormalities and potential biomarkers in eutopic endometria of women with endometriosis: A preliminary study. Biomed. Rep.; 2017; 6, pp. 654-662. [DOI: https://dx.doi.org/10.3892/br.2017.902]

31. Fagerberg, L.; Hallström, B.M.; Oksvold, P.; Kampf, C.; Djureinovic, D.; Odeberg, J.; Habuka, M.; Tahmasebpoor, S.; Danielsson, A.; Edlund, K. et al. Analysis of the human tissue-specific expression by genome-wide integration of transcriptomics and antibody-based proteomics. Mol. Cell Proteomics; 2014; 13, pp. 397-406. [DOI: https://dx.doi.org/10.1074/mcp.M113.035600]

32. Miyadera, H.; Tokunaga, K. Associations of human leukocyte antigens with autoimmune disease: Challenges in identifying the mechanism. J. Hum. Genet.; 2015; 60, pp. 697-702. [DOI: https://dx.doi.org/10.1038/jhg.2015.100]

33. Ramgopal, S.; Rathika, C.; Malini, R.P.; Murali, V.; Arun, K.; Balakrishnan, K. Critical amino acid variations in HLA-DQB1* molecules confers susceptibility to autoimmune thyroid disease in south India. Genes. Immun.; 2019; 20, pp. 32-38. [DOI: https://dx.doi.org/10.1038/s41435-017-0008-6]

34. Xu, J.; Chen, H.; Sun, C.; Wei, S.; Tao, J.; Jia, Z.; Chen, X.; Lv, W.; Lv, H.; Tang, G. et al. Epigenome-wide methylation haplotype association analysis identified HLA-DRB1, HLA-DRB5 and HLA-DRB1 as risk factors for rheumatoid arthritis. Int. J. Immunogenet.; 2023; 50, pp. 291-298. [DOI: https://dx.doi.org/10.1111/iji.12637] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37688529]

35. Seitz, H.K.; Moreira, B.; Neuman, M.G. Pathogenesis of alcoholic fatty liver a narrative review. Life; 2023; 13, 1662. [DOI: https://dx.doi.org/10.3390/life13081662] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37629519]

36. Wang, N.; Zhao, Y.; Wu, M.; Li, N.; Yan, C.; Guo, H.; Li, Q.; Li, Q.; Wang, Q. Gemfibrozil alleviates cognitive impairment by inhibiting ferroptosis of astrocytes via restoring the iron metabolism and promoting antioxidant capacity in type 2 diabetes. Mol. Neurobiol.; 2024; 61, pp. 1187-1201. [DOI: https://dx.doi.org/10.1007/s12035-023-03589-0]

37. Pergialiotis, V.; Frountzas, M.; Fasoulakis, Z.; Daskalakis, G.; Chrisochoidi, M.; Kontzoglou, K.; Perrea, D.N. Peroxisome proliferator-activated receptor alpha (PPAR-α) as a regulator of the angiogenic profile of endometriotic lesions. Cureus; 2022; 14, e22616. [DOI: https://dx.doi.org/10.7759/cureus.22616] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35371629]

38. Lin, Y.; Wang, Y.; Li, P. PPARα: An emerging target of metabolic syndrome, neurodegenerative and cardiovascular disease. Front. Endocrinol.; 2022; 13, 1074911. [DOI: https://dx.doi.org/10.3389/fendo.2022.1074911] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36589809]

39. Ishii, K.; Takakuwa, K.; Kashima, K.; Tamura, M.; Tanaka, K. Associations between patients with ednometriosis and HLA class II; the analysis of HLA-DQB1 and HLA-DPB1 genotypes. Hum. Reprod.; 2003; 18, pp. 985-989. [DOI: https://dx.doi.org/10.1093/humrep/deg192]

40. Zhu, L.; Tu, D.; Li, R.; Li, L.; Zhang, W.; Jin, W.; Li, T.; Zhu, H. The diagnostic significance of the ZNF gene family in pancreatic cancer: A bioinformatics and experimental study. Front. Genet.; 2023; 14, 1089023. [DOI: https://dx.doi.org/10.3389/fgene.2023.1089023] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37396042]

41. Jia, D.; Li, L.; Wang, P.; Feng, Q.; Pan, X.; Lin, P.; Song, S.; Yang, L.; Yang, J. ZNF24 regulates the progression of KRAS mutant lung adenocarcinoma by promoting SLC7A5 translation. Front. Oncol.; 2022; 12, 1043177. [DOI: https://dx.doi.org/10.3389/fonc.2022.1043177]

42. Wang, Z.; Sun, A.; Yan, A.; Yao, J.; Huang, H.; Gao, Z.; Han, T.; Gu, J.; Li, N.; Wu, H. et al. Circular RNA MTCL1 promotes advanced laryngeal squamous cell carcinoma progression by inhibiting C1QBP ubiquitin degradation and mediating beta-catenin activation. Mol. Cancer; 2022; 21, 92. [DOI: https://dx.doi.org/10.1186/s12943-022-01570-4]

43. Santiago, J.A.; Potashkin, J.A. Evaluation of RNA blood biomarkers in individuals at risk of Parkinson’s disease. J. Parkinson’s Dis.; 2017; 7, pp. 655-660. [DOI: https://dx.doi.org/10.3233/JPD-171155] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28922168]

44. Li, Y.; Du, X.; Li, F.; Deng, Y.; Yang, Z.; Wang, Y.; Pen, Z.; Wang, Z.; Yuan, W.; Zhu, C. et al. A novel zinc-finger protein ZNF436 suppresses transcriptional activities of AP-1 and SRE. Mol. Biol. Rep.; 2006; 33, pp. 287-294. [DOI: https://dx.doi.org/10.1007/s11033-006-9019-5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17089209]

45. Krasnopolsky, S.; Kuzmina, A.; Taube, R. Genome-wide CRISPR knockout screen identifies ZNF304 as a silencer of HIV transcription that promotes viral latency. PLoS Pathog.; 2020; 16, e1008834. [DOI: https://dx.doi.org/10.1371/journal.ppat.1008834] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32956422]

46. Zhang, W.; Liu, H.T. MAPK signal pathways in the regulation of cell proliferation in mammalian cells. Cell Res.; 2002; 12, pp. 9-18. [DOI: https://dx.doi.org/10.1038/sj.cr.7290105]

47. Strande, N.T.; Riggs, E.R.; Buchanan, A.H.; Ceyhan-Birsoy, O.; DiStefano, M.; Dwight, S.S.; Goldstein, J.; Ghosh, R.; Seifert, B.A.; Sneddon, T.P. et al. Evaluating the clinical validity of gene-disease associations: An evidence-based framework developed by the clinical genome resource. Am. J. Hum. Genet.; 2017; 100, pp. 895-906. [DOI: https://dx.doi.org/10.1016/j.ajhg.2017.04.015]