Identification and verification of diagnostic biomarkers for deep infiltrating endometriosis based on machine learning algorithms

Abstract

This study addresses the challenges in the early diagnosis of deep infiltrating endometriosis (DIE) by exploring the potential role of the deubiquitinating enzyme USP14. By analyzing the GSE141549 dataset from the Gene Expression Omnibus (GEO) database, using bioinformatics methods and three machine learning algorithms (LASSO, Random Forest, and Support Vector Machine), the key feature gene USP14 was identified. The results indicated that USP14 is significantly upregulated in DIE and exhibits good predictive value (AUC = 0.786). Further analysis revealed the important role of USP14 in muscle function, cellular growth factor response, and maintenance of chromosome structure, and its close association with various immune cell functions. Immunohistochemical staining confirmed the high expression of USP14 in DIE tissues. This study provides a new molecular target for the early diagnosis of DIE, which holds significant clinical implications and potential application value.

Introduction

Endometriosis is a common gynecological disease characterized by the presence and growth of functional endometrial-like tissue outside the uterus, primarily in the pelvic peritoneum, ovaries, and rectovaginal septum, leading to pelvic pain and infertility [1, 2]. Endometriosis lesions are classified into three types: superficial endometriosis, ovarian endometriosis, and deep infiltrating endometriosis (DIE). DIE is the most aggressive phenotype among the three forms of endometriosis [3, 4]. Deep endometriosis is arbitrarily defined as endometriosis infiltrating the peritoneum > 5 mm. It is an external uterine adenomyosis, mainly manifested as single nodules larger than 1 cm in diameter, located in vesicular muscle folds or near the lower 20 cm of the intestine. In > 95% of cases, deep endometriosis is associated with pelvic pain, including dysmenorrhea, deep pain, severe chronic pain, headache, and difficulty in defecation [5, 6].

Current drug therapy for endometriosis is symptomatic treatment rather than cytoreductive therapy [7]. Medical treatment for endometriosis is primarily aimed at pain relief, with no single drug capable of resolving the disease. This treatment is essentially palliative, and the disease often recurs after drug withdrawal [8]. Interestingly, the growth of deep infiltrating endometriosis lesions can also be explained by genetic and epigenetic changes, which account for lesion development and risk [9]. Therefore, this study explores the immunity and treatment of DIE based on newly discovered biomarkers.

Ubiquitin-specific protease 14 (USP14) is a deubiquitinating enzyme (DUB) associated with the proteasome. It is a dual-domain protein that plays a regulatory role in proteasome degradation. USP14 belongs to the USPs family [10, 11]. It is known to be broadly involved in different typical cellular signaling pathways, including nuclear factor κB (NF-κB) and Wnt/β-catenin signaling pathways [12]. However, there are currently no studies reporting the role and function of USP14 in DIE. This article aims to reveal the immune-related research of USP14 in DIE and explore the biological functions of cells. This study primarily focuses on the expression of USP14 in DIE, immune infiltration, correlation analysis, and its expression in individual cells.

Data was downloaded and mined from the Genomics Expression Omnibus Database (GEO), and bioinformatics analysis was conducted to confirm the expression profile of USP14 in DIE. GO and KEGG enrichment analyses were performed to identify potential functions and pathways among differentially expressed genes in endometriosis. Machine learning, including supervised and unsupervised learning, can be applied to clinical datasets to develop risk models and redefine patient categories and genomic analysis of tumors [13]. Least Absolute Shrinkage and Selection Operator (LASSO), Random Forest, and Support Vector Machine (SVM) can all handle high-dimensional data such as gene expression matrices, helping us analyze biomarkers of diseases and advance disease research. This confirms that USP14 has good predictive value in DIE and explores the relationship between USP14 and immune cell infiltration, individual cells, and non-coding RNA. The protein expression level of USP14 was verified by immunohistochemical staining. Finally, we found that high expression of USP14 is correlated with the occurrence of DIE, and its expression is significantly associated with immune cells.

Materials and methods

Dataset download and differential analysis

