Integration of Single Cell and Bulk RNA-Sequencing Reveals Key Genes and Immune Cell Infiltration to Construct a Predictive Model and Identify Drug Targets in Endometriosis

Introduction

Endometriosis is a chronic neuroinflammatory disease strongly associated with chronic pelvic pain in women. It is estimated to affect 6–10% of women of reproductive age,¹ with its incidence rising to 50% among infertile women.² Additionally, epidemiological research indicates that women with endometriosis have a higher likelihood of developing ovarian cancer,³ breast cancer, melanoma, asthma, rheumatoid arthritis, and cardiovascular diseases.⁴

The pathogenesis of endometriosis remains incompletely understood. The primary current treatments include surgical excision of lesions and pharmacotherapy to suppress ovarian hormone production.⁵ However, more than half of the patients who undergo surgery require additional procedures within 5.5 years.⁶ Moreover, endometrioma cystectomy may have long-term and detrimental effects on ovarian reserve.⁷ Regarding pharmacological treatment, current guidelines recommend hormonal therapies, such as combined oral contraceptives, progestins, and GnRH agonists and antagonists.⁸ However, it is important to note that these treatments may cause potential side effects, including bone loss and low estrogen symptoms induced by GnRH analogs.^9,10 Therefore, further elucidation of the pathogenesis of endometriosis and the development of more effective therapeutic strategies that alleviate symptoms without compromising fertility are critical research priorities.^5,11

Various hypotheses, including retrograde menstruation, intracavitary dissemination, and vascular or lymphatic metastasis, have been proposed to explain the mechanisms of endometriosis. However, these theories do not fully capture the diverse pathological subtypes of the disease.¹² Emerging evidence suggests that the cellular and molecular characteristics of the eutopic endometrium in endometriosis patients are pivotal in the disease’s pathogenesis. Studies have documented significant alterations in gene and protein expression within this tissue,^13,14 yet analyses focusing on single markers often fail to capture the full spectrum of regulatory changes. To address this limitation, high-throughput techniques such as microarray and RNA sequencing (RNA-seq) have been employed to provide a comprehensive analysis of transcriptional changes and to identify critical genes and pathways involved in endometriosis.^15,16 Nonetheless, inconsistencies in identifying differentially expressed genes across RNA-seq studies may arise from variations in sample collection, including differences in menstrual cycle phase, endometriosis subtype, and hormonal treatment backgrounds.^17,18 Therefore, it is crucial to standardize these variables in RNA-seq analyses to ensure consistent and reliable results.

Additionally, the endometrium is a dynamic, multicellular tissue that responds to steroid hormones. Its epithelial cells are supported by a stromal layer, which is rich in vasculature and populated by a diverse and fluctuating array of immune cells. During menstruation, the luminal layer undergoes breakdown and is subsequently shed.⁵ The shed endometrial cells include epithelial cells, stromal fibroblasts, vascular cells, and various immune cells, particularly neutrophils, monocytes/macrophages, and uterine natural killer (uNK) cells.¹⁹ While all these cell types may potentially contribute to peritoneal lesions, contingent on their ability to survive in the peritoneal cavity and evade immune surveillance, there remains significant debate regarding which specific cell types play a pivotal role in the pathogenesis of endometriosis. Changes in immune cells within the eutopic endometrium (eg, macrophages, NK cells, T cells, B cells, etc.) have garnered significant attention. Although findings regarding immune cell alterations in the eutopic endometrium in endometriosis are not yet fully consistent across studies,^5,20 existing evidence supports the involvement of immune cell dysfunction in the pathophysiological mechanisms of endometriosis. This dysfunction contributes to the establishment of an unfavorable environment for embryo implantation, the formation of ectopic lesions, fibrosis, and pain.²⁰

Since bulk RNA-seq is typically performed at the tissue or multi-cellular level, the genetic data obtained represents an average from multiple cells or cell types, which overlooks intercellular heterogeneity.²¹ In contrast, single-cell RNA sequencing (scRNA-seq), which has emerged over the past decade, enables the acquisition of genetic information at the single-cell level, revealing the gene structure and expression state of each individual cell. This approach captures the heterogeneity and genetic characteristics of rare cells that traditional methods fail to detect. Additionally, scRNA-seq facilitates the more precise discovery and validation of disease-related biomarkers, particularly those that may be diluted or obscured in conventional techniques.²¹ The advent of scRNA-seq technologies has also provided unprecedented insights into the unique roles of different cell populations in the development of endometriosis.²² By employing scRNA-seq, we can conduct an in-depth analysis of the transcriptomic profiles of individual cell types, thereby enhancing our understanding of their interactions and their collective contribution to the progression of endometriosis.^23,24

In this study, we systematically downloaded and analyzed Bulk RNA-seq and scRNA-seq data from the Gene Expression Omnibus (GEO) database, focusing on proliferative phase endometrial samples from endometriosis patients and healthy controls. We developed a predictive model to assess the risk of endometriosis onset, which was validated using external cohorts to confirm its predictive accuracy. The model integrates a combination of genetic biomarkers capable of predicting the development of endometriosis, with the goal of supporting early diagnosis and guiding personalized treatment and care. Additionally, we characterized the immune cell infiltration landscape within the eutopic endometrium of endometriosis, identified novel regulatory genes strongly associated with endometriosis pathogenesis, and validated these findings through in vitro experiments. Furthermore, we conducted drug target prediction analyses for potential therapeutic interventions. Collectively, our findings offer novel insights that may enhance the clinical management of endometriosis.

