Identification and Validation of a Four-Gene Signature as a Novel Potential Diagnostic Biomarker Panel for Endometriosis: Causal and Mechanistic Insights from Mendelian Randomization and Bioinformatics

Endometriosis (EMS) is an estrogen-dependent disease that can lead to chronic pain and infertility. Small ubiquitin-related modifier (SUMO)ylation, a key post-translational modification, has been demonstrated to be associated with endometrial dysfunction, but its role in EMS remains unclear. This study aims to identify SUMOylation-related genes as potential biomarkers for EMS utilizing bioinformatics approaches.

The GSE51981, GSE7305, and GSE214411 datasets and 194 SUMOylation-related genes were assessed. The differential expression analysis and weighted gene co-expression network analysis were performed to identify candidate genes associated with EMS and SUMOylation. Subsequently, the univariable and multivariable Mendelian randomization analyses were used to investigate causal relationships. The receiver operating characteristic curve and nomogram were constructed for assessing the performance of candidate genes. Furthermore, the biological roles of candidate biomarkers were evaluated through gene set enrichment analysis. Additionally, immune infiltration levels were analyzed and pseudo-time trajectory analysis was conducted to assess cell differentiation dynamics.

CRYZ, FRMD4B, FZD6, and TPST1 were identified as biomarkers associated with SUMOylation in EMS, all with area under the curve values above 0.75, suggesting good diagnostic performance. The constructed nomogram, which incorporated these biomarkers, demonstrated exceptional predictive capabilities for EMS risk. Furthermore, these biomarkers functioned through common enriched pathways, including the cell cycle, neuroactive ligand-receptor interaction, and ubiquitin-mediated proteolysis. Notably, immune cell analysis revealed significant infiltration of regulatory T cells and memory B cells in EMS. Pseudo-time trajectory analysis indicated that epithelial cells displayed distinct expression patterns for the biomarkers across five differentiation stages.

This study identifies CRYZ, FRMD4B, FZD6, and TPST1 as potential diagnostic and therapeutic targets associated with SUMOylation in EMS. These findings provide potential clinical insights for the non-invasive detection and risk stratification of EMS, which might help improve the efficiency of early diagnosis and the development of personalized treatment strategies.

endometriosis, small ubiquitin-like modifier ylation, biomarkers, Mendelian randomization, bioinformatics analysis

Endometriosis (EMS) is a prevalent gynecological disorder classified as an estrogen-dependent disease and characterized by the presence of normal endometrial tissue outside the uterine cavity, severely affecting women’s physical and mental well-being.1 Common sites for ectopic endometrial tissues include the ovaries, sacro-ligament, rectouterine pouch, pelvic peritoneum, as well as the rectum, bladder, and nerves. Approximately 10% of women of reproductive age are affected by EMS.2 Patients may experience complications such as infertility, chronic pelvic pain, and dysmenorrhea, with varying clinical manifestations depending on the ectopic locations of the endometrial tissue.2 For example, EMS involving the ovaries can lead to the formation of ovarian cysts, which may rupture and result in acute gynecological abdominal symptoms.3 When EMS affects the sacro-ligament and rectouterine pouch, patients may experience dyspareunia, as well as palpable nodules in the sacro-ligament and tenderness in the rectouterine pouch.4,5 Endometrial tissue located in the rectum can lead to intestinal masses, resulting in symptoms such as intestinal perforation.6 When nerves are involved, localized pain may manifest.2,7 Currently, laparoscopy remains the gold standard for diagnosing EMS; however, it is invasive, costly, and operator-dependent, making it difficult to use for large-scale screening and early diagnosis.8 Treatment strategies typically involve surgical intervention and conservative pharmacological management. Unfortunately, approximately half of patients experience recurrence within five years post-treatment, coupled with a potential risk for malignant transformation into ovarian endometrioid adenocarcinoma.9,10 Therefore, exploring non-invasive or minimally invasive biomarkers associated with EMS is crucial for early diagnosis, enhancing the understanding of its etiology, and refining treatment strategies, highlighting the urgent need for developing reliable biomarkers in the near future. Small ubiquitin-related modifier (SUMO)ylation is a post-translational modification that serves as a crucial molecular regulatory mechanism in cell death processes, including apoptosis, autophagy, and senescence.11 It has been shown that SUMOylation can stabilize the transcription factor (TF) 21.12 The SUMOylation of TF21 may enhance its interaction with upstream stimulatory factor 2, thereby increasing the binding ability of upstream stimulatory factor 2 with the SF-1 and estrogen receptor beta promoters and influencing the proliferation of embryonic stem cells, which contributes to the development of EMS. However, research on SUMOylation in the context of EMS remains limited, and its mechanisms are not yet fully understood. The causal relationship between key SUMOylation-related genes (SRGs) and EMS remains unclear. Moreover, the roles of key SRGs in this condition remain unclear. This study primarily utilizes single-cell RNA sequencing (scRNA-seq) data and conventional transcriptomic data, employing various bioinformatics approaches, including Mendelian randomization (MR) analysis, differential expression analysis, and weighted gene co-expression network analysis (WGCNA), to identify SUMOylation-related biomarkers in EMS. Additionally, gene set enrichment analysis (GSEA) was performed on these biomarkers, and a long non-coding RNA (lncRNA)-microRNA (miRNA)-biomarker networks and TF-biomarker-miRNA networks were constructed. Immune infiltration analysis was conducted to assess the levels of immune cell infiltration in EMS. Furthermore, pseudo-time analysis and cell communication analysis were based on identified key cell types to elucidate their significance in EMS. Our findings may provide new insights for the diagnosis and treatment of patients with EMS.