We obtained the dataset GSE141549 from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) as our data source, which included sample information from 71 non-DIE patients and 77 DIE patients. Samples from non-DIE patients were derived from peritoneal lesions, while DIE patient samples were sourced from 28 intestinal, 24 sacral ligament, 22 rectovaginal, and 3 bladder origins. We read the gene expression data, converted the read data frame into a matrix format, and extracted the expression data portion. Subsequently, we used the avereps function to average duplicate samples to reduce data redundancy. After obtaining the processed expression matrix, we examined the data distribution characteristics to determine whether logarithmic transformation (log2) was necessary to stabilize variance and improve data distribution, thereby ensuring the accuracy and analyzability of the gene expression data. Ultimately, we identified 22 genes with corrected P-values < 0.05 and a 1-fold difference (|LogFC| < 0.585), all of which were upregulated genes in DIE.

Clinical samples and immunohistochemical staining

Immunohistochemical (IHC) analysis was used to validate the protein lesel. These tissues were enrolled between 2022 and 2023, with all diagnoses based on pathological and/or cytological evidence. Tissues were fixed in 4% formaldehyde and embedded in paraffin. The paraffin-embedded samples were sectioned into 6µm-thick slices, and then incubated with anti-human USP14 (HPA001308, Sigma). The results were visualized under a white light scanner (Pannoramic SCAN II, 3DHistech) and fluorescent scanner(NanoZommer S360, Hamamatsu), and the slices were sealed with xylene.

GO and KEGG enrichment analysis of DIE differential genes

Using comprehensive bioinformatics analysis methods, this study conducted an in-depth functional and pathway analysis of differentially expressed genes. First, we successfully converted gene names to Entrez IDs using the org.Hs.eg.db and org.Hs.egSYMBOL2EG packages, which is a necessary step for subsequent Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. Next, we performed GO and KEGG enrichment analysis separately using the enrich-GO and enrich-KEGG functions to explore the potential roles of these differential genes in biological processes and metabolic pathways. During the enrichment analysis, we set a screening criterion of a q-value threshold less than 0.05 to capture GO and KEGG entries with statistical significance, thereby ensuring the reliability of the analysis results. This step is crucial for revealing potential biological pathways and mechanisms. To visually present the enrichment analysis results and the relationship between genes and entries, we employed various graphical methods. Specifically, we drew bar charts and bubble charts to display the statistical information and gene proportions of enriched pathways, which helped quickly identify key biological processes and metabolic pathways. Additionally, we further analyzed the complex relationship between genes and entries through circle plots and chord plots. The circle plot showed the most significant GO and KEGG entries and their associated genes in a circular layout, while the chord plot visually displayed the interactions between genes and entries through connecting lines.

Identification of diagnostic markers using three machine learning algorithms

To avoid overfitting of the model, we first randomly divided all samples into a training set and a testing set in a 7:3 ratio. The training set samples participate in the feature selection and model training process. The test set participates in the performance evaluation of the model. To ensure the reproducibility and portability of the model, we obtained a set of independent validation data (GSE193928) from the GEO database. The validation data includes 37 samples: 19 control patients, 6 with adenomyosis, 8 with ovarian endometriosis, and 4 with deep infiltrating endometriosis. We applied the same normalization preprocessing to the independent validation data as we did to the training set. We then trained the model using the discovery cohort involved in this study and predicted the samples from the independent cohort. In the study of feature gene screening, we employed three powerful machine learning algorithms—LASSO, Random Forest, and SVM-RFE—to identify feature genes closely related to DIE from high-dimensional gene expression data. First, LASSO, as an advance-d dimension reduction technique, exhibited unique advantages in processing high-dimensional data [14]. We implemented LASSO analysis using the glmnet package and utilized 10-fold cross-validation to optimize the turning/penalty parameter, ensuring model accuracy and stability. This process not only helped us screen out genes that significantly impact the prediction of target variables (such as disease status) but also achieved effective dimension reduction. Next, we utilized the recursive feature elimination (RFE) function of the Random Forest algorithm to further explore gene importance. As an ensemble learning method, Random Forest improves overall prediction performance by constructing multiple decision trees and integrating their prediction results [15, 16]. During the RFE process, we assessed the predictive performance of the model through 10-fold cross-validation and determined feature genes based on each gene's relative importance in the model, typically measured by the reduction in Gini impurity. Specifically, we considered genes with a relative importance greater than 0.25 as feature genes, which played a crucial role in disease prediction. We tested 20 different machine learning models, including logistic regression, random forest, and XGBoost. By comparing the performance of each model on the training and testing sets, we ultimately evaluated the most robust model. Based on metrics such as accuracy, precision, recall, and F1 score, we finally chose SVM-RFE as the core algorithm for the diagnostic model. The performance of all involved models can be found from Figure S1. Finally, we employed the Support Vector Machine Recursive Feature Elimination (SVM-RFE) method for training diagnostic model. SVM-RFE combines the powerful classification ability of SVM with the feature selection strategy of RFE, enabling the selection of features that contribute most to classification and the removal of redundant features. Through 10-fold cross-validation, we ensured the stability and reliability of the SVM-RFE process [17, 18]. During the screening process, we assessed the impact of each gene on the performance of the SVM classifier and determined the feature gene set accordingly. We successfully screened out a series of feature genes closely related to DIE and took their intersection genes as key genes for DIE occurrence. To evaluate the effectiveness of these feature genes in disease diagnosis, we used the receiver operating characteristic (ROC) curve and area under the curve (AUC) as evaluation metrics. These metrics intuitively reflect the diagnostic performance and overall accuracy of the model at different thresholds, providing compelling evidence to support the potential application value of feature genes in disease diagnosis [19, 20].