In this study, we systematically downloaded and analyzed Bulk RNA-seq and single-cell RNA-seq data from the Gene Expression Omnibus (GEO) database, focusing on proliferative phase endometrial samples from endometriosis patients and healthy controls. First, through in-depth analysis of scRNA-seq data and the calculation of contributions from various cell subtypes, we identified the key role of mesenchymal cells in the pathogenesis of endometriosis. Next, we developed a predictive model that, for the first time, integrated a combination of eight genes (SYNE2, TXN, NUPR1, CTSK, GSN, MGP, IER2, and CXCL12) to predict the onset of endometriosis, and validated its predictive accuracy using external cohorts. The model integrates key genes capable of predicting the development of endometriosis, aiming to support early diagnosis and guide personalized treatment and care. In addition, we characterized the immune cell infiltration landscape in the eutopic endometrium of endometriosis patients, identified novel regulatory genes closely associated with the pathogenesis of endometriosis, and validated these findings through in vitro experiments. Furthermore, we conducted drug target prediction analyses for potential therapeutic interventions. Overall, this study provides novel insights that could enhance the clinical management of endometriosis and holds significant potential for practical applications.

Materials and Methods

Dataset Source

Single-cell RNA sequencing (scRNA-seq) datasets and bulk RNA-seq datasets were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/info/datasets.html). The selection criteria for the datasets were as follows: (1) completeness and quality of sequencing data; (2) comprehensive patient information, with clear documentation of the menstrual cycle phase at the time of sample collection, and exclusion of patients who had used oral contraceptives or GnRH analogs prior to sampling; (3) availability of corresponding in vitro validation for the sequencing data; (4) inclusion of a healthy control group; (5) sufficiently large sample sizes. For scRNA-seq, we selected GSE179640 and GSE213216, which are currently the largest available scRNA datasets for endometriosis, with total sample sizes of 59 and 51, respectively. For bulk RNA-seq, we selected GSE25628 and GSE153739 due to the completeness of their data, comprehensive sample collection information, extensive in vitro validation, and frequent citation in other studies.^25,26 As previously discussed, dynamic changes in gene expression and immune cell infiltration occur in the eutopic endometrium across different phases of the menstrual cycle.^5,20 However, many sequencing studies have not accounted for the menstrual cycle phase in their analyses,^22,27 and there has been limited focus on immune infiltration and gene expression alterations in the proliferative phase of the eutopic endometrium. In this study, we specifically selected sequencing data from the proliferative phase of the eutopic endometrium in the aforementioned datasets to better explore the gene expression profiles and immune cell infiltration during this phase.

The GSE179640 dataset includes scRNA-seq data from the endometrial tissue of 3 healthy controls and the eutopic endometrial tissue of 3 endometriosis patients, while the GSE213216 dataset includes scRNA-seq data from the endometrial tissue of 3 healthy controls and the eutopic endometrial tissue of 7 endometriosis patients, all of which were incorporated into the subsequent analysis. Additionally, the bulk RNA-seq data from the endometrial tissue of 6 healthy controls and the eutopic endometrial tissue of 8 endometriosis patients in the GSE25628 dataset (based on the GPL571 platform), as well as from the endometrial tissue of 3 healthy controls and the eutopic endometrial tissue of 4 endometriosis patients in the GSE153739 dataset (based on the GPL18573), were included in this study for further analysis. A review of the relevant literature confirmed that all eutopic endometrial samples from the above-mentioned patients were in the proliferative phase. Detailed information for the aforementioned datasets was provided in Table S1.

ScRNA-Seq Data Processing

To handle ScRNA-seq data, samples were integrated using the anchors technique in the R package “Seurat”.²⁸ Core cells were obtained by filtering the scRNA-seq data. Cells with significant mitochondrial content (>20%), as well as those with fewer than 200 or more than 6000 genes (filtering out ruptured cells and potential non-singlet cells, respectively), were removed from each sample. The gene expression of core cells was normalized using a linear regression model, and the top 2000 genes with highly variable properties were identified using analysis of ANOVA. Principal component analysis (PCA) was performed on single-cell samples, with the top 20 principal components (PC) chosen for further analysis, and the Harmony algorithm²⁹ was used for batch effect correction. The top 20 principal components were subjected to overall dimensionality reduction analysis using the uniform manifold approximation and projection (UMAP) algorithm,³⁰ which determined the positional relationships between clusters. Subsequently, The R package “singleR”³¹ was initiated to automatically annotate the clusters using reference data from comprehensive databases such as the Human Primary Cell Atlas, PanglaoDB, and ImmuneCellExpressionData. These databases offer detailed information on human cell types, including those relevant to the endometrium. Following this, manual annotation was performed by cross-referencing the marker genes with the CellMarker database³² and recent literature. Any inconsistencies between the automated and manual annotations were resolved by reviewing the marker gene profiles and adjusting the annotations to ensure biological accuracy. This combined approach enhanced the precision and reliability of the cell type identification in our study.

Contribution of Different Cell Subpopulations to Endometriosis

To assess the contribution of each cell subpopulation to endometriosis, this study evaluated changes in both cell abundance and gene expression. Differential gene expression analysis was first conducted to identify the characteristic genes, focusing on the top 100 genes most highly expressed in the control versus endometriosis groups. The differential expression levels (Fold Change, FC) and percentage proportion (PctProp) of these genes in each cell subpopulation were then quantified. The product of FC and PctProp for each gene within each subpopulation was calculated. To normalize the contribution, the square root of FC * PctProp was computed, considering both gene expression intensity and cell representation. This methodology, commonly used in single-cell transcriptomics studies, allows for quantification of the gene contributions from different cell subpopulations.^33,34

Identification and Functional Enrichment Analysis of Differentially Expressed Genes in Eutopic Endometrium of Endometriosis

RNA-seq data from the eutopic endometrium of 6 healthy controls and 8 endometriosis patients in the GSE25628 dataset were processed using the limma tool to compute differentially expressed genes (DEGs) across the groups. The DEGs selection criteria were P.Value < 0.05 and |log2FC| > 1. The DEGs volcano plots and heatmaps were generated using the ggplot2³⁵ and pheatmap³⁶ programs, respectively.

