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 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.

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 ). Figure 1 Cellular composition and activated pathways in EMS based on scRNA-seq data. ( A ) Visualization of the top 2,000 highly variable genes (highlighted in red). ( B ) UMAP dimensionality reduction plot based on the top 50 principal components. ( C ) UMAP plot showing 30 initial cell clusters. ( D ) UMAP plot annotating the 9 distinct cell types. ( E ) Stacked bar chart displaying the relative proportions of the 9 cell types in EMS and control samples. ( F ) Box plot comparing the abundance of cells between EMs and controls. *P<0.05; ns, not significant. GSVA reveals biological pathways significantly activated in epithelial ( G ) cells and fibroblasts ( H ) of EMS samples, such as Notch signaling, angiogenesis, and epithelial-mesenchymal transition. The image A showing a scatter plot of Average Expression (x-axis) vs Standardized Variance (y-axis). Text: Non-variable count: 24386. Variable count: 2000. Labeled points include HBB, TPSB2, SCGB1D4, HBA1, HBA2, TPSAB1, PPBP, S100A8, S100A9, IGKC. Variable genes appear as a highlighted subset above the main point cloud. The image B showing a scatter plot of PC (x-axis) vs Standard Deviation (y-axis). PC ranges 0 to 50. Standard Deviation ranges about 1.0 to 10.0. Points decrease from about (1, 10.0) toward about (50, 1.0). The image C showing a UMAP scatter plot with umap1 (x-axis) and umap2 (y-axis). Cluster labels in the legend run from 0 to 29, indicating 30 clusters. The image D showing a UMAP scatter plot with umap1 (x-axis) and umap2 (y-axis) annotated by cell type. Legend labels: Endothelialcells, Epithelialcells, Fibroblasts, Monocytes, MSC, Neutrophils, NKcell, Smoothmusclecells, Tissuestemcells. Text labels appear on the map for several groups including Endothelialcells, Epithelialcells, Fibroblasts, Neutrophils, NKcell, MSC, Smoothmusclecells, Tissuestemcells. The image E showing a stacked bar chart of Proportion (y-axis, 0.00 to 1.00) vs Celltype (x-axis categories: Normal, Tumor, All). Legend labels: Neutrophils, Monocytes, Smoothmusclecells, Tissuestemcells, Endothelialcells, NKcell, MSC, Epithelialcells, Fibroblasts. The Tumor bar contains a larger Epithelialcells segment than the Normal bar, while the Normal bar contains a larger Fibroblasts segment than the Tumor bar. The image F showing box plots of Cell percentage (y-axis) by cell type (x-axis categories: Endothelialcells, Epithelialcells, Fibroblasts, Monocytes, MSC, Neutrophils, NKcell, Smoothmusclecells, Tissuestemcells) with group legend: EMS and Normal. Multiple comparisons are labeled ns and one comparison is marked with an asterisk. The Epithelialcells box for EMS is higher than Normal and the Fibroblasts box for EMS is lower than Normal. The image G showing a clustered heatmap with columns labeled EMS and Normal and a scale bar from negative 0.4 to positive 0.4. Row labels include HALLMARKMY CTARGETSV2, HALLMARKOXIDATIVEPHOSPHORYLATION, HALLMARKESTROGENRESPONSEEARLY, HALLMARKESTROGENRESPONSELATE, HALLMARKTGFBETASIGNALING, HALLMARKNOTCHSIGNALING, HALLMARKANGIOGENESIS, HALLMARKEPITHELIALMESENCHYMALTRANSITION and additional HALLMARK pathways. Many rows show higher values in EMS than Normal for pathways including HALLMARKNOTCHSIGNALING, HALLMARKANGIOGENESIS and HALLMARKEPITHELIALMESENCHYMALTRANSITION. The image H showing a second clustered heatmap with columns labeled EMS and Normal and a scale bar from negative 0.4 to positive 0.4. Row labels include HALLMARKANGIOGENESIS, HALLMARKEPITHELIALMESENCHYMALTRANSITION, HALLMARKTNFALPHASIGNALINGVIANFKB, HALLMARKNOTCHSIGNALING, HALLMARKAPICALSURFACE, HALLMARKCHOLESTEROLHOMEOSTASIS, HALLMARKGLYCOLYSIS, HALLMARKOXIDATIVEPHOSPHORYLATION, HALLMARKHYPOXIA, HALLMARKTGFBETASIGNALING and additional HALLMARK pathways. Multiple pathways show higher values in EMS than Normal, including HALLMARKANGIOGENESIS and HALLMARKEPITHELIALMESENCHYMALTRANSITION. Across images C to F, the 30 UMAP clusters in image C correspond to the 9 annotated cell types in image D and the same 9 cell types are compared by relative proportion in image E and by cell percentage distributions between EMS and Normal in image F. A multi-graph figure showing scRNA-seq gene variability, UMAP cell types, cell proportions and pathways in EMS. Cellular composition and activated pathways in EMS based on scRNA-seq data. ( A ) Visualization of the top 2,000 highly variable genes (highlighted in red). ( B ) UMAP dimensionality reduction plot based on the top 50 principal components. ( C ) UMAP plot showing 30 initial cell clusters. ( D ) UMAP plot annotating the 9 distinct cell types. ( E ) Stacked bar chart displaying the relative proportions of the 9 cell types in EMS and control samples. ( F ) Box plot comparing the abundance of cells between EMs and controls. *P<0.05; ns, not significant. GSVA reveals biological pathways significantly activated in epithelial ( G ) cells and fibroblasts ( H ) of EMS samples, such as Notch signaling, angiogenesis, and epithelial-mesenchymal transition. 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 ). Figure 2 Identification and functional analysis of 49 candidate genes in EMS. ( A ) Volcano plot of DEGs in the scRNA-seq dataset ( GSE214411 ). ( B ) Volcano plot of DEGs in the bulk transcriptomic dataset ( GSE51981 ). ( C ) Sample clustering and trait heatmap indicating no outlier samples in the GSE51981 dataset. ( D ) Analysis of network topology for various soft-thresholding powers (β). The red line indicates the selected β value of 16. ( E ) Gene clustering dendrogram and module assignment from WGCNA. ( F ) Module–trait relationships heatmap. The MEturquoise module shows the highest positive correlation with the SUMOylation score. ( G ) Scatterplot of Gene Significance for the SUMOylation score versus Module Membership in the MEturquoise module. GO ( H ) and KEGG ( I ) enrichment analysis of the 49 candidate genes. ( J ) Protein-protein interaction network of the candidate genes, highlighting key interactions. A multi-panel scientific infographic arranged in a grid and labeled A to J, summarizing DEG and WGCNA results and downstream analyses for 49 candidate genes. The image A showing a volcano plot for scRNA-seq dataset GSE214411 with x-axis label log2FC and y-axis label minus log10(adj.P.Val) and side labels DOWN, NOT, UP. The image B showing a volcano plot for bulk transcriptomic dataset GSE51981 with x-axis label log2(Fold Change) and y-axis label minus log10(adj.P.Val) and side labels DOWN, NOT, UP. The image C showing sample clustering dendrogram with y-axis label Height and a trait heatmap strip labeled score. The image D showing two plots: left plot with x-axis label Soft Threshold (power) and y-axis label Scale Free Topology Model Fit, signed R superscript 2; right plot with x-axis label Soft Threshold (power) and y-axis label Mean Connectivity. The image E showing a gene clustering dendrogram with y-axis label Height and a strip labeled Module colors. The image F showing a module–trait relationships heatmap with row labels MEbrown, MEturquoise, MEblue, MEyellow, column label score and a scale from minus 1 to 1. Cell text includes 0.77 (5.2e-27), 0.8 (5.3e-26), minus 0.69 (9e-17) and minus 0.42 (4.6e-09). The image G showing a scatterplot with x-axis label Module Membership in turquoise module and y-axis label Gene Significance for score. The image H showing a GO enrichment circular plot with labels BP biological process, CC cellular component and MF molecular function. The image I showing a KEGG enrichment dot plot with y-axis label Pathway and x-axis tick marks 1.25, 1.50, 1.75, 2.00 and pathway labels: Cushing syndrome, Hippo signaling pathway, Protein processing in endoplasmic reticulum, Cytokine-cytokine receptor interaction, Pathways of neurodegeneration - multiple diseases, SNARE interactions in vesicular transport, Cholesterol metabolism, Glutathione metabolism, Basal cell carcinoma, Cortisol synthesis and secretion. The image J showing a protein-protein interaction network diagram with labeled gene nodes connected by lines. A multi-panel infographic of DEG and WGCNA results linking MEturquoise to SUMOylation score. Identification and functional analysis of 49 candidate genes in EMS. ( A ) Volcano plot of DEGs in the scRNA-seq dataset ( GSE214411 ). ( B ) Volcano plot of DEGs in the bulk transcriptomic dataset ( GSE51981 ). ( C ) Sample clustering and trait heatmap indicating no outlier samples in the GSE51981 dataset. ( D ) Analysis of network topology for various soft-thresholding powers (β). The red line indicates the selected β value of 16. ( E ) Gene clustering dendrogram and module assignment from WGCNA. ( F ) Module–trait relationships heatmap. The MEturquoise module shows the highest positive correlation with the SUMOylation score. ( G ) Scatterplot of Gene Significance for the SUMOylation score versus Module Membership in the MEturquoise module. GO ( H ) and KEGG ( I ) enrichment analysis of the 49 candidate genes. ( J ) Protein-protein interaction network of the candidate genes, highlighting key interactions. 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 Symbol Exposure Outcome Method nSNP Beta-value SE P value FRMD4B eqtl-a-ENSG00000114541 Endometriosis Inverse variance weighted 13 −0.000169123 8.36E-05 0.042973949 CRYZ eqtl-a-ENSG00000116791 Endometriosis Inverse variance weighted 13 −0.000129299 5.08E-05 0.010912239 PLK2 eqtl-a-ENSG00000145632 Endometriosis Inverse variance weighted 9 0.001665824 0.000415086 5.99E-05 FZD6 eqtl-a-ENSG00000164930 Endometriosis Inverse variance weighted 8 −0.000417723 0.000150319 0.005454279 ZNF606 eqtl-a-ENSG00000166704 Endometriosis Inverse variance weighted 7 −0.000705308 0.000298719 0.018220561 TPST1 eqtl-a-ENSG00000169902 Endometriosis Inverse variance weighted 78 −0.000201089 2.60E-05 1.05E-14 PPA1 eqtl-a-ENSG00000180817 Endometriosis Inverse variance weighted 17 0.000207698 6.02E-05 0.000563279 Table 2 Results of the Heterogeneity Test Exposure Outcome Method Q Q_df Q_P value eqtl-a-ENSG00000164930 Endometriosis MR Egger 4.404766397 6 0.62207478 eqtl-a-ENSG00000164930 Endometriosis Inverse variance weighted 4.579476873 7 0.711126991 Table 3 Results of the Horizontal Pleiotropy Test Exposure Outcome Egger_Intercept SE P value eqtl-a-ENSG00000164930 Endometriosis 3.25E-05 7.77E-05 0.690505481 Table 4 Results of Horizontal Pleiotropy Test (Presso) Exposure Outcome RSSobs P value eqtl-a-ENSG00000164930 Endometriosis 5.37804293472148 0.819 Figure 3 Mendelian Randomization analysis identifies key causal genes for EMs. ( A ) Scatter plots of the univariable MR analysis for seven candidate genes, showing the causal effect of gene expression on EMS risk using five MR methods. ( B ) Forest plot displaying the MR effect size and 95% confidence interval for each candidate gene. ( C ) Funnel plots assessing the symmetry of MR estimates. ( D ) Leave-one-out sensitivity analysis plots. ( E ) Forest plot from the MVMR analysis, showing the odds ratio and 95% confidence interval for each gene after adjusting for other exposures. The image A showing seven scatter plots titled MR Scatter. Each plot has x axis label SNP effect on exposure and y axis label SNP effect on outcome. Each plot contains points with vertical and horizontal error bars and five fitted lines labeled MR Egger, Weighted median, Inverse variance weighted, Simple mode and Weighted mode. The fitted line slopes vary by gene, with several negative slopes and some positive slopes. The image B showing seven forest plots. Each plot has x axis label MR effect size and y axis label SNP. Each plot shows multiple horizontal confidence interval lines with a central point estimate per SNP and a vertical reference line at 0. The image C showing seven funnel plots titled MR Method and Inverse variance weighted. Each plot has x axis label MR effect size and y axis label Standard error. Each plot shows scattered points around a vertical reference line. The image D showing seven leave one out plots. Each plot has x axis label MR effect size and y axis label SNP. Each plot shows horizontal confidence intervals for leave one out estimates and a vertical reference line. The image E showing a forest plot table with columns Exposure, P value and Odd Ratio (95 percent CI). The x axis label is Odd Ratio and the y axis lists exposures FRMD4B, CRYZ, PLK2, FZD6, ZNF606, TPST1 and PPA1. The x axis range is 0.98 to 1.1 with ticks at 0.98, 0.99, 1 and 1.1. Values shown: FRMD4B P value less than 0.001, odd ratio 0.9813 (0.97048 to 0.99226); CRYZ P value 0.02728, odd ratio 1.0155 (1.00162 to 1.02772); PLK2 P value 0.76265, odd ratio 0.99571 (0.96829 to 1.0239); FZD6 P value 0.6884, odd ratio 0.99805 (0.98856 to 1.00763); ZNF606 P value 0.00275, odd ratio 0.99812 (0.99689 to 0.99935); TPST1 P value 0.000303, odd ratio 0.9987 (0.99812 to 0.99962); PPA1 P value 0.00182, odd ratio 0.99436 (0.99084 to 0.9979). A mixed figure showing seven scatter plots and multiple forest, funnel and leave one out plots for genes. Univariate MR Analysis of Exposure Factors Results of the Heterogeneity Test Results of the Horizontal Pleiotropy Test Results of Horizontal Pleiotropy Test (Presso) Mendelian Randomization analysis identifies key causal genes for EMs. ( A ) Scatter plots of the univariable MR analysis for seven candidate genes, showing the causal effect of gene expression on EMS risk using five MR methods. ( B ) Forest plot displaying the MR effect size and 95% confidence interval for each candidate gene. ( C ) Funnel plots assessing the symmetry of MR estimates. ( D ) Leave-one-out sensitivity analysis plots. ( E ) Forest plot from the MVMR analysis, showing the odds ratio and 95% confidence interval for each gene after adjusting for other exposures. 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 Symbol Exposure Outcome P value OR 95% CI (Lower Limit) 95% CI (Upper limit) FRMD4B eqtl-a-ENSG00000114541 Endometriosis 0.00086 0.981307914 0.970480682 0.992255942 CRYZ eqtl-a-ENSG00000116791 Endometriosis 0.02728 1.014587661 1.001624579 1.027718512 PLK2 eqtl-a-ENSG00000145632 Endometriosis 0.76265 0.995706825 0.968288017 1.023902045 FZD6 eqtl-a-ENSG00000164930 Endometriosis 0.6884 0.998047464 0.988559248 1.007626747 ZNF606 eqtl-a-ENSG00000166704 Endometriosis 0.00275 0.99812145 0.996893585 0.999350828 TPST1 eqtl-a-ENSG00000169902 Endometriosis 0.00303 0.998869867 0.998123357 0.999616934 PPA1 eqtl-a-ENSG00000180817 Endometriosis 0.00182 0.994364048 0.990837817 0.997902828 Multivariate MR 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. Figure 4 Diagnostic value of candidate key genes and construction of a predictive nomogram. ( A ) ROC curves for the five candidate key genes (CRYZ, FRMD4B, FZD6, TPST1, ZNF606) in the GSE51981 dataset. Expression validation of the five identified biomarkers in the GSE51981 ( B ) and GSE7305 ( C ) datasets. *P<0.05; **P<0.01; ****P<0.0001; ns, not significant. ( D ) Nomogram for predicting the risk of EMs based on the expression levels of the four biomarkers (CRYZ, FRMD4B, FZD6, TPST1). ( E ) Calibration curve of the nomogram. The dashed line represents the ideal prediction, and the solid line represents the performance of the nomogram. ( F ) Decision curve analysis for the nomogram and individual biomarkers. ( G ) The ROC curve of the nomogram model for predicting the risk of endometriosis. The image contains multiple graphs analyzing diagnostic and model performance for EMs versus control. Image A shows ROC curves for genes CRYZ, FRMD4B, FZD6, TPST1, ZNF606 with AUCs: CRYZ 0.791, FRMD4B 0.798, FZD6 0.824, TPST1 0.824, ZNF606 0.869. Axes: 1-Specificity vs Sensitivity. Image B presents box plots of gene expression for CRYZ, FRMD4B, FZD6, TPST1, ZNF606 in control and EMs groups, with significance marked as . Expression is higher in controls. Image C includes box plots with significance levels:, *, ns, showing similar trends. Image D features a nomogram with scales for Points (0-100), CRYZ, FRMD4B, FZD6, TPST1, Total Points (0-220) and Risk (0.1, 0.5, 0.9), used to predict EMs risk. Image E shows a calibration curve with axes: Nomogram-Predicted Probability vs Actual disease proportion, including Apparent, Bias-corrected, Ideal. Hosmer-Lemeshow test: B=1000, boot Mean absolute error 0.038, n=111. Image F displays decision curve analysis with axes: High Risk Threshold vs Net Benefit, including CRYZ, FRMD4B, FZD6, TPST1, pathologic_nomogram, All, None. Image G presents a ROC curve with AUC 0.877, axes: 1-Specificity vs Sensitivity. The graphs collectively demonstrate gene expression differences, model prediction accuracy and clinical utility. Graphs: ROC curves, boxplots, nomogram, calibration, decision analysis, model ROC for EMs vs control. Diagnostic value of candidate key genes and construction of a predictive nomogram. ( A ) ROC curves for the five candidate key genes (CRYZ, FRMD4B, FZD6, TPST1, ZNF606) in the GSE51981 dataset. Expression validation of the five identified biomarkers in the GSE51981 ( B ) and GSE7305 ( C ) datasets. *P<0.05; **P<0.01; ****P<0.0001; ns, not significant. ( D ) Nomogram for predicting the risk of EMs based on the expression levels of the four biomarkers (CRYZ, FRMD4B, FZD6, TPST1). ( E ) Calibration curve of the nomogram. The dashed line represents the ideal prediction, and the solid line represents the performance of the nomogram. ( F ) Decision curve analysis for the nomogram and individual biomarkers. ( G ) The ROC curve of the nomogram model for predicting the risk of endometriosis. 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). 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. Figure 5 Functional enrichment of the biomarkers. GSEA results for CRYZ ( A ), FRMD4B ( B ), FZD6 ( C ), and TPST1 ( D ), showing commonly enriched pathways. ( E ) Gene-gene interaction network for the four biomarkers, constructed using GeneMANIA. A) A multi-line enrichment plot. The x-axis is labeled Rank in Ordered Dataset (no unit), ranging from 0 to about 20000. The y-axis is labeled Running Enrichment Score (no unit), ranging from about minus 0.50 to plus 0.50, with a dashed baseline at 0. Multiple curves rise to peaks near ranks about 3000 to 6000 (highest curve near about plus 0.50), then decline toward 0 by about rank 20000; two curves dip below 0 to about minus 0.30 around ranks about 9000 to 14000 before returning toward 0. A middle band shows many vertical tick marks across the rank range. A lower subplot is labeled Ranked List Metric (no unit), showing a filled shape decreasing from about plus 1.0 at rank 0 to about minus 1.0 at rank 20000. B) Same plot structure and axes as A. Curves peak near ranks about 3000 to 6000 at about plus 0.45 to plus 0.55, then trend down toward 0 by about rank 20000; one curve drops below 0 to about minus 0.30 near ranks about 9000 to 14000. Tick-mark band and the Ranked List Metric subplot again span about plus 1.0 to minus 1.0. C) Same plot structure. The x-axis is Rank in Ordered Dataset (no unit) from 0 to about 20000. The y-axis is Running Enrichment Score (no unit) from about minus 0.50 to plus 0.60. Several curves peak near about plus 0.55 around ranks about 4000 to 7000, then decline toward 0; one curve stays below 0 around minus 0.20 to minus 0.30 through mid ranks before approaching 0 near the end. Tick-mark band and Ranked List Metric subplot run from about plus 1.0 to minus 1.0. D) Same plot structure. The x-axis is Rank in Ordered Dataset (no unit) from 0 to about 20000. The y-axis is Running Enrichment Score (no unit) from about minus 0.50 to plus 0.60. Curves peak near ranks about 3000 to 6000 at about plus 0.55 to plus 0.60, then decrease toward 0; one curve remains negative near about minus 0.20 across much of the range. Tick-mark band and Ranked List Metric subplot again span about plus 1.0 to minus 1.0. E) A circular gene interaction network diagram with four central nodes labeled CRYZ, FRMD4B, FZD6 and TPST1 connected to many surrounding gene nodes by numerous lines. Two legends list Network edge types (Physical Interactions, Co-expression, Predicted, Co-localization, Genetic Interactions, Pathway, Shared protein domains) and Functions (non-canonical Wnt signaling pathway, ciliary membrane, photoreceptor outer segment, detection of visible light, detection of light stimulus, 9 plus 0 non-motile cilium, photoreceptor cell cilium). Five-panel line plots and a network diagram showing enriched pathways for CRYZ, FRMD4B, FZD6 and TPST1. Functional enrichment of the biomarkers. GSEA results for CRYZ ( A ), FRMD4B ( B ), FZD6 ( C ), and TPST1 ( D ), showing commonly enriched pathways. ( E ) Gene-gene interaction network for the four biomarkers, constructed using GeneMANIA. 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 ). 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. Figure 6 Immune infiltration landscape and its correlation with biomarkers in EMs. ( A ) Bar plot showing the relative proportions of 22 immune cell types in each sample from the GSE51981 dataset. ( B ) Correlation heatmap of the 22 immune cell types. Red indicates positive correlations, and blue indicates negative correlations. *P<0.05; **P<0.01; ***P<0.001. ( C ) Box plot comparing the infiltration levels of immune cell types between EMS and control samples. *P<0.05; **P<0.01; ns, not significant. ( D ) Correlation heatmap between the four biomarkers and differentially infiltrated immune cells. ( E ) Drug-biomarker interaction network, depicting potential therapeutic drugs targeting the four biomarkers. Image A shows a stacked bar chart of 22 immune cell types in GSE51981 samples, with relative fractions on the Y-axis and samples on the X-axis, highlighting the mix of immune cells. Image B displays a correlation heatmap for these cell types, with significance markers: *P<0.05, **P<0.01, ***P<0.001. Image C compares immune cell infiltration between EMS and control using box plots, showing some cells are more prevalent in EMS, others less, marked by *P<0.05, **P<0.01 and ns. Image D features plots for CRYZ, FRMD4B, FZD6 and TPST1, with correlation coefficients on the X-axis and immune cell types on the Y-axis, including Tregs, follicular helper T cells, monocytes, resting mast cells, activated dendritic cells and memory B cells. Points and lines indicate correlation direction and magnitude, with p-values <0.001, <0.01, 0.05. Image E depicts a drug-biomarker interaction network linking drugs to biomarkers CRYZ, FRMD4B, FZD6 and TPST1, suggesting multiple drug candidates. An infographic on immune infiltration, biomarker correlations and drug interactions in EMs. Immune infiltration landscape and its correlation with biomarkers in EMs. ( A ) Bar plot showing the relative proportions of 22 immune cell types in each sample from the GSE51981 dataset. ( B ) Correlation heatmap of the 22 immune cell types. Red indicates positive correlations, and blue indicates negative correlations. *P<0.05; **P<0.01; ***P<0.001. ( C ) Box plot comparing the infiltration levels of immune cell types between EMS and control samples. *P<0.05; **P<0.01; ns, not significant. ( D ) Correlation heatmap between the four biomarkers and differentially infiltrated immune cells. ( E ) Drug-biomarker interaction network, depicting potential therapeutic drugs targeting the four biomarkers. 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 ). Figure 7 Expression dynamics of biomarkers in epithelial cells and intercellular communication. ( A ) Box plots showing the expression levels of the four biomarkers across different cell types. Epithelial cells were identified as the key cell type. ***P<0.001; ns, not significant. ( B ) UMAP plot of epithelial cells from the scRNA-seq dataset. ( C ) Pseudotime trajectory of epithelial cells, colored by differentiation state. ( D ) Pseudotime trajectory of epithelial cells, colored by the inferred pseudotime. ( E ) Expression dynamics of the four biomarkers along the pseudotime trajectory of epithelial cell differentiation. ( F ) The number of inferred interactions between epithelial cells and other cell types. ( G ) The strength of inferred interactions between epithelial cells and other cell types. ( H ) Circle plot illustrating the comprehensive cell-cell communication network between major cell types, with a highlight on the significant interactions between epithelial and endothelial cells. Image A shows box plots of expression levels for CRYZ, FRMD4B, FZD6, TPST1 by group (EMS and Normal). TPST1 has the highest expression. Image B is a scatter plot of epithelial cells with a pseudotime color scale (0-50), showing a curved trajectory. Image C uses the same scatter plot, colored by group (EMS and Normal). Image D colors the scatter plot by state (1-5). Image E presents four scatter plots of gene expression versus pseudotime with fitted curves: CRYZ rises then falls, FRMD4B rises then falls, FZD6 trends downward and TPST1 declines early and remains low. Image F displays a circle network of cell types (Fibroblasts, Epithelial cells, etc.) with connected nodes. Image G shows a similar network with dense connections. Image H is a directed interaction diagram centered on Epithelial cells, linking to other cell types. Trajectory plots B-D define pseudotime and states used in E, while networks F-H summarize interactions among cell types. A multi-graph figure showing epithelial cell pseudotime, biomarker expression and cell communication networks. Expression dynamics of biomarkers in epithelial cells and intercellular communication. ( A ) Box plots showing the expression levels of the four biomarkers across different cell types. Epithelial cells were identified as the key cell type. ***P<0.001; ns, not significant. ( B ) UMAP plot of epithelial cells from the scRNA-seq dataset. ( C ) Pseudotime trajectory of epithelial cells, colored by differentiation state. ( D ) Pseudotime trajectory of epithelial cells, colored by the inferred pseudotime. ( E ) Expression dynamics of the four biomarkers along the pseudotime trajectory of epithelial cell differentiation. ( F ) The number of inferred interactions between epithelial cells and other cell types. ( G ) The strength of inferred interactions between epithelial cells and other cell types. ( H ) Circle plot illustrating the comprehensive cell-cell communication network between major cell types, with a highlight on the significant interactions between epithelial and endothelial cells.

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 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. 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. 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. 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. 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. 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, r 2 = 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 (r 2 .exposure > r 2 .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 ( \documentclass[12pt]{minimal} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \DeclareFontFamily{T1}{linotext}{} \DeclareFontShape{T1}{linotext}{m}{n} {linotext }{} \DeclareSymbolFont{linotext}{T1}{linotext}{m}{n} \DeclareSymbolFontAlphabet{\mathLINOTEXT}{linotext} \begin{document}$F = {\beta ^2}/S{E^2}$\end{document} ). 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. 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. 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. 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_r 2 = 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. 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. 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. 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. 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 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.

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.794 53 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.