Validation of feature gene expression in disease and data processing

We performed comprehensive preprocessing on the original gene expression dataset using the affy package in R, including background correction, normalization, and log2 transformation [21]. For cases where multiple probes corresponded to the same gene, we summarized the expression values using the average method to reduce noise and enhance accuracy. Subsequently, we conducted differential expression analysis using the limma package, identifying genes with significantly altered expression under different conditions by setting screening criteria of p-value < 0.05 and |log2 fold change (FC)| > 0.2. In specific discrimination analysis, we further tightened the criteria to require an FDR value < 0.05 and |logFC| > 0.585, ensuring that the screened differentially expressed genes (DEGs) were both statistically significant and biologically meaningful. Finally, we used visualization tools in R to display the expression patterns of DEGs, further validating the differential status of feature genes in DIE.

GSEA and GSVA analysis of co-expressed gene sets for signature genes

To calculate the correlation between genes, we utilized the 'cor' function to obtain the correlation matrix of gene expression data, visually demonstrating whether the expression patterns of genes are similar or opposite. Subsequently, we explored the biological significance within the gene expression data. First, we read and organized the gene expression data files, removing duplicate samples. Then, based on the median expression level of a target gene, the remaining samples were divided into high and low expression groups. Next, we calculated the average expression levels of other genes in these two groups and obtained their log fold changes (logFC) to identify genes with significant differences between the two groups. Gene Set Enrichment Analysis (GSEA) is a widely used technique in transcriptome data analysis that ranks gene lists in microarray studies using pre-defined gene set databases to identify significant and coordinated changes in gene expression data [22]. We employed GSEA to identify biological pathways significantly enriched in the high and low expression groups. By setting a threshold of p < 0.05, we filtered out the most significant enrichment results. Finally, to visually present these enrichment results, we plotted graphs of significantly enriched pathways in both the high and low expression groups. GSVA is a non-parametric unsupervised analysis method primarily used to evaluate the results of gene set enrichment in microarrays and transcriptomes. It assesses whether different metabolic pathways are enriched among samples by transforming the gene expression matrix between samples into a gene set expression matrix between samples [23, 24]. The GSVA package and its ssGSEA function were used to obtain the GSVA score for each gene set, representing the absolute enrichment of each gene set. The Limma package was used to compare the differences in GSVA scores between subtypes for each gene set.

Immune characteristic analysis of signature genes

To investigate immune-related functions in the gene expression data, we first read the gene expression information and then removed samples from the control group, retaining only data from the treatment group. Subsequently, we performed Single-Sample Gene Set Enrichment Analysis (ssGSEA) scoring on the preprocessed expression data using these gene sets to assess the activity of immune-related functions in each sample [25]. To more intuitively compare the scores between different samples, we normalized the scores to the same scale. Next, we divided the samples into high and low expression groups based on the median expression level of a target gene. Then, we combined the normalized immune function scores with the target gene grouping information.

To analyze the infiltration of immune cells in the samples, we used the CIBERSORT algorithm and conducted differential analysis between the disease group and the control group [26]. When executing CIBERSORT, we performed 1000 permutation tests to enhance the reliability of the results. After completing the calculation, we first filtered the results, retaining only immune cell types with P-values less than 0.05 to ensure the statistical significance of these results. To further understand the relationships between the infiltration of these different types of immune cells, we used the Pearson correlation coefficient to identify correlations among the differentially expressed immune cell types [27].

ceRNA network construction