Construction and Validation of Disease Risk Model

To create a novel disease risk prediction model, the machine learning algorithm Least Absolute Shrinkage and Selection Operator (Lasso)³⁷ was used to identify genes indicative of disease risk. The analysis steps are summarized as follows. The intersection of DEGs from the GSE25628 dataset and mesenchymal cell cluster genes obtained from the above scRNA-seq data analysis was defined as candidate genes. The GSE25628 dataset served as the training set (13 samples), while the GSE153739 dataset was used as the validation set (7 samples). Univariate Cox proportional hazards regression analysis was performed on candidate genes in the training set to identify characteristic genes associated with the risk of developing endometriosis. Variables with p-values < 0.05 were included in the least absolute shrinkage and selection operator (LASSO) regression analysis, performed using the R package “glmnet”³⁸ to reduce the number of genes in the final risk model. The disease risk model was constructed using the formula: risk score = gene expression value1 × β1 + gene expression value2 × β2 + ... + gene expression valuen × βn (where gene expression value represents the gene expression level, and β represents the corresponding LASSO regression coefficient). The “survROC” package³⁹ was used to plot receiver operating characteristic (ROC) curves to evaluate the performance of the risk score in predicting the occurrence of endometriosis. Furthermore, the validity of the disease risk model was verified using the internal validation set and external datasets GSE25628 and GSE153739, respectively.

Immuno-Infiltration Analysis