Data Source The common transcriptome sequencing datasets related to EMS, namely GSE51981 and GSE7305, along with the scRNA-seq dataset GSE214411, were obtained from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/). The GSE51981 dataset included endometrial tissue samples from 77 EMS patients (including 28 with minimal/mild [early-stage] and 49 with moderate/severe [late-stage] disease) and 34 healthy controls (Supplementary Table S1). Similarly, for expression validation, the GSE7305 dataset contained endometrial tissue samples from 10 patients with EMS and 10 healthy controls. The GSE214411 dataset comprised 6 EMS endometrial tissue samples and 7 normal control tissue samples. A total of 194 SRGs were sourced from the Molecular Signatures Database (https://www.gsea-msigdb.org/).13,14 scRNA-Seq Data Pre-Processing and Cell Annotation First, the raw data were normalized using the “PercentageFeatureSet” function. Both the gene/feature counts (nCount/nFeature_RNA) and the percentages of mitochondrial genes across EMS and controls within the GSE214411 dataset were calculated. To ensure the selection of high-quality cells, quality control was performed by creating a “Seurat” object from the scRNA-seq data in the GSE214411 dataset using the R package “Seurat” (Version 5.0.1),15 with gene counts between 200 and 10,000, a total gene expression sum in each cell of less than 50,000, and a mitochondrial gene proportion of less than 20%. These thresholds were selected based on the biological characteristics of endometrial tissue and are consistent with previous studies.16–18 Specifically, the mitochondrial gene proportion threshold (<20%) accounts for the high metabolic activity of endometrial cells and potential cellular stress during single-cell suspension preparation; the gene count range (200–10,000) accommodates cell type diversity and pathological heterogeneity while removing low-quality cells and potential doublets; and the total expression threshold (<50,000) is primarily used to exclude doublet cells. Subsequently, the “LogNormalize” function was applied to normalize feature expression values in each cell by scaling them based on total expression. The “FindVariableFeatures” function was used to filter the data, selecting the top 2,000 highly variable genes following quality control. Furthermore, the principal components were obtained after dimensionality reduction and clustering using the “RunPCA” and “RunUMAP” functions. To identify the optimal dimensions for cell clustering, the scatter plots for the dimensions of linear dimensionality reduction were generated using the “ElbowPlot” function. UMAP was used to delineate cell clusters based on a latitude value of 50 using the “RunMap” function. An unsupervised cluster analysis was conducted on all cells in the GSE214411 dataset using the “FindNeighbors” and “FindClusters” functions to obtain cell clusters (resolution = 0.4). Finally, marker genes for each cell cluster were identified using the “FindAllMarkers” function, applying a filter of |Log2 fold-change| > 0.5. Based on these marker genes, the cell clusters were annotated to specific cell types in both EMS and control samples within the GSE214411 dataset using the R package “SingleR”. Additionally, the relative abundance of each cell type in EMS and controls was presented in a stacked bar chart. Differences in cell type abundance between EMS and controls were assessed using the Wilcoxon rank-sum test. Gene Set Variation Analysis (GSVA) GSVA was conducted to investigate the variation in biological functions between EMS and controls within the GSE214411 dataset. Specifically, the “msigdbr” package (Version 7.5.1) was used to load the “hallmark gene sets”, subsequently assigning pathway activity scores to individual cells. Following this, the “limma” package in R (Version 3.54.0)19 was utilized to calculate the discrepancies in pathway activity scores across each cell type. Identification of Differentially Expressed Genes (DEGs) in EMS at the Single-Cell Level To identify DEGs between EMS and controls at the single-cell level, the gene expression profiles within specific cell types in the GSE214411 dataset were analyzed using the “FindMarkers” function. Subsequently, likelihood ratio tests were conducted to detect DEGs. The p-values were adjusted using the Benjamini-Hochberg method for multiple testing, and significant DEGs were identified at the single-cell level with a statistical threshold of false discovery rate 1 and p-value < 0.05. The intersecting DEGs from the GSE214411 dataset and the GSE51981 dataset were considered common DEGs. Analysis of Key Module Genes Associated with SUMOylation in EMS Using WGCNA The single-sample GSEA algorithm within the R package “GSVA” was first used to calculate SUMOylation scores for EMS and controls in the GSE51981 dataset. The differences between EMS and controls were compared using the Wilcoxon rank-sum test. Subsequently, WGCNA was conducted using the R package “WGCNA” (Version 1.70.3).21 Outlier samples in the GSE51981 dataset were identified and removed through cluster analysis. The optimal soft threshold power (β) was determined when the scale-free topological fitting index exceeded 0.8 and the mean connectivity approached 0, ensuring that the constructed network conformed to scale-free distribution principles. After that, adjacency and similarity among genes were calculated to infer dissimilarity coefficients, facilitating the construction of the co-expression network according to the dynamic tree-cutting algorithm. The minimum number of genes per module was set to 100 for merging clustering modules. Using the SUMOylation score as the phenotype, Spearman correlation analysis was conducted to assess the relationship between each module and the scores of PRGs (Prognostic Risk Genes Score) and RMCRGs (RNA Modification Consensus Prognostic Risk Score). The clustering module exhibiting the highest correlation (|R| > 0.3 and p-value 0.8 and |Gene Significance| > 0.2. The overlap between DEGs and key module genes, analyzed using the R package “VennDiagram” (Version 1.7.3),22 was defined as candidate genes. Functional Analysis and Protein-Protein Interaction Evaluation of Candidate Genes To investigate the biological functions of candidate genes, functional annotation analyses, including Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG), were performed using the R package “clusterProfiler” (Version 4.7.1.003).23 To explore the protein-level interactions of candidate genes, they were uploaded to the Search Tool for the Retrieval of Interacting Genes database (STRING, http://www.string-db.org/) for constructing a protein-protein interaction network (interaction score > 0.15). The results were then imported into Cytoscape software (Version 3.7.2)24 for visualization. Screening of Instrumental Variables (IVs) Genome-wide association studies (GWAS) data for expression quantitative trait loci of candidate genes and the trait ID for the outcome of EMS (ebi-a-GCST90018839) were obtained from the IEU Open GWAS database (https://gwas.mrcieu.ac.uk/). A total of 24,089,752 single-nucleotide polymorphisms (SNPs) from 231,771 samples (4,511 patients with EMS and 227,260 healthy controls) were analyzed. These data were derived from the European population. The IVs were selected based on three fundamental assumptions: Assumption 1, IVs must have a strong and consistent association with the candidate genes; Assumption 2, the causal relationship between the candidate genes and EMS must be independent of confounding factors; Assumption 3, IVs should influence the risk of EMS directly through the candidate genes without affecting other pathways. Specifically, exposure reading and SNP filtering were performed using the “extract instruments” function in the R package “TwoSampleMR” (Version 0.5.8).25 With p < 5 × 10−8, SNPs that were significantly correlated with candidate genes were identified, employing parameters clump = TRUE, r2 = 0.001, and kb = 10, to remove SNPs with linkage disequilibrium. Next, the “extract_outcome_data” function was used to read the outcome SNPs in conjunction with those corresponding to the candidate genes, while IVs unrelated to the outcome were filtered (proxies = TRUE; rsq = 0.8). The exposure and outcome data were then combined using the “harmonise_data” function to prepare data for MR analysis. This step harmonized effect allele directions between exposure and outcome datasets to avoid sign errors in causal estimates. Then, a Steiger directionality test was conducted to verify the direction of causality. This test confirmed the causal direction by ensuring the instrumental variables explained more variance in the exposure than in the outcome (r2.exposure > r2.outcome, Steiger test p < 0.05); otherwise, reverse causality could not be excluded. Furthermore, to assess the potential bias of weak instrumental variables on the causal estimation, we calculated the F-statistic for each SNP (). SNPs with F-values 10, indicating no bias from weak instrumental variables (Supplementary Table S2). If fewer than three SNPs were available for candidate genes, those genes were excluded from subsequent analyses. Handling of Potential Biases in MR Analysis To address potential false positives, the MR analysis was hypothesis-driven based on 49 pre-selected candidate genes; thus, no strict multiple comparison correction was applied. Population bias was minimized as both exposure and outcome GWAS data were derived from European populations. Allele mismatch was resolved using the harmonise_data function to align effect alleles and remove palindromic SNPs. Univariable MR (UVMR) Analysis UVMR analysis was conducted to elucidate the causal relationships between candidate genes and EMS. The effect alleles and effect sizes were harmonized using the “harmonise_data” function. To perform the UVMR analysis, the “mr function” was applied in combination with five different algorithms, including MR Egger,26 weighted median,27 simple mode,28 inverse variance weighted (IVW) method,29 and weighted mode.30 The outcomes of the IVW method were predominantly selected among the five algorithms evaluated. Risk factors for EMS were identified when candidate exposure factors had a p-value below 0.05 and a beta-value greater than 0, while protective factors were recognized with a beta-value below 0. The scatter plots, forest plots, and funnel plots were used for visualizing the results of the UVMR analysis. Multivariable MR (MVMR) Analysis Based on the candidate exposure factors identified through the UVMR analysis, SNPs significantly associated with multiple exposures were extracted using the “mv_extract_exposures” function, with a p-value threshold of 5×10−8. This was followed by the removal of SNPs exhibiting linkage disequilibrium, applying the parameters of clump = TRUE, clump_r2 = 0.001, and clump_kb = 10. SNPs significantly associated with EMS were excluded using the “textract_outcome_data” function (proxies = TRUE; rsq = 0.8). The effect alleles and effect sizes were standardized using the “mv_harmonise_data” function. Subsequently, the “mv_lasso_feature_selection” function was employed to eliminate collinear screening variables. MVMR analysis was conducted using the “mv_multiple” function with five algorithms. Ultimately, candidate key genes were identified based on their odds ratio (OR) values. Genes with an OR greater than 1 and upregulated in EMS, as well as those with an OR less than 1 and downregulated in EMS, were selected as candidate key genes. The reliability and robustness of the MR analysis were assessed using sensitivity analyses. In detail, the heterogeneity test was conducted using the “mr_heterogeneity” function. The IVW and MR-Egger regression methods were used to determine the presence of heterogeneity by calculating Cochran’s Q statistic. Typically, a p-value greater than 0.05 indicates the absence of heterogeneity. Furthermore, the “mr_pleiotropy_test” and “mr_presso” functions were used to examine horizontal pleiotropy among SNPs, where p-values greater than 0.05 suggest an absence of horizontal pleiotropy. Additionally, the leave-one-out test was performed using the “mr_leaveoneout” function to evaluate whether the MR results were influenced by any individual SNP. Nomogram Modeling and Performance Evaluation The receiver operating characteristic (ROC) curves were plotted using the R package “pROC” (version 1.0–11)31 to explore the diagnostic value of candidate key genes in EMS, including their performance for early-stage (minimal/mild) and late-stage (moderate/severe) disease based on the GSE51981 dataset. Moreover, the expression patterns of candidate key genes in EMS and controls from the GSE51981 and GSE7305 datasets were examined, with candidate key genes showing disparate expression levels in EMS and consistent expression trends across the two datasets being identified as biomarkers (p-value < 0.05). To investigate the specific roles of these biomarkers in diagnosing EMS, a nomogram was constructed based on the biomarkers from the GSE51981 dataset using the “rms” package in R (version 6.5.0).32 Additionally, a calibration curve was generated using the “rms” package to assess the accuracy of the nomogram model predictions, along with a decision curve analysis. The diagnostic performance of the nomogram was tested in the GSE51981 dataset and externally validated in the independent GSE7305 dataset, with its ROC curve plotted using the “pROC” package. Exploration of Biological Pathways and Regulatory Mechanisms Associated with Biomarkers in EMS To investigate the biological pathways associated with biomarkers in EMS, GSEA was conducted for each biomarker in the GSE51981 dataset using the R package “clusterProfiler” and the KEGG gene set (c2.cp.kegg.v7.5.1.symbols.gmt). Specifically, Spearman correlation analysis was performed between each biomarker and other genes, with the correlation coefficient serving as a ranking criterion, and the ranked genes were selected for GSEA (False Discovery Rate < 0.05). Furthermore, to explore interactions between biomarkers and other genes sharing similar functions, a gene-gene interaction network was developed using the GeneMANIA database (https://genemania.org/). To clarify the regulatory mechanisms of biomarkers in EMS, an lncRNA-miRNA-biomarker network and a TF-biomarker-miRNA network were constructed. First, crucial miRNAs targeting biomarkers were identified by intersecting miRNAs from the Miranda database with those from the MicroCosm database. Subsequently, lncRNAs targeting these crucial miRNAs were predicted using the TSTARBASE database (https://www.mirnet.ca). Simultaneously, biomarker-associated TFs were predicted in the NetworkAnalyst database (https://www.networkanalyst.ca/). Finally, both the lncRNA-miRNA-biomarker network and the TF-biomarker-miRNA network were visualized using Cytoscape software. Assessment of Immune Cell Infiltration and Identification of Potential Therapeutic Targets in EMS To assess the level of immune cell infiltration during the development of EMS, the infiltration scores of 22 immune cell types in each sample from the GSE51981 dataset (p-value < 0.05) were calculated using the CIBERSORT (version 0.1.0).33,34 The correlations among the 22 immune cells were estimated using Spearman correlation analysis. The Wilcoxon rank-sum test was utilized to compare the infiltration levels of the 22 immune cells between the EMS and control samples. Furthermore, the correlations between differing immune cell types and biomarkers were evaluated through Spearman correlation analysis. Additionally, potential therapeutic drugs targeting the biomarkers were searched for in the Drug Signatures Database (https://ngdc.cncb.ac.cn/databasecommons/database/id/4603), and a drug-biomarker network was subsequently created using Cytoscape software. Cell Communication and Pseudo-Temporal Trajectory Analyses of Key Cell Types Box plots were plotted to display the expression of biomarkers in various cell types from the GSE214411 dataset. Cell types exhibiting significant differences (p-value < 0.05) in biomarker expression were identified as key cell types. To investigate the key differentiation stages critical to the process of EMS, the pseudo-temporal trajectory analysis of these key cell types was performed using the “Monocle2” package (version 2.26.0).35 Furthermore, cell communication networks were constructed to understand the communication between key cell types and other cell types using the “CellChat” package (version 1.6.1).36 Statistical Analysis The R software (version 4.2.3) was employed for bioinformatics analysis. The Wilcoxon rank-sum test was conducted to identify inter-group differences, with a p-value of less than 0.05 considered statistically significant.

Cellular Composition and Activated Pathways in EMS Based on scRNA-Seq Data We first examined the cellular composition and identified activated biological pathways associated with EMS based on the scRNA-seq dataset. The violin plots in the Supplementary Figure S1A and S1B illustrate the number of nFeature_RNA, nCount_RNA, and the percentage of mitochondrial genes both before and after quality control. The top 2,000 highly variable genes were highlighted with red dots (Figure 1A). The top 50 principal components were selected for UMAP dimensionality reduction (Figure 1B). Following UMAP analysis and cell annotation, 30 cell clusters were merged into 9 distinct cell types based on marker gene expression, including endothelial cells, epithelial cells, fibroblasts, monocytes, mesenchymal stem cells, neutrophils, natural killer cells, smooth muscle cells, and tissue stem cells (Figure 1C and D). The top 5 marker genes for each cell cluster effectively distinguished the various cell types (Supplementary Figure S1C). In terms of cell type proportions, fibroblasts, tissue stem cells, natural killer cells, and epithelial cells constituted a substantial portion of the nine cell types present in EMS (Figure 1E). Notably, epithelial cells and fibroblasts displayed distinct characteristics, as epithelial cells were more abundant while fibroblasts were relatively sparse in EMS (Figure 1F). Furthermore, GSVA revealed that pathways co-enriched by epithelial cells and fibroblasts, such as Notch signaling, angiogenesis, and epithelial-mesenchymal transition, were significantly activated in EMS (Figure 1G and H). Identification and Functional Analysis of 49 Candidate Genes in EMS In the GSE214411 dataset, a total of 2,482 DEGs were identified, with 1,246 exhibiting upregulation and 1,236 showing downregulation in EMS (Figure 2A and Supplementary Figure S2A). In the GSE51981 dataset, 4,959 DEGs were identified between EMS and control samples, comprising 1,883 upregulated genes and 3,076 downregulated genes (Figure 2B and Supplementary Figure S2B). An overlap analysis revealed 566 common DEGs across both datasets (Supplementary Figure S2C). The SUMOylation score was significantly lower in EMS compared to controls (p = 3.5×10−7) (Supplementary Figure S2D). Cluster analysis indicated the absence of outlier samples in the GSE51981 dataset (Figure 2C). Additionally, the parameter β was determined to be 16 when R2 approached 0.806, with mean connectivity close to 0 (Figure 2D). Furthermore, a co-expression network was constructed based on systematic clustering criteria, resulting in four clustering modules (Figure 2E). Among these, the MEturquoise module was identified as the key module due to its highest positive correlation with the SUMOylation score (R = 0.8, p = 5.3×10−26) (Figure 2F). Subsequently, using |Module Membership| > 0.8 and |Gene Significance| > 0.2 as thresholds, 2,871 key module genes were identified (Figure 2G). These genes overlapped with the 2,482 DEGs, yielding 49 candidate genes (Supplementary Figure S2E). The candidate genes were found to be associated with 108 Gene Ontology entries, including 67 terms in biological processes, 17 in cellular components, and 24 in molecular functions. The specific terms included signal transduction of p53, RNA polymerase II specificity, and protein dephosphorylation (Figure 2H). Regarding the KEGG enrichment analysis of these candidate genes, they were associated with 39 pathways, including protein processing in the endoplasmic reticulum, cholesterol metabolism, and glutathione metabolism (Figure 2I). Overall, the functions of these candidate genes were linked to the regulation of cell signaling, maintenance of metabolic balance, and response to stress. The protein-protein interaction network revealed the interactions among 42 candidate genes, including FRMD4B-ZKSCAN4, Protein tyrosine sulfotransferase 1 (TPST1)-TMEM59L, Frizzled 6 (FZD6)-DSG2, CRYZ-PPA1, and others (Figure 2J). CRYZ (Zeta-Crystallin), FRMD4B (FERM Domain Containing 4B Gene), FZD6 (Frizzled Class Receptor 6), TPST1 (Protein-Tyrosine Sulfotransferase 1), and ZNF606 (Zinc Finger Protein 606) are Identified as Key Protective Factors for EMS The screening of SNPs yielded 40 candidate genes for UVMR analysis. Among the seven candidate exposure factors examined, five exhibited a beta-value below 0, indicating protective factors for EMS (FRMD4B, CRYZ, FZD6, ZNF606, and TPST1), while two showed a beta-value greater than 0, signifying risk factors for EMS (PLK2 and PPA1) (Table 1). According to the IVW method, FRMD4B, CRYZ, FZD6, ZNF606, and TPST1 were protective factors for EMS (slope 0), with the results largely unaffected by confounding effects (Figure 3A). The MR effect sizes for FRMD4B, CRYZ, FZD6, ZNF606, and TPST1 were less than 0 in the forest plot, indicating their potential to reduce the risk of EMS, while PLK2 and PPA1 could increase the risk of EMS due to their MR effect sizes of greater than 0 (Figure 3B). The evenly distributed points on the funnel plots indicated that the UVMR analysis adhered to Mendel’s second law (Figure 3C). The heterogeneity test confirmed no significant heterogeneity in the UVMR analysis, as indicated by p-values greater than 0.05 (Table 2). Additionally, the horizontal pleiotropy test demonstrated an absence of horizontal pleiotropy among SNPs (Tables 3 and 4). Furthermore, in the leave-one-out test, sequentially eliminating SNPs revealed no SNPs sensitive to EMS, suggesting that no single SNP significantly influenced causality (Figure 3D). Overall, the sensitivity tests verified the reliability and robustness of the UVMR results. | Table 1 Univariate MR Analysis of Exposure Factors | | Table 2 Results of the Heterogeneity Test | | Table 3 Results of the Horizontal Pleiotropy Test | | Table 4 Results of Horizontal Pleiotropy Test (Presso) | To investigate the causality of candidate exposure factors at a multivariate level, MVMR analysis was performed. The results emphasized that, in addition to PLK2 as a risk factor for EMS (OR > 1), the other factors acted as protective elements for EMS (OR < 1) (Table 5 and Figure 3E). CRYZ, FRMD4B, FZD6, TPST1, and ZNF606 were identified as key candidate genes due to their risk/protective characteristics aligning with their patterns of upregulation/downregulation. | Table 5 Multivariate MR | Nomogram Constructed Based on CRYZ, FRMD4B, FZD6, and TPST1 Has Good Predictive Capacity for the Risk of EMS Through ROC analysis, we observed that the area under the curve (AUC) for all five candidate key genes was greater than 0.7 in the GSE51981 (Figure 4A) (with CRYZ at 0.791, FRMD4B at 0.798, FZD6 at 0.824, TPST1 at 0.824, and ZNF606 at 0.869) and GSE7305 (Supplementary Figure S3) datasets, indicating high diagnostic accuracy for EMS. Additionally, their expression patterns were analyzed in the GSE51981 and GSE7305 datasets. CRYZ, FRMD4B, FZD6, and TPST1 were significantly downregulated in EMS samples in both datasets (p < 0.05) (Figure 4B and C). In contrast, ZNF606 was significantly downregulated in GSE51981 (p 0.05), indicating a lack of reproducibility across datasets. Therefore, ZNF606 was excluded from further analysis, and CRYZ, FRMD4B, FZD6, and TPST1 were retained as candidate biomarkers. Multivariate logistic regression analysis showed that TPST1 and FZD6 had larger absolute coefficients (TPST1: −1.220; FZD6: −0.981) compared to CRYZ (−0.170) and FRMD4B (−0.194), indicating that TPST1 and FZD6 carry higher weights in the predictive model (Supplementary Table S3). Based on the expression levels and regression coefficients of these biomarkers, a nomogram was constructed (Figure 4D). The total score was calculated by summing the individual scores for each biomarker. A higher total score indicated an increased risk of EMS. Importantly, the slope of the calibration curve was close to 1, demonstrating that the nomogram’s predictive capacity was highly accurate (Hosmer-Lemeshow p-value = 0.23) (Figure 4E). Decision curve analysis revealed that the net benefit of the nomogram was higher than that of the individual factors, emphasizing its superior predictive capacity (Figure 4F). Furthermore, the AUC value of this model was 0.877, indicating that the model has a good discriminative ability (Figure 4G). Finally, we conducted an external validation of the nomogram model using the independent dataset GSE7305. The results showed that the AUC was 0.940, and the DCA also indicated its favorable clinical net benefit (Supplementary Figure S5). However, the sample size of this dataset was small, and further validation in larger external cohorts is warranted in the future. We further evaluated the diagnostic performance of these four biomarkers for different disease stages. As shown in Supplementary Figure S4A, the four biomarkers demonstrated significant diagnostic potential in differentiating late-stage EMS from normal tissues (with AUC values all greater than 0.8), but showed poor diagnostic performance for early-stage EMS. Additionally, as shown in Supplementary Figure S4B, all four genes were significantly downregulated in late-stage EMS compared to both early-stage and control groups (p < 0.05 to p 0.05). The Function of Biomarkers Might Involve Coordinating Gene Expression, Protein Metabolism, and Cell Function GSEA was conducted to elucidate the functions of the biomarkers. The four biomarkers (CRYZ, FRMD4B, FZD6, and TPST1) were enriched in identical pathways, including the spliceosome, neuroactive ligand-receptor interaction, ubiquitin-mediated proteolysis, the cell cycle, protein export, RNA degradation, and oocyte meiosis (Figure 5A–D). Consequently, the functions of these biomarkers may involve the coordination of gene expression, protein metabolism, and cellular function. The gene-gene interaction network revealed 20 additional genes with functional similarities to the four biomarkers, resulting in a total of 155 interactions. For example, FZD6 was linked to SFRP1 through the non-canonical Wnt signaling pathway (Figure 5E). A total of six crucial miRNAs targeting the biomarkers were predicted by intersecting 115 miRNAs from the Miranda database with 82 from the MicroCosm database, resulting in 17 miRNA-mRNA relationship pairs, including two biomarkers (TPST1 and CRYZ) and six crucial miRNAs (hsa-let-7g-3p, hsa-miR-302a-5p, hsa-miR-549a-3p, hsa-miR-23a-3p, hsa-miR-23b-3p, and hsa-miR-656-3p). In total, 204 miRNA-lncRNA relationship pairs were obtained. Finally, a lncRNA-miRNA-biomarker network comprising one biomarker, three miRNAs, and 25 lncRNAs was constructed (Supplementary Figure S6A). The transcription of TPST1 was simultaneously regulated by hsa-miR-23a-3p, hsa-miR-656-3p, and hsa-miR-23b-3p (Supplementary Figure S6A). In the NetworkAnalyst database, 19 TFs were found to target CRYZ, while 5 TFs targeted TPST1. The TF-biomarker-miRNA network revealed that Zinc Finger Protein 2 (ZFP2) could simultaneously regulate the expression of CRYZ and TPST1 (Supplementary Figure S6B). Regulatory T Cells (Tregs), Monocytes, Resting Mast Cells, and Memory B Cells Might Each Play Significant Roles in EMS Figure 6A presents the relative percentages of 22 immune cell types. Correlation analysis among these immune cells indicated that memory B cells and Tregs exhibited a negative correlation with activated memory CD4 T cells, while a positive correlation was observed between Tregs and memory B cells, as well as between memory B cells and plasma cells (p-value < 0.001) (Figure 6B). Memory B cells, monocytes, T follicular helper cells, and Tregs demonstrated high infiltration levels in EMS; however, activated dendritic cells and resting mast cells displayed lower infiltration levels in EMS (Figure 6C). Additionally, Tregs, monocytes, and memory B cells showed a negative correlation with all four biomarkers, whereas resting mast cells exhibited a positive correlation with all four biomarkers (Figure 6D). The Drug Signatures Database provided a total of 143 drugs targeting these biomarkers (Figure 6E). For instance, drugs targeting CRYZ, FRMD4B, FZD6, and TPST1 included clindamycin, medrysone, daunorubicin, and rimexolone, respectively. Key Role of Epithelial Cells in Biomarker Expression and Cell Communication Epithelial cells were identified as key cell types due to the significant differences in the expression levels of the four biomarkers within these cells (Figure 7A). Pseudo-temporal trajectory analysis revealed that the differentiation of epithelial cells could be categorized into five stages (Stage 1 to Stage 5), with Stage 5 marking the onset of differentiation and Stage 1 representing the endpoint (Figure 7B–D). The expressions of the four biomarkers exhibited substantial changes throughout the differentiation process. Specifically, the expression of CRYZ initially increased and then decreased, while FRMD4B demonstrated a similar pattern with an initial rise followed by a decline. FZD6 expression decreased during the later stages of epithelial cell differentiation, whereas TPST1 expression exhibited a downward trend (Figure 7E). Furthermore, cell communication analysis indicated that the frequency and intensity of interactions between epithelial cells and endothelial cells were significantly more pronounced (Figure 7F–H).

In this study, we identified 49 SUMOylation-related DEGs in EMS through integrated analysis of single-cell and transcriptomic data. To mitigate potential confounding factors and reverse causality, MR analysis was further conducted, ultimately revealing four key genes—CRYZ, FRMD4B, FZD6, and TPST1—as significantly associated with EMS. To our knowledge, there have been no previous reports linking these four genes to EMS, making this study the first to establish their association with the disease. GeneMANIA analysis did not indicate strong evidence of direct interactions among these genes. ROC analysis demonstrated that FZD6 and TPST1 exhibited AUC values exceeding 0.8, indicating high diagnostic accuracy, while CRYZ and FRMD4B both showed AUC values above 0.7, suggesting considerable diagnostic value. A nomogram model constructed based on these four genes effectively predicted the incidence of EMS, with the calibration curve slope approaching 1, indicating high predictive accuracy. Indirect evidence suggests that the four biomarkers may be biologically relevant to EMS. CRYZ encodes zeta-crystallin, a nicotinamide adenine dinucleotide phosphate-dependent quinone reductase that promotes insulin resistance via ubiquitination.37 Insulin resistance has been associated with chronic pelvic pain similar to that experienced by EMS patients, potentially linked to inflammatory responses and endocrine dysregulation.38 In ovarian cancer cells, CRYZ acts as a post-transcriptional regulator of Bcl-2 mRNA, inhibiting apoptosis and promoting chemoresistance; pharmacological or genetic inhibition of CRYZ restores chemosensitivity,39 suggesting that CRYZ may play a critical regulatory role in gynecological diseases. FRMD4B interacts with GRP1 and participates in insulin receptor and insulin-like growth factor-mediated PIP3-dependent signaling pathways.40 In EMS, activation of the PI3K/AKT signaling pathway significantly enhances pyroptosis and inflammatory factor levels, thereby contributing to disease pathogenesis.41 FZD6 is closely associated with the efficacy of immunotherapy. Moreover, it can promote melanoma invasion and metastasis by modulating the Wnt signaling and epithelial–mesenchymal transition pathways.42 The canonical Wnt signaling pathway is critically involved in embryonic development, cell proliferation, epithelial–mesenchymal transition, and carcinogenesis.43,44 Crosstalk between Wnt and TLR4/NF-κB signaling pathways has been implicated in chronic inflammation, disease progression, and tumorigenesis.45,46 In endometrial cells, it is shown that inhibiting the interaction between Wnt7a and FZD6 inactivates the Wnt/β-catenin signaling pathway,47 suggesting that FZD6 plays an important role in endometrial cell proliferation and survival. Knockout of TPST1 or TPST2 in animals severely impairs growth, reproductive function, and immune responses.48 Specifically, TPST1 knockout mice exhibit reduced litter sizes,49 indicating an important role in reproductive function. TPST1 may modulate immune and inflammatory responses by catalyzing sulfation,50,51 and impaired adhesion and migration of peritoneal macrophages may contribute to the development of EMS.52 Although no direct in vitro or in vivo functional validation of these four genes in EMS currently exists, evidence from other disease models, animal experiments, and endometrial tissue studies provides important clues for their potential roles. We next compared the diagnostic performance of our four-gene signature with CA-125, a widely studied serum biomarker for EMS. Previous studies have reported AUC values for CA-125 ranging from 0.79453 to 0.938.54 In our study, the four-gene signature demonstrated good diagnostic performance, with the nomogram achieving an AUC of 0.877 in the training set and 0.940 in external validation. These values are within the range reported for CA-125. However, these four genes are derived from the SUMOylation pathway, a novel post-translational modification mechanism that has not been previously explored for EMS diagnosis. Thus, our gene signature offers complementary biological information beyond traditional biomarkers like CA-125, providing new perspectives for early diagnosis and mechanistic insights into the disease. Future studies integrating clinical parameters with this gene signature may further improve diagnostic performance. GSEA of the four biomarkers revealed three commonly enriched biological pathways: neuroactive ligand–receptor interaction, cell cycle, and ubiquitin-mediated proteolysis. The neuroactive ligand–receptor pathway was known to mediate pain perception and was linked to EMS-associated pain via neuropeptide Y and VEGF.55,56 Studies have shown that upregulation of Cyclin D1 and CDK4 accelerates the G1 to S phase transition,57 while decreased expression of p21 and p27 leads to dysregulated cell cycle progression, promoting the formation and progression of endometriotic lesions.58 The ubiquitin–proteasome system influences the growth and maintenance of ectopic endometrial tissue.59 Reduced expression of the E3 ubiquitin ligase TRIM33 upregulates proteins associated with cell proliferation, such as TGFBR1 and α-SMA, thereby promoting cellular proliferation and fibrosis.60 Additionally, the E3 ubiquitin ligase MDM2 (murine double minute 2) regulates the stability of estrogen receptors, affecting the balance between cell proliferation and apoptosis. Overexpression of estrogen is closely associated with the development of EMS.61 These three biological pathways may provide theoretical support for understanding the pathogenesis of EMS. Immune infiltration analysis showed that Tregs, memory B cells, and monocytes were negatively correlated with the four biomarkers, suggesting their involvement in EMS pathogenesis. Conversely, resting mast cells were positively correlated with the biomarkers, indicating a potential protective role. At the single-cell level, all four biomarkers were differentially expressed in epithelial cells, pointing to active epithelial differentiation during disease progression. Pseudotime trajectory analysis revealed that epithelial cell differentiation could be divided into five stages, with significant changes in biomarker expression over time: CRYZ expression initially increased, then decreased, and later rose again; FRMD4B expression increased initially and then declined; FZD6 expression decreased in later stages; and TPST1 expression showed a downward trend. These patterns suggest that all four genes are involved throughout epithelial cell differentiation. Cell communication analysis indicated stronger and more frequent interactions between epithelial and endothelial cells. Epithelial cells may secrete cytokines that act on endothelial cells to regulate angiogenesis or inflammatory responses, while endothelial cells may release signaling molecules such as nitric oxide that influence epithelial cell metabolism or differentiation, forming a bidirectional regulatory network. This interplay may play a significant role in the development and progression of EMS. However, the specific molecular mechanisms and causal relationships involved still require experimental validation. Using the Miranda and MicroCosm databases, we predicted miRNAs targeting the biomarkers and constructed a lncRNA–miRNA–biomarker regulatory network. Specifically, hsa-miR-23a-3p, hsa-miR-656-3p, and hsa-miR-23b-3p were found to regulate TPST1 transcription. Both hsa-miR-23a-3p and hsa-miR-23b-3p are small non-coding RNAs involved in cell cycle regulation. miR-23b-3p can downregulate p21/CDKN1A expression, leading to G1 phase arrest and thus inhibiting cell proliferation.62,63 hsa-miR-656-3p may influence the synthesis and function of cell adhesion molecules by modulating the Wnt signaling pathway.64,65 These findings provide theoretical support for the role of the ceRNA network (lncRNA–miRNA–mRNA) in the pathogenesis of EMS. Prediction of TFs associated with the biomarkers using the NetworkAnalyst database revealed that ZFP2 simultaneously regulated the expression of CRYZ and TPST1. ZFP2 plays important roles in cell proliferation, differentiation, and response to oxidative stress.66 CRYZ could bind with ZFP2 and may influence gene expression by modulating the transcriptional activity of ZFP2.67 This interaction may be affected by changes in oxidative stress levels.68 ZFP2 may also regulate TPST1 expression by binding to its promoter region, thereby influencing transcriptional activity.69 ZFP2 expression is regulated by multiple signaling pathways, such as Wnt/β-catenin, which further affects TPST1 expression and function,70 implicating it in cell growth and development.71 ZFP2 may be a key factor in the pathogenesis of EMS, potentially influencing the disease through regulation of CRYZ and TPST1 expression, though further experimental validation is required. Several limitations of this study should be acknowledged. First, all conclusions are based on transcriptomic data from public databases and lack systematic wet-lab validation. The four SUMOylation-related biomarkers (CRYZ, FRMD4B, FZD6, and TPST1) were identified only at the transcriptional level, which cannot fully capture protein SUMOylation status given that SUMOylation is a post-translational modification. Furthermore, all samples in the datasets were derived from eutopic endometrium, precluding assessment of gene expression changes in ectopic lesions; thus, the tissue-specific diagnostic value of these biomarkers remains unclear. Second, the public datasets lack detailed clinical information, including age, hormonal status, menstrual cycle phase, and disease subtypes, making it impossible to adjust for potential confounding factors. Third, the diagnostic performance of the four biomarkers for early-stage EMS is limited, and the effect sizes observed in the MR analysis are small, suggesting that their clinical predictive utility may be constrained. Fourth, the GWAS data used for MR analysis were derived exclusively from European populations. Given that the prevalence of EMS varies globally, the generalizability of our conclusions to non-European populations requires further validation. To address these limitations, future studies should: (1) perform gene overexpression/knockdown and co-immunoprecipitation experiments in endometrial cell models to validate gene function and SUMOylation status; conduct proteasome activity assays and qPCR for cell-cycle regulators (Cyclin D1, CDK4, p21, p27) to validate GSEA findings; and elucidate the molecular pathways involved; (2) collect large-scale, well-phenotyped cohorts with paired eutopic and ectopic endometrial tissues to evaluate the diagnostic value and net benefit of the biomarkers for early-stage and different disease subtypes; and (3) conduct cross-population validation as multi-ancestry GWAS data become available.

In this study, we identified four SUMOylation-related genes (CRYZ, FRMD4B, FZD6, and TPST1) that are downregulated in EMS based on bioinformatics analyses. These genes may participate in the pathogenesis of EMS through pathways such as the cell cycle, neuroactive ligand-receptor interaction, and ubiquitin-mediated proteolysis. From a clinical perspective, these four SUMOylation-related genes demonstrate good diagnostic accuracy in distinguishing EMS patients from controls, suggesting their potential as auxiliary diagnostic biomarkers, particularly for the development of non-invasive diagnostic strategies. However, this study is entirely based on public transcriptomic data and computational analyses, lacking in vitro, in vivo, and prospective clinical cohort validation. Therefore, the above conclusions remain exploratory. Future studies combining functional experiments (eg, gene overexpression/knockdown) and independent large-scale clinical studies are needed to further validate the specific roles and diagnostic value of these genes in EMS. Abbreviations DEGs, differentially expressed genes; EMS, endometriosis; GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; GWAS, genome-wide association studies; IVs, instrumental variables; lncRNA, long non-coding RNA; miRNA, microRNA; MR, Mendelian randomization; MVMR, multivariable MR; OR, odds ratio; ROC, receiver operating characteristic; scRNA-seq, single-cell RNA sequencing; SNPs, single-nucleotide polymorphisms; SRGs, SUMOylation-related genes; SUMO, small ubiquitin-related modifier; TF, transcription factor; UVMR, univariable MR; WGCNA, weighted gene co-expression network analysis. Data Sharing Statement All data generated or analysed during this study are included in this published article and its Supplementary Information Files. Ethics Approval and Informed Consent All data used in this study were obtained from publicly available databases that have been anonymized, and the study involved only secondary analysis. In accordance with item 1 and 2 of Article 32 of the Measures for the Ethical Review of Life Science and Medical Research Involving Human Subjects (issued on February 18, 2023), this study is exempt from ethical review. Author Contributions All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work. Funding This study was supported by the 2021 Hainan Province Basic and Applied Basic Research Program (Hainan Province Natural Fund Youth Fund) (No. 821QN423), the Hainan Provincial Natural Science Foundation of China (No. 824RC561), and the Joint Program on Health Science & Technology Innovation of Hainan Province (No. WSJK2024QN067 and No. WSJK2026QN047). Disclosure The authors declare no conflicts of interest in this work.

1. Taylor HS, Kotlyar AM, Flores VA. Endometriosis is a chronic systemic disease: clinical challenges and novel innovations. Lancet. 2021;397(10276):839–22. doi:10.1016/S0140-6736(21)00389-5 2. Wei Y, Liang Y, Lin H, Dai Y, Yao S. Autonomic nervous system and inflammation interaction in endometriosis-associated pain. J Neuroinflammation. 2020;17(1):80. doi:10.1186/s12974-020-01752-1 3. Rahman L, Anwar R, Zulvayanti Z, Tjandraprawira KD. Rupture endometriomas presenting as acute abdomen infection in hasty and limited resources setting: a Pitfall not to miss - A case report. Int Med Case Rep J. 2024;17:635–641. doi:10.2147/IMCRJ.S472024 4. Gruber TM, Mechsner S. Pathogenesis of endometriosis: the origin of pain and subfertility. Cells. 2021;10(6):1381. doi:10.3390/cells10061381 5. Ferrero S, Esposito F, Abbamonte LH, Anserini P, Remorgida V, Ragni N. Quality of sex life in women with endometriosis and deep dyspareunia. Fertil Sterility. 2005;83(3):573–579. doi:10.1016/j.fertnstert.2004.07.973 6. Forouhar F, Mesbah N, Esmailpour S, Bastani P, Salimi M. Endometriosis presenting as a rare cause of intestinal perforation: a case report with literature review. Clin Case Rep. 2025;13(2):e70226. doi:10.1002/ccr3.70226 7. Somigliana E, Vigano P, Barbara G, Vercellini P. Treatment of endometriosis-related pain: options and outcomes. Front Biosci. 2009;1(2):455–465. doi:10.2741/e41 8. Zela-Coila F, Quispe-Vicuna C, Nunez-Lupaca JN, Aparicio-Curazi M, Goicochea-Lugo S. Comparison of clinical practice guidelines methods to reach diagnostic test recommendations regarding diagnostic laparoscopy for endometriosis: a scoping review. PLoS One. 2024;19(12):e0310593. doi:10.1371/journal.pone.0310593 9. Shafrir AL, Farland LV, Shah DK, et al. Risk for and consequences of endometriosis: a critical epidemiologic review. Best Pract Res Clin Obstet Gynaecol. 2018;51:1–15. doi:10.1016/j.bpobgyn.2018.06.001 10. Vercellini P, Vigano P, Somigliana E, Fedele L. Endometriosis: pathogenesis and treatment. Nat Rev Endocrinol. 2014;10(5):261–275. doi:10.1038/nrendo.2013.255 11. Gareau JR, Lima CD. The SUMO pathway: emerging mechanisms that shape specificity, conjugation and recognition. Nat Rev Mol Cell Biol. 2010;11(12):861–871. doi:10.1038/nrm3011 12. Zhu J, Wu P, Zeng C, Xue Q. Increased SUMOylation of TCF21 improves its stability and function in human endometriotic stromal cellsdagger. Biol Reprod. 2021;105(1):128–136. doi:10.1093/biolre/ioab038 13. Li SC, Yan LJ, Wei XL, Jia ZK, Yang JJ, Ning XH. A novel risk model of three SUMOylation genes based on RNA expression for potential prognosis and treatment sensitivity prediction in kidney cancer. Front Pharmacol. 2023;14:1038457. doi:10.3389/fphar.2023.1038457 14. Xia QD, Sun JX, Xun Y, et al. SUMOylation pattern predicts prognosis and indicates tumor microenvironment infiltration characterization in bladder cancer. Front Immunol. 2022;13:864156. doi:10.3389/fimmu.2022.864156 15. Yu L, Shen N, Shi Y, et al. Characterization of cancer-related fibroblasts (CAF) in hepatocellular carcinoma and construction of CAF-based risk signature based on single-cell RNA-seq and bulk RNA-seq data. Front Immunol. 2022;13:1009789. doi:10.3389/fimmu.2022.1009789 16. Mareckova M, Garcia-Alonso L, Moullet M, et al. An integrated single-cell reference atlas of the human endometrium. Nat Genet. 2024;56(9):1925–1937. doi:10.1038/s41588-024-01873-w 17. Liu Y, Gao G, Tian W, Lv Q, Liu D, Li C. Uncovering potential biomarkers of endometriosis: transcriptomic and single-cell analysis. Front Med Lausanne. 2025;12:1528434. doi:10.3389/fmed.2025.1528434 18. Shao W, Ju H, Xiahou Z, et al. Fibroblast heterogeneity and FN1-mediated signaling in endometriosis revealed by single-cell and spatial transcriptomics. Front Immunol. 2025;16:1680849. doi:10.3389/fimmu.2025.1680849 19. Liu S, Xie X, Lei H, Zou B, Xie L. Identification of key circRNAs/lncRNAs/miRNAs/mRNAs and pathways in preeclampsia using bioinformatics analysis. Med Sci Monit. 2019;25:1679–1693. doi:10.12659/MSM.912801 20. Xin S, Liu X, Li Z, et al. ScRNA-seq revealed an immunosuppression state and tumor microenvironment heterogeneity related to lymph node metastasis in prostate cancer. Exp Hematol Oncol. 2023;12(1):49. doi:10.1186/s40164-023-00407-0 21. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf. 2008;9:559. doi:10.1186/1471-2105-9-559 22. Zhang W, Shi G, Wang H, et al. Molecular mechanism of Xingnao Kaiqiao pill for perioperative neurocognitive disorder and its correlation with immune and inflammatory signaling pathways based on network pharmacology and molecular docking. Front Aging Neurosci. 2022;14:925072. doi:10.3389/fnagi.2022.925072 23. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation. 2021;2(3):100141. doi:10.1016/j.xinn.2021.100141 24. Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8(Suppl 4):S11. doi:10.1186/1752-0509-8-S4-S11 25. Lu F, Wu B, Wang Y. Mendelian randomization indicates that atopic dermatitis contributes to the occurrence of diabetes. BMC Med Genomics. 2023;16(1):132. doi:10.1186/s12920-023-01575-y 26. Bowden J, Davey Smith G, Burgess S. Mendelian randomization with invalid instruments: effect estimation and bias detection through Egger regression. Int J Epidemiol. 2015;44(2):512–525. doi:10.1093/ije/dyv080 27. Bowden J, Davey Smith G, Haycock PC, Burgess S. Consistent estimation in Mendelian randomization with some invalid instruments using a weighted median estimator. Genet Epidemiol. 2016;40(4):304–314. doi:10.1002/gepi.21965 28. Hemani G, Zheng J, Elsworth B, et al. The MR-Base platform supports systematic causal inference across the human phenome. Elife. 2018;7:e34408. doi:10.7554/eLife.34408 29. Burgess S, Scott RA, Timpson NJ, Davey Smith G, Thompson SG, Consortium E-I. Using published data in Mendelian randomization: a blueprint for efficient identification of causal risk factors. EurJ Epidemiol. 2015;30(7):543–552. doi:10.1007/s10654-015-0011-z 30. Hartwig FP, Davey Smith G, Bowden J. Robust inference in summary data Mendelian randomization via the zero modal pleiotropy assumption. Int J Epidemiol. 2017;46(6):1985–1998. doi:10.1093/ije/dyx102 31. Sui Z, Wu X, Du L, et al. Characterization of the immune cell infiltration landscape in esophageal squamous cell carcinoma. Front Oncol. 2022;12:879326. doi:10.3389/fonc.2022.879326 32. Liu TT, Li R, Huo C, et al. Identification of CDK2-related immune forecast model and ceRNA in lung adenocarcinoma, a pan-cancer analysis. Front Cell Dev Biol. 2021;9:682002. doi:10.3389/fcell.2021.682002 33. Yu Z, Qiu B, Zhou H, Li L, Niu T. Characterization and application of a lactate and branched chain amino acid metabolism related gene signature in a prognosis risk model for multiple myeloma. Cancer Cell Int. 2023;23(1):169. doi:10.1186/s12935-023-03007-4 34. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–457. doi:10.1038/nmeth.3337 35. Wu R, Zhang X, Zhang X, Sun L, Xia T, Zhang LJ. Deciphering the age-dependent changes of pulmonary fibroblasts in mice by single-cell transcriptomics. Front Cell Dev Biol. 2023;11:1287133. doi:10.3389/fcell.2023.1287133 36. Fang Z, Li J, Cao F, Li F. Integration of scRNA-Seq and Bulk RNA-Seq reveals molecular characterization of the immune microenvironment in acute pancreatitis. Biomolecules. 2022;13(1):78. doi:10.3390/biom13010078 37. Wei L, Tian Y, Chen Y, et al. Identification of TYW3/CRYZ and FGD4 as susceptibility genes for amyotrophic lateral sclerosis. Neurol Genet. 2019;5(6):e375. doi:10.1212/NXG.0000000000000375 38. Akbaba E, Sezgin B, Edgunlu T. The role of adropin, salusin-alpha, netrin-1, and nesfatin-1 in endometriosis and their association with insulin resistance. Turk J Obstet Gynecol. 2021;18(3):175–180. doi:10.4274/tjod.galenos.2021.12080 39. Lulli M, Trabocchi A, Roviello G, et al. Targeting z-Crystallin by aspirin restores the sensitivity to cisplatin in resistant A2780 ovarian cancer cells. Front Pharmacol. 2024;15:1377028. doi:10.3389/fphar.2024.1377028 40. Klarlund JK, Holik J, Chawla A, Park JG, Buxton J, Czech MP. Signaling complexes of the FERM domain-containing protein GRSP1 bound to ARF exchange factor GRP1. J Biol Chem. 2001;276(43):40065–40070. doi:10.1074/jbc.M105260200 41. An M, Fu X, Meng X, et al. PI3K/AKT signaling pathway associates with pyroptosis and inflammation in patients with endometriosis. J Reprod Immunol. 2024;162:104213. doi:10.1016/j.jri.2024.104213 42. Dong B, Simonson L, Vold S, et al. FZD6 promotes melanoma cell invasion but not proliferation by regulating canonical Wnt signaling and epithelial‒mesenchymal transition. J Invest Dermatol. 2023;143(4):621–629e626. doi:10.1016/j.jid.2022.09.658 43. Yang YZ. Wnt signaling in development and disease. Cell Biosci. 2012;2(1):14. doi:10.1186/2045-3701-2-14 44. Sompel K, Elango A, Smith AJ, Tennis MA. Cancer chemoprevention through frizzled receptors and EMT. Discov Oncol. 2021;12(1):32. doi:10.1007/s12672-021-00429-2 45. Du Q, Geller DA. Cross-regulation between Wnt and NF-κB signaling pathways. Forum Immunopathol Dis Therap. 2010;1(3):155–181. doi:10.1615/ForumImmunDisTher.v1.i3.10 46. Ma B, Hottiger MO. Crosstalk between Wnt/beta-Catenin and NF-kappaB signaling pathway during inflammation. Front Immunol. 2016;7:378. doi:10.3389/fimmu.2016.00378 47. Chandra V, Fatima I, Manohar M, et al. Inhibitory effect of 2-(piperidinoethoxyphenyl)-3-(4-hydroxyphenyl)-2H-benzo(b)pyran (K-1) on human primary endometrial hyperplasial cells mediated via combined suppression of Wnt/beta-catenin signaling and PI3K/Akt survival pathway. Cell Death Dis. 2014;5(8):e1380. doi:10.1038/cddis.2014.334 48. Westmuckett AD, Hoffhines AJ, Borghei A, Moore KL. Early postnatal pulmonary failure and primary hypothyroidism in mice with combined TPST-1 and TPST-2 deficiency. Gen Comp Endocrinol. 2008;156(1):145–153. doi:10.1016/j.ygcen.2007.12.006 49. Ouyang YB, Crawley JT, Aston CE, Moore KL. Reduced body weight and increased postimplantation fetal death in tyrosylprotein sulfotransferase-1-deficient mice. J Biol Chem. 2002;277(26):23781–23787. doi:10.1074/jbc.M202420200 50. Wang E, Wang Y, Zhou S, et al. Identification of three hub genes related to the prognosis of idiopathic pulmonary fibrosis using bioinformatics analysis. Int J Med Sci. 2022;19(9):1417–1429. doi:10.7150/ijms.73305 51. Nakamura N, Shimaoka Y, Tougan T, et al. Isolation and expression profiling of genes upregulated in bone marrow-derived mononuclear cells of rheumatoid arthritis patients. DNA Res. 2006;13(4):169–183. doi:10.1093/dnares/dsl006 52. Soni UK, Tripathi R, Jha RK. MCP-1 exerts the inflammatory response via ILK activation during endometriosis pathogenesis. Life Sci. 2024;353:122902. doi:10.1016/j.lfs.2024.122902 53. Szubert M, Suzin J, Wierzbowski T, Kowalczyk-Amico K. CA-125 concentration in serum and peritoneal fluid in patients with endometriosis - preliminary results. Arch Med Sci. 2012;8(3):504–508. doi:10.5114/aoms.2012.29529 54. Tuten A, Kucur M, Imamoglu M, et al. Copeptin is associated with the severity of endometriosis. Arch Gynecol Obstet. 2014;290(1):75–82. doi:10.1007/s00404-014-3163-2 55. Varga J, Reviczka A, Hakova H, Svajdler P, Rabajdova M, Ostro A. Predictive factors of endometriosis progression into ovarian cancer. J Ovarian Res. 2022;15(1):5. doi:10.1186/s13048-021-00940-8 56. Gyorfi M, Rupp A, Abd-Elsayed A. Fibromyalgia pathophysiology. Biomedicines. 2022;10(12):3070. doi:10.3390/biomedicines10123070 57. Song J, Ham J, Park S, et al. Alpinumisoflavone activates disruption of calcium homeostasis, mitochondria and autophagosome to suppress development of endometriosis. Antioxidants. 2023;12(7):1324. doi:10.3390/antiox12071324 58. Poli-Neto OB, Meola J, Rosa ESJC, Tiezzi D. Transcriptome meta-analysis reveals differences of immune profile between eutopic endometrium from stage I-II and III-IV endometriosis independently of hormonal milieu. Sci Rep. 2020;10(1):313. doi:10.1038/s41598-019-57207-y 59. Cassidy K, Zhao H. Redefining the scope of targeted protein degradation: translational opportunities in hijacking the autophagy-lysosome pathway. Biochemistry. 2023;62(3):580–587. doi:10.1021/acs.biochem.1c00330 60. Yang M, Jiang H, Ding X, et al. Multi-omics integration highlights the role of ubiquitination in endometriosis fibrosis. J Transl Med. 2024;22(1):445. doi:10.1186/s12967-024-05245-0 61. Chen LJ, Hu B, Han ZQ, et al. BAG2-mediated inhibition of CHIP expression and overexpression of MDM2 contribute to the initiation of endometriosis by modulating estrogen receptor status. Front Cell Dev Biol. 2020;8:554190. doi:10.3389/fcell.2020.554190 62. Zhang YS, Wang MY, Zhang WL, Tang CH. Proliferation, migration and apoptosis of acute myeloid leukemia cells regulated by mir-23a-3p targeting SMC1A and the mechanism. Zhonghua Zhong Liu Za Zhi. 2019;41(10):753–759. doi:10.3760/cma.j.issn.0253-3766.2019.10.006 63. Lataster L, Huber HM, Bottcher C, Foller S, Takors R, Radziwill G. Cell cycle control by optogenetically regulated cell cycle inhibitor protein p21. Biology. 2023;12(9):1194. doi:10.3390/biology12091194 64. Xiao D, Xiong M, Wang X, et al. Regulation of the function and expression of EpCAM. Biomedicines. 2024;12(5):1129. doi:10.3390/biomedicines12051129 65. Carter JK, Tsai MC, Venturini N, et al. Stellate cell-specific adhesion molecule protocadherin 7 regulates sinusoidal contraction. Hepatology. 2024;80(3):566–577. doi:10.1097/HEP.0000000000000782 66. Alam G, Luan Z, Gul A, et al. Activation of farnesoid X receptor (FXR) induces crystallin zeta expression in mouse medullary collecting duct cells. Pflugers Arch. 2020;472(11):1631–1641. doi:10.1007/s00424-020-02456-4 67. Nasher F, Taylor AJ, Elmi A, et al. MdaB and NfrA, two novel reductases important in the survival and persistence of the major Enteropathogen Campylobacter jejuni. J Bacteriol. 2022;204(1):e0042121. doi:10.1128/JB.00421-21 68. Alves F, Lane D, Nguyen TPM, Bush AI, Ayton S. In defence of ferroptosis. Signal Transduct Target Ther. 2025;10(1):2. 69. Chen Y, Zhou D, Qian X, Ge S, Shuai Z. Characteristic changes in the mRNA expression profile of plasma exosomes from patients with MPO-ANCA-associated vasculitis and its possible correlations with pathogenesis. Clin Exp Med. 2024;24(1):222. doi:10.1007/s10238-024-01457-2 70. Xu T, Xu W, Zheng Y, et al. Comprehensive FGFR3 alteration-related transcriptomic characterization is involved in immune infiltration and correlated with prognosis and immunotherapy response of bladder cancer. Front Immunol. 2022;13:931906. doi:10.3389/fimmu.2022.931906 71. Luo M, Li X, Zhang J, Miao Y, Liu D. The C3H gene PtZFP2-like in Pinellia ternata acts as a positive regulator of the resistance to soft rot caused by Pectobacterium carotovorum. Physiol Plant. 2025;177(1):e70121. doi:10.1111/ppl.70121 © 2026 The Author(s). This work is published and licensed by Dove Medical Press Limited. The full terms of this license are available at https://www.dovepress.com/terms and incorporate the Creative Commons Attribution - Non Commercial (unported, 4.0) License. By accessing the work you hereby accept the Terms. Non-commercial uses of the work are permitted without any further permission from Dove Medical Press Limited, provided the work is properly attributed. For permission for commercial use of this work, please see paragraphs 4.2 and 5 of our Terms. Recommended articles Identification of Critical Modules and Biomarkers of Ulcerative Colitis by Using WGCNA Yuan Y, Li N, Fu M, Ye M Journal of Inflammation Research 2023, 16:1611-1628 Published Date: 17 April 2023 Bioinformatics Analysis and Experimental Validation of Mitochondrial Autophagy Genes in Knee Osteoarthritis Tang K, Sun L, Chen L, Feng X, Wu J, Guo H, Zheng Y International Journal of General Medicine 2024, 17:639-650 Published Date: 23 February 2024 Causal Relationship Between Endometriosis and Pelvic Inflammatory Diseases: Mendelian Randomization Study Liu K, Liu X, Cao T, Cui X, Sun P, Zhang L, Wu X International Journal of Women's Health 2024, 16:727-735 Published Date: 24 April 2024 The Effect of Circulating Inflammatory Proteins on Endometriosis: A Mendelian Randomization Study Wei Y, Zhao X, Li L ImmunoTargets and Therapy 2024, 13:585-593 Published Date: 1 November 2024 Endometriosis Severity and Risk of Preeclampsia: A Combined Mendelian Randomization and Observational Study Zu Y, Xie Y, Zhang H, Chen L, Yan S, Wang Z, Fang Z, Lin S, Yan J International Journal of Women's Health 2025, 17:923-935 Published Date: 27 March 2025