We used the StarBase database (http://starbase.sysu.edu.cn/) to predict target genes (mRNA), upstream circRNA, and upstream lncRNA for signature miRNAs. We cross-referenced the predicted mRNAs, circRNAs, lncRNAs, and corresponding DE-mRNAs, DE-circRNAs, DE-lncRNAs using the miRanda, miRDB, and TargetScan databases, based on the ceRNA regulatory mechanism. Finally, we created a lncRNA (circRNA)-miRNA-mRNA network using "Cytoscape" software.

Single-cell RNA-seq data analysis

ScRNA-seq has revealed the expression patterns of the USP14 gene in different cell types, showcasing the breadth and complexity of its biological functions [28]. First, we downloaded and extracted information from GSE203191 in the GEO database. We used the Seurat package to perform rigorous data quality control (QC), including checking sequencing depth, removing low-quality cells (such as those with high mitochondrial gene expression or low gene detection counts), normalizing gene expression levels, and identifying and correcting potential batch effects [29, 30]. Subsequently, to visualize high-dimensional single-cell data in a low-dimensional space, we employed the t-Distributed Stochastic Neighbor Embedding (t-SNE) algorithm. The t-SNE algorithm creates a map of cell clusters in two-dimensional or three-dimensional space by calculating the similarity between cells, allowing similar cells to be close together and different cells to be far apart in the graph [31, 32]. This enabled us to effectively explore the expression patterns of the USP14 gene in different cell types or states and the distribution and relationships of these cells within the overall cell population.

Result

Screening of DEGs in DIE and non-DIE

Figure 1 showed a flowchart. In the GSE141549 dataset, differential expression analysis was conducted on 77 DIE samples and 71 NON-DIE samples. Among the 148 identified DEGs, 38 genes were significantly upregulated according to the screening criteria of |log2 Fold Change (FC)| > 0.585 and P-value < 0.05, while no genes were observed to be significantly downregulated. Heatmaps and volcano plots of DEGs between DIE and NON-DIE in the GSE141549 dataset were generated. The heatmap revealed that most DEGs were upregulated in the DIE group, particularly USP14, LRRFIP1, C20ORF127, HSPA1A, and FTHL12 (Fig. 2A). The volcano plot displayed 22 upregulated genes, among which 10 genes, including USP14, HSPA1A, FTHL12, USP49, A2M, LRRFIP1, C20ORF127, SEMA3C, RHOB, and NR4A2, were significantly upregulated (Fig. 2B).

Fig. 1

The flowchart of the study

Fig. 2

Differential Expression Analysis in DIE and non-DIE (A) Shows upregulated DEGs in DIE, especially USP14, LRRFIP1, C20ORF127, HSPA1A,and FTHL12. B Displays 22 upregulated genes; 10 are significantly elevated in DIE, including USP14, HSPA1A, and FTHL12

GO and KEGG enrichment analysis of differential genes

GO enrichment analysis indicated that the differential genes significantly impacted biological processes, influencing steroid hormone response, chemotaxis regulation, cellular response to external stimuli, response to tumor necrosis factor, DNA binding transcription activator activity, RNA polymerase II specificity, and related functions. Furthermore (Fig. 3A-C). KEGG pathway enrichment analysis revealed that these differential genes' signal pathways regulate human T-cell leukemia virus 1 infection, which may activate or disrupt immune response pathways and apoptosis pathways during viral infection (Fig. 3D-E).

Fig. 3

Functional Enrichment Analysis of Differential Genes (A-C) GO Enrichment: Shows BP, CC, MF functions and gene enrichment areas. D-E KEGG Enrichment: Highlights pathways regulating human T-cell leukemia virus 1 infection, immune response and apoptosis

Identification of USP14

Three different machine learning methods, the LASSO, Random Forest (RF), and SVM, were employed for decision-making and variable selection. By combining these methods, meaningful feature variables among the differential genes were identified, and the optimal diagnostic biomarker was selected using a Venn diagram. The LASSO method successfully identified seven meaningful feature variables: USP14, HSPA1A, FTHL12, USP49, C20ORF127, NR4A2, and HLA-DRB5(Fig. 4A-B). The RF method, an ensemble learning approach consisting of multiple decision trees, efficiently selected feature variables with an important score > 4, identifying six feature genes: USP14, HSPA1A, A2M, USP49, FTHL12, and LRRFIP1 (Fig. 4C-D). The SVM method identified three feature genes from the differential genes: FTHL12, USP14, and DUSP1(Fig. 4E-F). Ultimately, through the Venn diagram, USP14 and FTHL12 were selected as common feature genes and key molecules in DIE(Fig. 5A).

Fig. 4

Machine Learning Methods for Gene Identification. A-B LASSO Regression: Verifies accuracy via cross-validation curve. C-D Random Forest: Selects features based on importance scores. E-F SVM: Fits data and identifies key genes

Fig. 5

Key Gene Identification and Validation (A) Intersection of genes identified by three ML methods. B Line plots show gene expression trends. C-D Violin plots highlight FTLH12 & USP14 differences in DIE & non-DIE. E-F ROC AUC indicates prediction accuracy of USP14

Validation of USP14's predictive ability and accuracy in DIE

The AUC of the ROC curve was used to validate the diagnostic marker's accuracy in predicting disease survival. In the GSE141549 dataset, violin plots and line graphs visualized the expression levels of USP14 and FTHL12 in the DIE and non-DIE groups. Both plots revealed that USP14 expression was higher in the DIE group than in the non-DIE group (Fig. 5B-D). To further clarify USP14's true predictive ability in the disease prognosis model, we used the ROC curve to validate USP14's accuracy in DIE diagnosis. The results showed that the AUC of USP14 in DIE was 0.786, indicating that the USP14 model had high authenticity in predicting DIE prognosis (Fig. 5E-F). Combined with the ROC curve results, this suggests that the USP14 gene can accurately distinguish between DIE and non-DIE and has good value in prognosis prediction. In the independent validation cohort, we found that when only healthy controls and DIE samples were retained, the model could predict all 4 DIE samples with 100% accuracy, with no false positives. However, when we included samples with adenomyosis and ovarian endometriosis, the model incorrectly predicted one ovarian endometriosis sample as DIE, resulting in a false positive. The results indicate that our trained diagnostic model can accurately differentiate DIE patients from the healthy population. At the same time, we can still accurately identify all DIE patients even when we increase sample heterogeneity by introducing patients with other similar diseases (FigureS2).

GSEA and GSVA analysis of USP14 co-expressed gene sets

USP14 single-gene co-expression analysis aimed to reveal gene interactions and regulatory patterns by detecting expression correlations between genes, further annotating and predicting USP14's functions. We screened 177 USP14 co-expressed genes, including 58 negatively correlated genes such as WNT4, CGN, MAL2, TOP2A, and 119 positively correlated genes such as SFRP2, MYH11, ACTG2, DES (Fig. 6A-C). To better reveal the functions of these co-expressed genes within and outside cells, we conducted GSEA and GSVA enrichment analyses. GSEA revealed that USP14's co-expressed genes were significantly enriched in high-sample sets related to muscle contraction and muscle system processes, as well as the regulation of cellular response to growth factors. These findings highlighted their important roles in muscle function, structure, and cellular growth factor response. Additionally, these genes were associated with specific structures in myofibrils, such as the I-band, further indicating their crucial roles in muscle fine structure and functional regulation. This comprehensive analysis provided a new perspective for exploring the specific mechanisms of USP14 and its co-expressed genes in muscle biology. USP14's co-expressed genes were significantly enriched in the low-expression group and tightly related to biological processes such as chromosome organization, nuclear chromosome separation, and organelle fission. This discovery revealed that USP14 and its related genes may play key roles in maintaining chromosome structure, regulating cell division, and ensuring proper organelle distribution, providing new research avenues for further exploring the fine regulation mechanisms of cell division and chromosome dynamics. GSEA showed significant enrichment in muscle system functions, regulation of cellular responses, and muscle fiber structure, relating to immune-related pathways, inflammatory mediator signaling pathways, calcium ion channels, cardiac muscle energy metabolism, apoptosis, and oxidative stress (Fig. 6D). GSVA functional enrichment analysis revealed significant differences in gene functions in cellular aldehyde synthesis and metabolism, aldosterone metabolism, nucleotide monophosphate metabolism, glucocorticoid biosynthesis, chromatin structure synthesis, and primary alcohol biometabolism processes. Subsequently, we analyzed the impact of differential genes on disease pathways through GSVA and found increased expression of differential genes in non-small cell lung cancer, valine/leucine/isoleucine biosynthesis, sulfur metabolism, lysosomes, and steroid biosynthesis (Fig. 6E-F). We further found that the biological functions primarily regulated by USP14 include enabling cysteine-type endopeptidase activity [33], endopeptidase inhibitor activity [34], protein binding [35], and peptidase activity. These functions suggest that USP14 is likely to induce DIE by intervening in endopeptidase activity.

Fig. 6

USP14 Coexpression Analysis. A-B Heatmap and volcano plot show coexpressed genes. C Correlation analysis with USP14 coexpressed genes. D GSEA reveals enriched functions and pathways. E-F GSVA highlights metabolic & disease pathway impacts

Impact of USP14 on immune cells and their functions in DIE and non-DIE

Using CIBERSORT to calculate immune cell infiltration, we observed differences in the proportion of immune cell infiltration between the DIE and non-DIE groups. B cells and plasma cells accounted for a minimal proportion in both diseases, while mature dendritic cells, T cells, NK cells, and macrophages comprised a larger proportion (Fig. 7A). Furthermore, we evaluated the correlation among the above 21 immune cell types and calculated correlation coefficients. Most cells showed a negative correlation, with the highest negative correlation coefficient found in one group, including activated natural killer cells (correlation coefficient R = 0.37) (Fig. 7B). Additionally, two groups of cells had a relatively pronounced positive correlation: M1 macrophages (R = 0.43) and memory B cells (R = 0.35), with M1 macrophages having the highest correlation coefficient with USP14. To further investigate USP14's impact on 22 immune cell types, we used dumbbell plots and correlation scatter plots to display the correlation between immune cells and feature genes. M1 macrophages (p = 0.004) and memory B cells (p = 0.023) showed a significant positive correlation with USP14, while activated natural killer cells (p = 0.015) exhibited a significant negative correlation (Fig. 7C-D). In the differential analysis of immune cells, box plots show that cytotoxic activity and NK cells are highly expressed in non-DIE conditions, while mast cells are highly expressed in DIE conditions (Fig. 7E). Based on this, after scoring the relationship between USP14 and 29 immune functions, we found that USP14 played a role in pro-inflammatory responses, mast cells, and NK cell functions. USP14's role was statistically significant in cell lytic activity (crucial in immune responses for clearing pathogens and abnormal cells), immature dendritic cells, and type II interferon function (Fig. 7F).

Fig. 7

Immune Function Correlation Analysis. A Immune cell infiltration and percentage in DIE and non-DIE. B Heatmap of immune cell correlations. C-D Correlation coefficients between immune cells and DIE. E Immune cell difference analysis. F Analysis of immune-related functions affected by USP14

ceRNA regulatory network and single-cell RNA-seq data analysis

Using RNA-seq data from the GEO database, single-cell transcriptome sequencing revealed the expression patterns of the USP14 gene in different cell types, showcasing its biological functions' breadth and complexity. Through color distinction, we could visually see the distribution of T cells, B cells, and other cell types including fibroblasts, NK cells, endothelial cells, smooth muscle cells, epithelial cells, and monocytes, highlighting their expression differences and potential biological function distinctions (Fig. 8A-B). The ISNE_1 axis and cell size information provided additional dimensions for analyzing cell distribution patterns, further enhancing our understanding of cell type diversity (Fig. 8C). We understood that USP14 was expressed at different levels in various cell types, reflecting its activity levels in different cells (Fig. 8D-E). To further validate USP14's diagnostic and predictive value in DIE, we constructed a three-dimensional regulatory network consisting of DIE lncRNA, miRNA, and mRNA using bioinformatics methods. From the ceRNA regulatory network, we found 16 miRNAs and 32 lncRNAs interacting with USP14 mRNA (Fig. 8F).

Fig. 8

Single-cell Data and ceRNA Network Analysis. A-D Cellular expression of USP14. E ceRNA-related networks interacting with USP14

Immunohistochemical staining of USP14 in endometriosis

We first screened DIE and non-DIE tissue sections containing endogeneity lesions, and determined the level of USP14 expression in two different tissue sites by immunohistochemical staining. The results further confirmed the existence of high level of USP14 expression in DIE tissues. In IHC, it can be seen that USP14 is highly expressed in glandular epithelial cells of ectopic lesions, and its expression level is significantly higher than that in surrounding cells, no matter in deep endometriosis or non-deep endometriosis (Fig. 9). The difference is that in deep endometriosis, the expression of USP14 in other cells except the focal glandular epithelium is significantly higher than that in non-deep endometriosis.

Fig. 9

Immunohistochemistry of DIE and non-DIE clinical samples reveals significantly higher USP14 expression in DIE tissues compared to peritoneal tissues

Discussion

Endometriosis is a debilitating disease with chronic inflammatory characteristics. It is estimated that 10–15% of women of reproductive age suffer from pelvic endometriosis, and the recurrence rate among treated patients ranges from 5% to > 60% [36, 37]. So far, none of the proposed pathogenic theories (retrograde menstruation, metaplasia, and remnants of the Müllerian system) can explain all different types of endometriosis. Current research demonstrates that there is no relationship between the extent of the disease and its symptoms. There is no blood test available for diagnosing endometriosis, and little is known about its true prevalence, distribution in the population, or risk factors [38]. For example, a reduction in the mean maximum ovarian endometriosis diameter and a 30% decrease in the need for surgical treatment were observed after 12 months in the dienogest and estrogen-progestin groups. Long-term dienogest treatment has a greater effect on alleviating menstrual pain and discomfort. Progestins and low-dose oral contraceptives are ineffective in one-third of symptomatic women worldwide. When first-line drugs are ineffective, oral gonadotropin-releasing hormone (GnRH) antagonists are an effective and tolerable treatment option to optimize and personalize endometriosis care [39]. Many molecular differences between endometriotic lesions and normal endometrium pose difficulties in the development of new drug therapies and treatment approaches. Surgery remains the gold standard for definitive diagnosis but must be balanced against the risks of surgical complications and potential reduction in ovarian reserve [40]. Different types of endometriosis have distinct biological and clinical features. Unlike superficial endometriosis, DIE often causes severe clinical symptoms, including significant pain and gastrointestinal and urinary tract abnormalities. For women with DIE who do not respond to medical treatment, surgical resection may require extensive surgical procedures, including segmental bowel resection and pelvic tissue resection [41]. The primary goals of treatment are to reduce pain, correct infertility, and potentially avoid or delay the onset of long-term endometriosis-related sequelae such as fibrosis, adhesions, and malignant transformation. Although advancements in technology (minimally invasive diagnostic tools, magnetic resonance imaging, high-resolution vaginal ultrasound, etc.) and a better understanding of the pathophysiology of endometriosis for the development of new treatment strategies are ongoing for the diagnosis and management of endometriosis, many controversial issues remain. It has been reported that many genetic, endocrine, immune, and other factors are significantly associated with the occurrence of endometriosis, but so far, the pathogenic factors of endometriosis have not been identified [42]. Therefore, this study mainly explores the value of biomarkers in the prevention and treatment of the most invasive deep infiltrating endometriosis (DIE).

In this study, we address the challenge of early diagnosis of DIE by employing machine learning algorithms and bioinformatics analysis to explore the value of the deubiquitinating enzyme USP14 as a potential biomarker. Currently, the diagnosis of DIE relies mainly on surgical pathological confirmation, and the lack of early non-invasive diagnostic methods greatly limits disease management and treatment. Through comprehensive analysis of the GSE141549 dataset in the GEO database, our study found that USP14 is significantly upregulated in DIE, suggesting its potential role in disease development. We used differential expression analysis to screen for genes significantly upregulated between DIE and non-DIE and further confirmed USP14 as a key feature gene through joint analysis using three machine learning algorithms: LASSO, Random Forest, and Support Vector Machine. ROC curve analysis showed that USP14 has good predictive value in DIE diagnosis (AUC = 0.786), verifying its potential as a biomarker. This finding provides a new molecular target for the early diagnosis of DIE and is expected to enable the development of non-invasive diagnostic methods in the future. Next, through GO and KEGG enrichment analysis, we revealed the important roles of USP14 co-expressed genes in steroid hormones, muscle function, cellular growth factor response, and chromosome structure maintenance. These results suggest that USP14 may participate in the development of DIE by regulating these key biological processes and pathways. The pathogenesis of endometriosis may be multifactorial, including anatomical, immune, inflammatory, hormonal (estrogen), oxidative stress, genetic, epigenetic, and environmental factors. Studies have shown that the evasion of endometrial cells from peritoneal immune surveillance contributes to the creation and maintenance of peritoneal endometriosis, but the specific mechanisms are not yet clear [43]. Therefore, our further immunocorrelation analysis indicated that USP14 is closely related to various immune cells and their functions, especially M1 macrophages and memory B cells. This finding not only deepens our understanding of the role of USP14 in cell biology but also provides a new perspective for explaining the complex pathological process of DIE. In particular, the regulatory role of USP14 in the immune microenvironment may have an important impact on the infiltration and progression of DIE, providing potential targets for future immunotherapy strategies. To gain a deeper understanding of the function of USP14 in DIE, we analyzed the expression pattern of USP14 in different cell types using single-cell RNA-Seq data. The results showed that USP14 is expressed in various cell types but with significant differences in expression levels. This finding reveals the cell-specific role of USP14 in DIE and helps us better understand its specific mechanisms during disease development. In this study, we found that M1 macrophages, memory B cells, and NK cells have the highest correlation coefficient with USP14. Numerous studies have also confirmed that these three types of cells are significantly related to the occurrence and progression of deep infiltrating endometriosis (DIE). M1NVs inhibit the development of endometriosis directly or by repolarizing macrophages from M2 to M1 phenotype. Therefore, the administration of M1NVs may represent a novel method for the treatment of endometriosis [44]. Shengnan Chen applied three endometriosis datasets and found that memory B cells were significantly associated with endometriosis in all three datasets [45]. Another group evaluated the relationship between Treg and NK cell-related cytokines in deep infiltrating endometriosis lesions and clinical symptoms of the disease. They eventually discovered that deep infiltrating rectosigmoid endometriosis displays alterations in Treg and NK cells [46]. These studies have confirmed, from different perspectives, the dysregulation of various immune cells, which may be related to the development of DIE. Therefore, on one hand, when abnormal levels of these immune cells are detected, it indicates a risk for the presence of DIE. On the other hand, if pharmacological intervention can adjust the abnormal levels of immune cells back to normal, it may have potential effects in alleviating DIE. Finally, through immunohistochemical staining, we further verified the high expression characteristics of USP14 in DIE tissues. The staining results showed that the expression level of USP14 in DIE tissues was significantly higher than that in peritoneal tissues. Simultaneously, the results of high expression in cells other than endoheterotopic lesions (including smooth muscle cells) further confirm the above views. This finding is not only consistent with our bioinformatics analysis results but also provides intuitive evidence for the clinical detection of USP14 as a DIE biomarker.

When exploring the potential role and function of the USP14 gene in deep infiltrating endometriosis (DIE), we should also pay attention to existing research on USP14 in other disease areas. As a deubiquitinating enzyme, USP14 plays an important role in various complex diseases such as cancer and neurodegenerative diseases [47, 48]. It has been widely studied for its regulatory role in multiple cellular signaling pathways, particularly in functions such as the NF-κB and Wnt/β-catenin pathways [49, 50]. For the invasiveness of DIE, the specific mechanism of USP14 in it has not been fully elucidated, but our study reveals that its high expression is closely related to the occurrence and development of DIE. In the field of oncology, USP14 has been reported to be closely related to the occurrence and development of various malignancies, promoting the proliferation and migration of tumor cells by regulating mechanisms such as protein degradation and cell cycle [51, 52]. These findings indicate that USP14 may act as a key regulatory factor involved in processes such as cell proliferation, apoptosis, and migration. This study, through the integrated use of machine learning algorithms and bioinformatics analysis methods, reveals the value of USP14 as a potential biomarker for DIE. Our research not only provides a new molecular target for the early diagnosis of DIE but also offers important clues for a deeper understanding of its pathogenesis. In the future, we will further explore the specific mechanism of USP14 in DIE and devote ourselves to developing early non-invasive diagnostic methods based on USP14, with the aim of improving the diagnosis and management of DIE patients. At the same time, targeted therapy against USP14 may also become a new strategy for treating other related diseases. We believe that with the deepening of research, USP14 will play an increasingly important role in the diagnosis and treatment of DIE and other related diseases, contributing new forces to improving patients' quality of life and increasing the cure rate of the disease.

Conclusion

This study has achieved significant progress in investigating the relationship between USP14 and DIE. It not only provides a new molecular target for the early diagnosis of DIE but also lays a solid foundation for further exploring the mechanism of USP14 in the disease and developing novel therapeutic approaches. Future research should continue to delve into the specific pathways of USP14 in DIE, explore its potential as a therapeutic target, and integrate research findings from other diseases to construct a more comprehensive disease network model. We believe that as research continues to deepen, USP14 will play an increasingly important role in the diagnosis and treatment of DIE and other related diseases, contributing new insights to improving patients' quality of life and enhancing disease cure rates.