The CIBERSORT algorithm⁴⁰ utilizes deconvolution analysis on immune cell subtype expression matrices through support vector regression to differentiate 22 immune cell types based on 547 biomarkers, including T cells, B cells, plasma cells, and myeloid cell subtypes. To evaluate the immune cell infiltration profiles in the endometrium of endometriosis patients, we uploaded the standardized GSE153739 gene expression matrix to the CIBERSORT website (https://cibersort.stanford.edu/) and applied the LM22 signature matrix to estimate the relative proportions of 22 immune cell types, retaining only data with a p-value < 0.05. Additionally, we applied immune cell signature gene sets derived from previous studies^41,42 and performed single sample gene set enrichment analysis (ssGSEA) to calculate the absolute immune cell infiltration scores for each sample in the GSE153739 dataset. The results were presented as bar plots, and immune cell correlations were calculated using the Psych package.

Functional Pathway Analysis in Endometriosis Using Gene Set Enrichment Analysis and Gene Set Variation Analysis

Using the clusterProfiler package (Version 4.12.0), we conducted GSEA enrichment analysis on target gene sets in endometrium samples from both the endometriosis and healthy control groups to investigate the functional pathway differences associated with endometriosis development. Target marker pathway genes related to immune response, cell signaling, and other relevant processes were obtained from the Molecular Signature Database (MSigDB) (http://www.gsea-msigdb.org/gsea/index.jsp). Subsequently, Gene set variation analysis (GSVA) enrichment analysis, a nonparametric and unsupervised algorithm, was performed on these target gene sets in the endometrial samples from both groups, and the differences in GSVA enrichment scores between the two groups were analyzed using the R package limma.

MicroRNAs Network Construction

MicroRNAs (miRNA) are small non-coding RNAs that have been shown to regulate gene expression by promoting the degradation of mRNAs or inhibiting the translation of mRNAs.⁴³ We further analyzed the key genes identified in the above analysis to determine if there are significant miRNAs involved in their transcription or degradation. Using the miRcode database (http://www.mircode.org), we identified miRNAs related to these key genes (Table S2) and visualized the gene-miRNA network using Cytoscape software (version 3.10.2).

Regulatory Network Analysis of Key Genes

This study used the R package “RcisTarget (version 1.23.1)”⁴⁴ to predict transcription factors for the selected DEGs set identified in the above analysis. All computations were based on motifs, and the normalized enrichment score (NES) for each motif depends on the total number of motifs in the database. In addition to motifs annotated by the source data, further annotations were inferred based on motif similarity and gene sequences. The process of estimating the overexpression of each motif in the gene set consists of two steps. First, the area under the curve (AUC) for each motif-gene set pair is calculated based on the recovery curve generated by ranking the gene set according to the motif. Next, the NES for each motif is determined using the AUC distribution of all motifs within the gene set. We used the RcisTarget.hg19.motifDBs.cisbpOnly.500bp as the motif-gne rankings database.

Connectivity Map Drug Prediction

The Connectivity Map (CMap) is a tool used to analyze the relationship between gene expression profiles and drug responses. It aims to explore the connection between gene expression and drug actions, helping to identify potential therapeutic candidates or elucidate drug mechanisms.^45,46 In this study, we uploaded 150 upregulated and 150 downregulated differentially expressed genes (DEGs) from ectopic endometrial tissues of endometriosis patients and healthy controls to CMap (version 1.1.1.43, https://clue.io).^45,46 CMap compares these gene expression profiles with those of over 1.5 million human cell line samples treated with more than 5,000 bioactive compounds, generating connectivity scores. These scores, ranging from −1 to 1, assess the degree of similarity between the query signature and CMap compounds: positive scores indicate stimulation of the query signature by the compound, while negative scores suggest inhibitory effects. As a result, small molecules with potential therapeutic effects were identified.

Sample Collection and Quantitative Polymerase Chain Reaction (qPCR)

A total of four patients with endometrioma treated by laparoscopic surgery were recruited from the Union Hospital, Tongji Medical College, Huazhong University of Science and Technology, between October 2023 and March 2024. All patients had a pathological diagnosis confirming endometriosis with clinical stages III–IV according to the revised American Fertility Society (rAFS) classification,⁴⁷ with an average age of 30.50 years (age range 25–36 years). Additionally, four non-endometriotic patients who underwent surgery for uterine fibroids, benign ovarian cysts, etc. were included as controls. These control patients had an average age of 32.25 years (age range 28–37 years). All participants had regular menstrual cycles (21–35 days) and had not received hormone therapy for three months prior to surgery. Specimen collection was performed during the proliferative phase of the menstrual cycle. The medical ethics committee of Union Hospital, Tongji Medical College, Huazhong University of Science and Technology (approval no. 2023–0349), and all patients provided written informed consent before surgery.

The collected tissue specimens were used to extract total RNA using TRIzol reagent (Takara, Japan). The concentration and purity of the extracted total RNA were measured using a Nanodrop 2000 spectrophotometer (Thermo Scientific). A total of 1 µg of RNA was used to synthesize the first strand of cDNA in a 20 µL reaction system using the PrimeScript™ RT kit (Takara, Japan). Real-time quantitative PCR was performed using the SYBR Green Master Mix (Vazyme, China) kit in a 25 µL reaction system. Primer sequences are listed in Table S3, with GAPDH used as an internal reference. The relative gene expression levels were calculated using the 2^(-ΔΔCT) method.

Statistical Analysis

All statistical analyses and presentations were carried out with open-source R (version 4.3.1). The cell proportions between two groups were compared using an unpaired, two-tailed Student’s t-test. To determine statistical significance in gene expression or gene signatures between two groups, an unpaired, two-tailed Wilcoxon rank-sum test with Bonferroni correction was used. A P-value of less than 0.05 indicated statistical significance.

The Wilcoxon rank-sum test was used to compare gene expression between groups due to the non-normal distribution of the data. This non-parametric test is appropriate for comparing independent groups with skewed data or outliers, ensuring robustness. To control for type I errors from multiple comparisons, Bonferroni correction was applied. However, this stringent approach may reduce statistical power and increase the risk of overlooking biologically or clinically relevant findings. Future research could explore alternative methods, such as false discovery rate (FDR) control, to improve sensitivity while maintaining appropriate error control.

Results

Cell Subpopulation Annotation and Contribution Analysis to Endometriosis in scRNA-Seq Data

The overall schematic outline of the current study is illustrated in Figure S1. We downloaded scRNA-seq data files for GSE179640 and GSE213216 from the GEO database. GSE179640 included 6 samples: proliferative endometrium from healthy controls (n=3) and endometriosis patients (n=3). GSE213216 included 10 samples: proliferative endometrium from healthy controls (n=3) and endometriosis patients (n=7). The combined data from both datasets were processed through standardization, normalization, PCA, and Harmony analysis. Initial filtering was based on nFeature_RNA (>6000) and percent.mt (<20) (Figure 1A). Ultimately, in the UMAP analysis, 39 clusters were identified based on the unlabeled gene expression profiles (Figure 1B). To link these clusters to known cell types or functional states, we further annotated these subpopulations using the “SingleR”³¹ package, the CellMarker database,³² and previously reported literature. After both automatic and manual annotation, a total of 39 clusters were classified into eight distinct cell types: B cells/plasma cells, endothelial cells, epithelial cells, mast cells, mesenchymal cells, myeloid cells, smooth muscle cells, and T/NKT cells (Figure 1C and D). The expression patterns of key marker genes across these cell types were visualized using bubble plots (Figure 1E).

Figure 1 Identification of eight cell clusters from core cells in the scRNA-seq dataset and calculation of their contribution to endometriosis. (A). Using nFeature RNA<6, 000 and>200, along with percent.mt<20 as filtering criteria, quality control and filtering were performed on the obtained scRNA-seq data, resulting in a total of 128,550 core cells. (B). Using the UMAP algorithm, the core cells were divided into 39 distinct cell clusters. The cellular annotation profiles of the 39 cell clusters were classified into 8 cell types: B/plasma cells, endothelial cells, epithelial cells, mast cells, mesenchymal cells, myeloid cells, smooth muscle cells, and T/NKT cells. (C). Annotation results for the core cells from the scRNA-seq datasets of endometrial samples from the control and endometriosis groups; (D). Annotation results after merging the core cells from the scRNA-seq datasets of both groups. (E). The expression levels of marker genes for each cell cluster were visualized using bubble plots. (F). The contribution levels of each cell cluster in the eutopic endometrium to the endometriosis were illustrated using a rose chart.

We then analyzed scRNA-seq data from endometrial tissue samples of both healthy controls and endometriosis patients. From these data, we identified the top 100 genes with the highest expression levels in eutopic endometrial tissue (see Table S4 for the detailed list of top 100 highest expression genes in eutopic endometrium vs healthy endometrium). Subsequently, we calculated the differential expression levels (FC) and PctProp of these genes within each cell subtype (refer to Table S4 for the detailed results). The contribution of each cell type to endometriosis was then assessed by calculating the square root of FC* PctProp. Our results indicate that mesenchymal cells have the highest contribution to the pathology of endometriosis (Figure 1F).

Identification of DEGs in Proliferative Phase Endometrial Bulk RNA-Seq Data

Additionally, we obtained the GSE25628 dataset from the GEO public database, which includes bulk RNA-seq data from 13 endometrial tissue samples: 6 from healthy controls and 7 from patients with endometriosis. Differential gene expression analysis between the two groups was conducted using the limma package. Genes were considered differentially expressed if they met the criteria of P.Value < 0.05 and |logFC| > 1. This analysis identified a total of 2,084 differentially expressed genes, with 810 genes upregulated and 1,274 genes downregulated. (Figure 2A, B and Table S5).

Figure 2 The identification of differentially expressed genes (DEGs) between the eutopic endometrium of healthy controls and endometriosis patients. (A). Volcano plot of DEGs between the endometrium of healthy controls and the eutopic endometrium of endometriosis patients in the GSE25628 dataset. P < 0.05 and |log2FoldChange| > 1 were considered significant DEGs. Red dots indicate upregulated genes, while blue dots indicate downregulated genes.(B). Heatmap of DEGs.

Construction and Validation of a Characteristic Gene-Based Predictive Model for Endometriosis

We intersected the DEGs from the bulk RNA-seq dataset GSE25628 with the genes of the mesenchymal cells cluster identified in the scRNA-seq dataset, resulting in 67 intersecting genes (Figure 3A and Box S1). We used the GSE25628 dataset as the training set and the GSE153739 dataset as the validation set, and performed feature selection using Lasso regression. LASSO is a regression technique that allows for many covariates in a model. It uniquely regularizes the absolute values of regression coefficients, potentially shrinking some to zero, thereby enabling variable selection.⁴⁸ The Lasso regression identified 8 genes as characteristic genes for endometriosis, which we selected as candidate genes for further study and used to construct a predictive model (Figure 3B–D). The model formula is as follows: RiskScore = SYNE2 x (−0.119478802623019) + TXN x (−0.0938631926715449) + NUPR1 x (−0.0241295692262499) + CTSK x 0.0227825216179987 + GSN x 0.0244531595971548 + MGP x 0.0595150231025177 + IER2 x 0.134667729309355 + CXCL12 x 0.18988373081461. The results showed that the predictive model constructed with these 8 genes demonstrated good diagnostic performance, with an AUC of 1 (Figure 3E). We further validated the predictive model using the GSE153739 dataset as the validation set, and the results indicated that the model had strong stability, with an AUC of 0.8125 (Figure 3F). Information on the performance metrics of the predictive model is detailed in Table S6.

Figure 3 Construction and validation of a characteristic gene-based predictive model for endometriosis. (A). The DEGs from the Bulk RNA-seq dataset GSE25628 were intersected with the genes of the mesenchymal cell cluster identified in the scRNA-seq dataset, resulting in 67 intersecting genes. (B). Tenfold cross-validation for tuning parameter selection in the LASSO model. (C). Distribution of LASSO coefficients for theDEGs. Positive coefficients (blue) indicate a positive association, while negative coefficients (yellow) indicate a negative association. Bar width represents the magnitude of each gene’s contribution to the predictive model. (D). Histogram of LASSO coefficients for the eight key genes. (E and F). Predictive performance of the training and validation sets. (G). The RT-PCR analysis results of key gene expression levels in endometrial samples from healthy controls and patients with endometriosis. *p < 0.05, **p < 0.01, ns: not significant.

In this study, we further analyzed endometrial tissue from healthy controls and eutopic endometrium from endometriosis patients to assess the expression levels of the 8 key genes identified by our machine learning model. RT-PCR results showed significant upregulation of CXCL12, IER2, MGP, and GSN and downregulation of NUPR1 and TXN in endometriosis. Although the expression levels of CTSK and SYNE2 did not show significant differences, their trends were consistent with model predictions (Figure 3G and Table S7). These RT-PCR results further confirm the reliability of our predictive model. These 8 genes were identified as key candidate genes for our research. Subsequently, we presented the expression of these 8 key genes in 8 different cell types: Mast cells, B/plasma cells, Smooth muscle cells, Mesenchymal cells, T/NKT cells, Myeloid cells, Endothelial cells, and Epithelial cells (Figure 4A and B).

Figure 4 Expression of the 8 key genes across cell clusters. (A). Bubble plots visualizing the expression levels of the 8 key genes across each cell cluster. (B). UMAP plots showing the distribution and expression levels of the 8 key genes (SYNE2, TXN, NUPR1, CTSK, GSN, MGP, IER2, CXCL12) in the scRNA-seq dataset. Colors range from gray to blue, with deeper blue indicating higher expression levels.

Analysis of Immuno-Infiltration in Eutopic Endometria of Endometriosis

The microenvironment, which was primarily composed of immune cells, extracellular matrix, various growth factors, inflammatory factors, and specific physicochemical characteristics, played a significant role in influencing disease diagnosis and clinical treatment sensitivity. Immunologic dysfunction played a crucial role in the implantation of endometriotic tissue. In recent years, several studies suggested that immune dysregulation was also present in the microenvironment of the eutopic endometrium in patients with endometriosis.^49,50 Our study delved into the potential molecular mechanisms of endometriosis by analyzing the relationship between key genes and immune infiltration. In this study, we used both CIBERSORT and ssGSEA algorithms to analyze the distribution of immune cells in each eutopic endometrial sample and the correlations between different immune cell types (Figures 5A, B, and S2A, B). CIBERSORT analysis revealed that, compared to the healthy control group, the levels of CD8 T cells and monocytes were significantly increased in endometriosis patients, while the levels of memory B cells were significantly decreased (Figure 5C). The study also highlighted a strong correlation between key genes and immune cells (Figure 5D–K), indicating that NUPR1, CTSK, GSN, MGP, IER2, and CXCL12 were significantly positively correlated with monocytes, while NUPR1, CTSK, MGP, IER2, and CXCL12 were significantly negatively correlated with memory B cells. Furthermore, the enrichment of 28 immune cell gene sets by ssGSEA further validated the immune cell infiltration distribution in the eutopic endometrium of endometriosis (Figure S2C–K). The study further explored the correlations between these 8 key genes and various immune factors, including immunomodulators, chemokines, and receptors, using data from the TISIDB (Tumor-Immune System Interactions Database)⁵¹ (Figure 5L–P). These findings suggested that these key genes were closely linked to immune cell infiltration and played critical roles within the immune microenvironment.

Figure 5 Analysis of Immuno-infiltration in Eutopic Endometria of Endometriosis. (A). The stacked bar chart shows the relative proportions of different immune cell types in the endometria of healthy controls and endometriosis patients. (B). The correlation heatmap displays the correlations between different immune cell types, with colors ranging from blue (negative correlation) to red (positive correlation). (C). The box plot compares the expression levels of different immune cell types in the endometria of healthy controls and endometriosis patients. (D-K). The correlation dot plot shows the correlation analysis between the eight key genes (TXN, CTSK, MGP, CXCL12, NUPR1,GSN, SYNE2, and IER2 respectively) and the infiltrating immune cell types in the endometrial tissues. (L-P). Correlation analysis of the eight key genes with various chemokines, immunoinhibitors, immunostimulators, MHC, and receptors. Dot color indicates the Pearson correlation coefficient (red for positive, blue for negative), and dot size reflects correlation significance, with larger dots indicating smaller p-values. *p < 0.05, **p < 0.01, ns: not significant.

GSEA Analysis of Specific Signaling Pathways Involved in Key Genes

In this study, we further analyzed the specific signaling pathways associated with the 8 key genes to explore the potential molecular mechanisms by which these genes may influence disease progression. The GSEA results revealed that CTSK is enriched in pathways such as the Calcium signaling pathway and the MAPK signaling pathway (Figure 6A); CXCL12 is enriched in pathways including the Calcium signaling pathway and the cGMP-PKG signaling pathway (Figure 6B); GSN is enriched in pathways such as the cGMP-PKG signaling pathway and the Relaxin signaling pathway (Figure 6C); IER2 is enriched in the cGMP-PKG signaling pathway and the Relaxin signaling pathway (Figure 6D); MGP is enriched in pathways such as the cGMP-PKG signaling pathway and the Relaxin signaling pathway (Figure 6E); NUPR1 is enriched in the Apelin signaling pathway and the cGMP-PKG signaling pathway (Figure 6F); SYNE2 is enriched in pathways such as the Chemokine signaling pathway and the IL-17 signaling pathway (Figure 6G); and TXN is enriched in the Fanconi anemia pathway and the mRNA surveillance pathway (Figure 6H). Detailed information on the GSEA analysis of signaling pathways related to key genes can be found in Table S8.

Figure 6 GSEA and Pathway Analysis of Key Genes. (A-H). Gene Set Enrichment Analysis (GSEA) and pathway analysis for eight key genes (CTSK, CXCL12, GSN, IER2, MGP, NUPR1, SYNE2, TXN). The upper panels show enrichment curves for top associated pathways, highlighting gene enrichment as they rank higher in the dataset. The lower panels display circular diagrams of key pathways, with segment width and color intensity indicating the significance of each pathway in the gene’s function.

GSVA and miRNA Interaction Analysis of Eight Key Genes

GSVA results show that high CTSK expression is associated with enrichment in IL6_JAK_STAT3_SIGNALING and IL2_STAT5_SIGNALING pathways (Figure 7A). High CXCL12 expression is linked to WNT_BETA_CATENIN_SIGNALING and HEDGEHOG_SIGNALING pathways (Figure 7B). High GSN expression enriches IL6_JAK_STAT3_SIGNALING and TGF_BETA_SIGNALING pathways (Figure 7C). High IER2 expression is associated with WNT_BETA_CATENIN_SIGNALING and IL6_JAK_STAT3_SIGNALING pathways (Figure 7D). High MGP expression enriches HEDGEHOG_SIGNALING and WNT_BETA_CATENIN_SIGNALING pathways (Figure 7E). High NUPR1 expression is linked to IL6_JAK_STAT3_SIGNALING and IL2_STAT5_SIGNALING pathways (Figure 7F). High SYNE2 expression enriches IL2_STAT5_SIGNALING and TGF_BETA_SIGNALING pathways (Figure 7G). High TXN expression is associated with MTORC1_SIGNALING and DNA_REPAIR pathways (Figure 7H).

Figure 7 Pathway Activity Comparison via GSVA and miRNA-mRNA Interaction Network for Key Genes. (A-H). Gene Set Variation Analysis (GSVA) comparing pathway activity between high-expression (HExp) and low-expression (LExp) groups for eight key genes (CTSK, CXCL12, GSN, IER2, MGP, NUPR1, SYNE2, TXN). Positive t-values (blue bars) indicate pathways more active in the HExp group, while negative t-values (green bars) indicate pathways more active in the LExp group. (I). miRNA-mRNA interaction network for key genes. The network illustrates the interaction between eight key genes (CTSK, CXCL12, GSN, IER2, MGP, NUPR1, SYNE2, TXN) and their predicted targeting miRNAs. The Orange nodes represent the key genes, while the blue nodes represent the miRNAs.

In addition, using the miRcode database (http://www.mircode.org), we performed reverse prediction on the eight key genes, identifying 83 miRNAs and a total of 253 mRNA-miRNA interactions. These interactions were visualized using Cytoscape (Figure 7I).

Transcription Factor Enrichment and Motif Analysis of Key Genes

In this study, we analyzed eight key genes and found that they are regulated by multiple transcription factors. Using GSEA, we evaluated the enrichment of these transcription factors. Motif-TF annotation and key gene selection results indicated that the motif with the highest normalized enrichment score (NES) was cisbp__M5609, with an NES value of 6.24 (Figure 8A and B). We also identified the motifs enriched in these key genes and their corresponding transcription factors, as detailed in Figure 8C, which illustrates the relationships between these genes and the associated motifs.

Figure 8 Transcription factor (TF) motif enrichment analysis for key genes. (A). Summary table of top enriched TF motifs associated with the key genes, including the normalized enrichment score (NES), area under the curve (AUC), and transcription factors annotated to motifs (TF_highConf). (B). Global mean and standard deviation (SD) estimation plots for the top TF motifs, illustrating the rank distribution of associated genes. (C) Interaction network showing the relationship between the enriched TFs (red hexagons) and the key genes (green circles), highlighting potential regulatory interactions.

Correlation Analysis Between Key Genes and Endometriosis-Related Disease Genes

This study obtained endometriosis-related disease genes through the GeneCards database (https://www.genecards.org/). The Relevance Score in the GeneCards database is used to measure the degree of association between a gene and a specific disease; the higher the score, the stronger the association with the target disease. We selected the top 20 genes with the highest Relevance Scores (Top 20) and analyzed their expression levels in endometriotic and healthy control endometrial tissues, exploring the differences in expression between the two groups (Figure 9A). The results showed that the expression of genes such as CYP17A1, ESR1, ESR2, and GNRH1 differed significantly between the two groups. Additionally, we found that the expression levels of eight key genes were significantly correlated with the expression levels of the Top 20 regulatory genes in endometriosis. Among them, CXCL12 and MMP2 were significantly positively correlated (r=0.963), while CTSK and ESR1 were significantly negatively correlated (r= −0.885) (Figure 9B).

Figure 9 Drug-Target Predictions for Endometriosis. (A). Box plots showing gene expression levels in control (blue) and endometriosis (pink) groups. Significance: *p < 0.05, **p < 0.01, ***p < 0.001. (B). Pearson correlation analysis of key genes with endometriosis-associated genes, with blue indicating negative and red indicating positive correlations. (C-F). The chemical structures of retinol, orantinib, piperacillin, and NECA, respectively. These four drugs were identified as potential treatments for endometriosis through the Connectivity Map (CMap) database.

Drug-Target Predictions Based on Common DEGs

By analyzing the Bulk RNA-seq datasets GSE25628 and GSE153739, we identified common DEGs between the two datasets and selected the top 150 upregulated and downregulated genes based on differential expression. Drug-target predictions for these differential genes were performed using the CMap.⁴⁶ CMap, a transcriptomics-based drug repurposing tool, has been used to identify potential therapeutic candidates for various conditions, including endometriosis. Several studies have demonstrated its capability to identify compounds with therapeutic potential for this disease.^52,53 In our study, the results indicated that the expression profiles perturbed by drugs such as Retinol (Figure 9C), Orantinib (Figure 9D), Piperacillin (Figure 9E), and NECA (Figure 9F) were significantly negatively correlated with the expression profiles perturbed by endometriosis, suggesting that these drugs may have the potential to alleviate or even reverse the disease state. (Table S9).

Discussion

As a prevalent chronic condition, endometriosis can significantly impact the well-being of women and their families.⁵⁴ Despite the lack of comprehensive understanding regarding its pathogenesis, it is postulated that multifactorial interactions contribute to the development of endometriosis, while current knowledge in this field remains limited, particularly regarding immune dysregulation.⁵⁵ Therefore, further research is urgently warranted to comprehensively elucidate the immunopathological mechanisms underlying endometriosis and to investigate potential diagnostic biomarkers and optimal immunotherapy strategies for this condition.

Immune Cell Profile Dysregulation in Eutopic Endometrium of Endometriosis Immune Dysfunction

Immune dysfunction has been identified as a key factor in the implantation of endometriotic tissue. Recent research indicates that immune dysregulation also occurs within the microenvironment of the eutopic endometrium in endometriosis patients.^49,50 Huang X et al⁵⁰ utilized the 10x Genomics platform for scRNA-seq to analyze endometrial cells from endometriosis patients and controls. They identified differences in immune cell proportions and IL-10 secretion during the secretory phase, with elevated proinflammatory cytokines in endometriosis. These findings suggest that immune microenvironment alterations in the eutopic endometrium may impair endometrial receptivity. In our current study, according to CIBERSORT and ssGSEA analysis, we observed a higher abundance of CD8+ T cells and monocytes and an increasing trend in M2 macrophages in patients with endometriosis compared to the control group, which aligns with previous research findings.^56,57 In 2021, Ma J et al performed scRNA-seq on endometrial tissue from three healthy controls and three endometriosis patients.⁵⁷ The study identified significant differences in cell subtypes and gene expression between the eutopic endometrium of endometriosis patients and healthy controls, with monocytes notably more prevalent in the eutopic endometrium of endometriosis. The mononuclear macrophage system is considered a pivotal factor in the promotion of endometriosis, wherein mononuclear cells are recruited to the lesion site and induced to differentiate into macrophages, which serve as the principal immune cells responsible for pro-inflammatory chemokine production.^58,59 Previous studies have reported a higher number of CD8+ T cells in the eutopic endometrium of endometriosis patients compared to normal controls.⁵⁶ Furthermore, increased levels of CD8+ T cells are associated with an elevated risk of infertility among individuals with endometriosis. In patients with endometriosis, CD8 MAIT cells represent the predominant subset in both peritoneal fluid and peripheral blood. These cells have the potential to induce a pro-inflammatory state and create an immunosuppressive microenvironment by secreting IL-17 and initiating a Th17-biased immune response.⁶⁰ Unfortunately, despite compelling evidence of dysregulation in key regulators of CD8+ T cell function such as PD-1, PD-L1, and CD4+ T regulatory cells in endometriosis, functional analyses of CD8+ T cells have not yet been conducted to date.^20,61,62

Identifying Key Genes and Models for Predicting Endometriosis

Historically, the lack of specificity and sensitivity in a single biomarker or group of biomolecules has hindered their utilization as diagnostic tests for endometriosis.⁶³ Using integrated scRNA-Seq and bulk RNA-Seq datasets with machine learning algorithms, we developed a predictive model for endometriosis and identified eight key genes (SYNE2, TXN, NUPR1, CTSK, GSN, MGP, IER2, and CXCL12). These findings were validated through qRT-PCR using clinical specimens. The identified genes were closely associated with immune cell infiltration.

Previous studies have shown that CXCL12 is significantly upregulated in endometriosis lesions compared to healthy endometrial tissue,⁶⁴ consistent with our findings. CXCL12 binds to CXCR4/7, activating multiple signaling pathways,^65,66 inhibiting autophagy in stromal cells, and promoting the proliferation, migration, and invasion of endometriotic cells. GSEA analysis revealed that CXCL12 is closely linked to monocyte infiltration, suggesting it recruits immune cells, causing immune dysregulation in the endometrial microenvironment and accelerating disease progression.

CTSK, a lysosomal cysteine protease, activates the TLR7/9 pathway in antigen-presenting cells, promoting Th17 differentiation.^67,68 CTSK has been proposed as a diagnostic marker for endometriosis and correlates with immune scores.⁶⁹ Our study showed that CTSK expression is elevated in endometriosis lesions, with GSEA analysis indicating its enrichment in the MAPK pathway, which is crucial for implantation, growth, invasion, proliferation, and apoptosis in endometriotic tissue.⁷⁰

We also found that SYNE2 expression is significantly reduced in the eutopic endometrium of endometriosis, which aligns with prior studies.⁷¹ GSEA analysis indicated that SYNE2 is enriched in the chemokine and IL-17 signaling pathways, both key in immune regulation.^72–74

TXN (Thioredoxin), a gene critical for cellular redox balance, plays a central role in immune regulation and oxidative stress.⁷⁵ Despite ongoing debate about its expression in endometriosis,^76,77 TXN may contribute to disease pathogenesis by modulating oxidative stress and immune pathways like NF-κB and MAPK.

MGP is low in normal endometrium but significantly elevated in endometriosis lesions, especially in glandular and endothelial cells.⁷⁸ While no studies have reported on GSN,⁷⁹ IER2, and NUPR1 expression in endometriosis, GSEA analysis shows these genes, along with MGP, are enriched in the cGMP-PKG and Relaxin pathways. The Relaxin pathway may regulate endometriosis by inhibiting IL-8 expression in stromal cells.⁸⁰

Collectively, these genes hold potential as novel diagnostic biomarkers for endometriosis; however, further investigation is warranted to elucidate their precise role in the pathogenesis of this condition. Moreover, our study primarily focused on endometriosis lesions, and additional clinical studies are necessary to determine whether these genes exhibit differential expression in peripheral blood and can serve as non-invasive biomarkers for diagnosis.

Drug Target Prediction for Endometriosis Based on Eutopic Endometrium Sequencing Data

Through drug-target prediction analysis, we identified Retinol, Orantinib, Piperacillin, and NECA as compounds with expression profiles significantly inversely correlated with those associated with endometriosis. These findings suggest their potential therapeutic value in mitigating or reversing the disease state. Retinol and its derivatives, collectively known as retinoids, play a crucial role in endometrial development, maintenance and stromal decidualization, and blastocyst implantation. Dysregulated retinoid metabolism may disrupt mitochondrial function and facilitate the implantation and growth of ectopic cells, potentially serving as a pivotal factor in the pathogenesis of endometriosis.^81–84 Retinoids have found widespread use in the treatment of various skin diseases, malignant tumors, and immune disorders. However, their poor solubility in aqueous solutions, increased catabolism upon intravenous injection, and potential for systemic side effects have prompted researchers to develop novel retinoid delivery agents with improved tolerability.⁸⁵ In the context of endometriosis, it is worth investigating strategies for targeted delivery of retinoids to the lesion site and specific immune cell populations. Orantinib (SU6668) exhibits multipotent inhibition of VEGF, FGF, and PDGF receptor tyrosine kinase activity, effectively suppressing the angiogenesis and vascularization processes in ectopic endometrial tissue.⁸⁶ However, further investigation is required to elucidate its specific regulatory mechanism. The primary adverse events associated with orantinib treatment were edema, ascites, and elevations in levels of aspartate aminotransferase and alanine aminotransferase.⁸⁷

Limitations

Although this study explored the correlation between immune cell infiltration in the endometrium and endometriosis occurrence, and developed a predictive model using bioinformatics based on identified key genes, there are still some limitations. First, this study analyzed sequencing data only from the proliferative phase of the endometrium, without investigating dynamic changes in the immune microenvironment and transcriptomic expression across the entire menstrual cycle. This limitation may restrict the ability to fully capture the molecular characteristics of the eutopic endometrium in endometriosis. Longitudinal studies are warranted to explore these dynamic changes and provide a more comprehensive understanding. Second, the relatively small sample size, particularly in the validation set (n = 7), and the lack of detailed clinical follow-up data may limit the generalizability and predictive accuracy of the model. Furthermore, the influence of comorbidities or other health conditions was not fully considered, which could impact biomarker identification and model stability. Future studies should include larger sample sizes and more comprehensive clinical datasets to improve model reliability and applicability. Third, the reliance on in silico analyses introduces methodological constraints, including potential biases in publicly available sequencing datasets and the inherent limitations of computational predictions. While in vitro validation using qPCR provided initial support for the findings, the small number of clinical samples analyzed may reduce the generalizability of these results. Expanding the validation cohort to include larger clinical sample sizes will be crucial for confirming these findings and enhancing their broader applicability. Lastly, drug-target predictions in this study identified potential therapeutic agents, such as Retinol and Orantinib, with expression profiles negatively correlated with those associated with endometriosis. However, their efficacy and underlying mechanisms remain to be validated in vivo. Rigorous experimental studies utilizing endometriosis models are essential to confirm their therapeutic potential and elucidate their mechanisms of action. Addressing these limitations in future research will help strengthen the study’s findings, advance understanding of endometriosis pathophysiology, and support the clinical translation of these insights.

Conclusion

This study elucidates the role of immune dysregulation in the pathogenesis of endometriosis and identifies potential diagnostic biomarkers (SYNE2, TXN, NUPR1, CTSK, GSN, MGP, IER2, CXCL12) as well as candidate drugs for this condition. These findings offer valuable insights for the development of early detection methods and targeted therapeutic strategies for endometriosis.