The active components of BSHQ were initially screened from the TCM Systems Pharmacology (TCMSP) database (https://www.tcmsp-e.com/tcmsp.php) using the criteria of oral bioavailability ≥20% and drug-likeness ≥0.10. Following this, the corresponding protein targets of the selected active components were identified within the TCMSP database. These protein targets were subsequently normalized using the UniProt database (http://www.uniprot.org) by selecting the ‘Homo sapiens’ taxonomy. Finally, Cytoscape 3.10.1 (https://cytoscape.org/) was employed to construct a network that illustrates the relationships between the active components and their target genes. Topological analysis of the network was performed using the NetworkAnalyst plugin (https://www.networkanalyst.ca/) embedded in Cytoscape to calculate and rank the degree of connectivity of each active ingredient.

The datasets GSE7305, GSE11691 and GSE58178 were sourced from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) database and normalization analyses were subsequently performed (16-18). The GEO datasets used in the present study, GSE7305, GSE11691 and GSE58178, were also included in our previous study on hypoxia-related mechanisms in endometriosis (15); however, the present study was independently designed to investigate BSHQ-related pharmacological mechanisms and quercetin-associated therapeutic targets. All data preprocessing, gene screening, downstream analyses and figure generation were independently performed for the current study. The GSE7305 dataset includes 10 normal endometrial samples (paired eutopic endometrium from the same patients) and 10 endometriosis samples, and was analyzed using the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. The GSE11691 dataset consists of nine endometrial samples and nine matched ectopic endometriotic lesion samples, which are paired tissues collected from the same nine women with endometriosis. All samples were analyzed with the GPL96 [HG-U133A] Affymetrix Human Genome U133A Array. The GSE58178 dataset contains six normal endometrial samples (from healthy women without endometriosis) and six endometriosis samples, analyzed on the GPL6947 Illumina HumanHT-12 V3.0 expression beadchip. Endometriosis-related gene sets were retrieved from GeneCards (https://www.genecards.org/) using ‘endometriosis’ as the search key word. The GeneCards relevance score is an integrated metric automatically computed by the database based on multi-source evidence including disease-gene annotations, peer-reviewed literature, genome-wide association study data and mutation records, where a higher score indicates a stronger disease-gene association. Genes with a relevance score >1 were screened as candidate endometriosis-related genes, yielding a total of 1,359 genes. The GSE7846 dataset served as an independent validation cohort to verify the diagnostic performance of the identified hub genes. This dataset was also retrieved from the GEO database, containing 5 samples of primarycultured human endometrial endothelial cells (HEECs) isolated from the eutopic endometrium of diseasefree premenopausal women and 5 HEEC samples from women with endometriosis, analyzed with the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. This dataset was initially reported in a previous study exploring the transcriptomic profile of endometrial endothelial cells in endometriosis (19).

The DEGs between endometriosis samples and normal endometrial control samples analysis was conducted using the limma package (version 3.66.0; https://bioconductor.org/packages/limma/) in R software (version 4.3.1; https://www.r-project.org/), which utilizes a linear model to compute fold change (FC) values and P-values. The Benjamini-Hochberg method was employed to control the false discovery rate. Genes with |log2FC| >0.3 and P<0.05 were considered significantly differentially expressed. To visualize the expression profiles of the identified DEGs, heatmaps and volcano plots were generated using the pheatmap (version 1.0.13; https://cran.r-project.org/package=pheatmapcran) and ggplot2 (version 4.0.3; https://ggplot2.tidyverse.org) packages, respectively, in R software. Principal component analysis (PCA) was additionally performed to visualize the overall distribution of samples after normalization, and to assess the separation tendency between normal and endometriosis samples. Receiver operating characteristic (ROC) curve analysis was performed to assess the diagnostic performance of the key genes in differentiating ectopic from eutopic endometrial tissues. ROC curves were plotted using the pROC package (v1.18.0; https://cran.r-project.org/package=pROC) in R (v4.3.1) and the area under the curve (AUC) was calculated to evaluate the predictive accuracy of the model.

Chromosomal localization analysis of candidate genes was performed using the RIdeogram R package (version 0.2.2; https://cran.r-project.org/web/packages/RIdeogram/index.html) to visualize the chromosomal distribution of the genes. GO and KEGG enrichment analyses were conducted to identify the biological functions and pathways associated with the DEGs. Using the R package clusterProfiler (version 4.8.3; https://bioconductor.org/packages/clusterProfiler/), DEGs were analyzed for enrichment in biological processes (BP), cellular components (CC), molecular functions (MF) and pathways. Pearson correlation analysis was performed to calculate the expression correlation among candidate genes. Correlation plots were generated using the corrplot R package (version 0.92; https://cran.r-project.org/web/packages/corrplot/index.html), with color intensity indicating the strength of correlation. Chord plots were constructed using the circlize R package (version 0.4.16; https://cran.r-project.org/web/packages/circlize/index.html) to visualize the pairwise correlation relationships. The protein-protein interaction (PPI) network of the candidate genes was constructed using the STRING database (version 12.0; https://string-db.org/), with the organism limited to Homo sapiens and a minimum interaction score of 0.4 as the inclusion threshold. The network was then imported into Cytoscape 3.10.1, and the top 10 hub genes were identified via the maximal clique centrality (MCC) algorithm using the cytoHubba plugin (https://apps.cytoscape.org/apps/cytohubba).

WGCNA was performed to identify gene modules that are significantly associated with endometriosis. The normalized gene expression matrix derived from GEO datasets served as the input, with low-variance genes excluded to improve network stability and minimize noise. Sample clustering was conducted using the flashClust R package (version 1.01-2; https://cran.r-project.org/web/packages/flashClust/index.html). Using the WGCNA R package (version 1.72-1; https://cran.r-project.org/package=WGCNA), pairwise gene correlations were calculated and subsequently transformed into an adjacency matrix. A soft-thresholding power of eight was selected according to the scale-free topology criterion, ensuring the construction of a robust and biologically relevant network. The topological overlap matrix (TOM) was then computed to assess gene connectivity and similarity within the network. Hierarchical clustering was employed to organize genes into distinct co-expression modules, with each module assigned a unique color for visualization purposes. Module eigengenes, which represent the first principal component of each module, were calculated to summarize the expression profiles of the modules.

The single-cell transcriptomic dataset GSE213216 for endometriosis was retrieved from the GEO database (20). All downstream analyses were performed using the Seurat R package (version 4.3.0; https://satijalab.org/seurat/). First, quality control was performed to filter low-quality cells. After quality control, the filtered gene expression matrix was normalized using the LogNormalize method, and highly variable genes were identified via the FindVariableFeatures function. PCA was then conducted based on highly variable genes for linear dimensionality reduction. Subsequently, the uniform manifold approximation and projection (UMAP) method was applied for non-linear dimensionality reduction and 2D visualization of all cells, with the first 20 principal components selected as input. Cell clustering was performed using the FindNeighbors and FindClusters functions, and cell type annotation was conducted based on canonical marker genes. The expression patterns of target genes across cell types were visualized via feature plots, violin plots and heatmaps generated using the Seurat R package.

Machine learning techniques, specifically least absolute shrinkage and selection operator (LASSO) and support vector machine-recursive feature elimination (SVM-RFE), were employed to identify key genes from the DEGs. LASSO regression was conducted using the glmnet package in R (version 4.1-7; https://cran.r-project.org/package=glmnet), aiming to select the most relevant genes by shrinking regression coefficients and excluding non-informative features. Genes with non-zero coefficients were retained for further analysis. Following this, SVM-RFE was applied to rank the selected genes according to their significance in classification. This method, implemented via the e1071 package (version 1.7-12; https://cran.r-project.org/package=e1071), recursively removes the least influential features to enhance the performance of the model. Cross-validation was utilized to verify the stability and reliability of the feature selection process.

A nomogram was constructed based on the expression levels of the three identified hub genes (PTGER3, SERPINE1 and VCAM1) using the rms R package (version 6.7-1; https://cran.r-project.org/web/packages/rms/index.html). Each biomarker was assigned a weighted point according to its regression coefficient, and the total point was calculated to predict the risk of endometriosis. DCA was performed to evaluate the clinical utility of the nomogram using the rmda R package (version 1.6; https://github.com/mdbrown/rmda), by calculating the net benefit of the model across different risk thresholds (21).

Small molecules were selected from the TCMSP database and their structures were prepared using AutoDockTools-1.5.7 (https://ccsb.scripps.edu/mgltools/). This preparation involves adding hydrogen atoms, calculating partial charges and exporting the structures as PDBQT-format files. Protein structures were obtained from the RCSB PDB database (version 2023; https://www.rcsb.org). These structures were processed with PyMOL (version 2.5.2; https://pymol.org/) and AutoDockTools-1.5.7, where water molecules were removed and hydrogen atoms were added. AlphaFold-predicted structures (version 4; https://alphafold.ebi.ac.uk/) underwent similar processing. Docking simulations were carried out using AutoDock Vina (version 1.2.5; https://vina.scripps.edu/), generating nine conformations for each ligand. The binding affinity was assessed based on the docking score, where values <-7 kcal/mol were considered indicative of strong binding. The conformations with the lowest binding free energy were selected for subsequent analysis. Non-covalent interactions were examined using Protein-Ligand Interaction Profiler (PLIP; version 2.1.7; https://plip-tool.biotec.tu-dresden.de/), and visualized with PyMOL. The docking models with a binding free energy of <-7 kcal/mol were chosen for further study.

GSEA was performed to identify enriched biological pathways and MF terms associated with DEGs. The normalized gene expression matrix of the merged GEO datasets were ranked based on their association with the phenotype and used to generate a ranked gene list. Predefined gene sets from the Molecular Signatures Database (version 7.4; https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) or custom gene sets were selected for analysis. GSEA was conducted using GSEA software (version 4.3.2; https://www.gsea-msigdb.org/gsea/index.jsp) and an enrichment score was calculated for each gene set to determine its significance. Statistical significance was assessed through permutation testing. Enrichment score plots and heatmaps were generated to visualize the results.

The VK2/E6E7 and 12Z endometriosis cell lines, obtained from Procell Life Science & Technology Co., Ltd., were cultured in high-glucose Dulbecco's Modified Eagle's Medium (cat. no. PM15021; Procell Life Science & Technology Co., Ltd.), supplemented with 10% fetal bovine serum (FBS; cat. no. 10099141C; Gibco; Thermo Fisher Scientific, Inc.). The cells were grown in a humidified incubator at 37˚C with 5% CO₂ to ensure optimal growth conditions and maintain their viability for experimental purposes.

Transfection was performed on both 12Z and VK2/E6E7 endometriosis cell lines. Cells were seeded at a density of 2x10⁵ cells/well in 6-well plates 24 h prior to transfection, reaching 70-80% confluence on the day of transfection. Prior to transfection, 1.5-ml microcentrifuge tubes were UV-irradiated and Opti-MEM medium (cat. no. 31985070; Gibco; Thermo Fisher Scientific, Inc.) was pre-warmed to 37˚C. Cells were washed twice with PBS and the culture medium was then replaced with Opti-MEM. The transfection mixture was prepared by combining 125 µl Opti-MEM, 2.5 µg SERPINE1 overexpression plasmid and 4 µl Lipo 8000 reagent (cat. no. C0533; Beyotime Biotechnology). The mixture was added dropwise to the wells and incubated at 37˚C in a humidified 5% CO₂ incubator for 6 h. After incubation, the medium was replaced with fresh medium containing 10% FBS. For SERPINE1 overexpression experiments, the plasmid backbone used was pcDNA3.1(+), an ampicillin-resistant vector purchased from Beijing Tsingke Biotech Co., Ltd. Data S1 provides the complete coding sequence (consensus coding sequence region) of the human SERPINE1 gene (transcript variant NM_000602.5), which was used as the insert fragment for plasmid construction, the complete plasmid map is shown in Fig. S1. The negative control group was transfected with the corresponding empty pcDNA3.1 vector, which was obtained from the same supplier and also carried ampicillin resistance. For cell proliferation rescue assays, quercetin (cat. no. HY18085; MedChemExpress) treatment at the respective half-maximal inhibitory concentration (IC₅₀) concentration for each cell line was initiated 24 h after transfection. Cells were incubated at 37˚C in a humidified 5% CO₂ incubator, and proliferation was measured at 0, 24, 48, 72 and 96 h after quercetin administration.

Total RNA was extracted from cells treated with quercetin at the respective IC₅₀ concentration for each cell line for 24 and 48 h at 37˚C in a 5% CO₂ humidified incubator (only data obtained at the 24 h timepoint are presented in the current study), with DMSO as the vehicle control or SERPINE1-overexpression transfected cells using TRIzol^® reagent (Invitrogen; Thermo Fisher Scientific, Inc.), following the standard RNA extraction protocol, and cDNA was synthesized using the PrimeScript RT Reagent Kit (Takara Bio, Inc.) in accordance with the manufacturer's recommended protocol. qPCR was performed on the Applied Biosystems 7500 Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.) using 2X SYBR Green qPCR Master Mix (cat. no. G3320-15; Wuhan Servicebio Technology Co., Ltd.). Each 20-µl reaction contained 4 µl 2X SYBR Green Mix, 0.2 µl SERPINE1 forward primer (25 µM, 5'-CGCCTCTTCCACAAATCAGACG-3') and 0.2 µl SERPINE1 reverse primer (10 µM, 5'-GGCAGTTCCAGGATGTCGTAGT-3'), or 0.2 µl GAPDH forward primer (25 µM, 5'-ACCCACTCCTCCACCTTTGACG-3') and 0.2 µl GAPDH reverse primer (25 µM, 5'-TCTCTTCCTCTTGTGCTCTTG-3'), 0.5 µl cDNA and nuclease-free water to adjust the final volume. The SERPINE1 target gene and GAPDH housekeeping gene were detected in independent reaction tubes, with only one pair of primers added per reaction. The cycling conditions were set as follows: Initial denaturation at 95˚C for 10 min, followed by 40 cycles at 95˚C for 15 sec and 60˚C for 1 min. Relative expression levels of SERPINE1 were normalized to GAPDH using the 2^-ΔΔCq method (22).Western blotting (WB). For WB, cells were lysed in radioimmunoprecipitation assay buffer (cat. no. G2002; Wuhan Servicebio Technology Co., Ltd.) supplemented with 1% phenylmethylsulfonyl fluoride, and protein concentrations were measured by bicinchoninic acid assay (cat. no. G2026; Wuhan Servicebio Technology Co., Ltd.). Equal amounts of protein (30 µg/lane) were separated by SDS-PAGE on 10% gels and transferred to polyvinylidene fluoride membranes. After blocking with 5% (w/v) non-fat skim milk (cat. no. D8340; Beijing Solarbio Science & Technology Co., Ltd.) for 1 h at room temperature, the membranes were incubated with anti-SERPINE1 (1:1,000; cat. no. ab28207; Abcam) and anti-GAPDH (1:5,000; cat. no. ab8245; Abcam) primary antibodies overnight at 4˚C, followed by horseradish peroxidase-conjugated secondary antibodies (1:5,000; cat. nos. ab205718 and ab205719; Abcam) for 1 h at room temperature. Protein bands were detected using an enhanced chemiluminescence kit (cat. no. G2014; Wuhan Servicebio Technology Co., Ltd.) and semi-quantified by densitometric analysis using ImageJ software (version 1.53t; National Institutes of Health).

Cell proliferation was assessed using the CCK-8 assay kit (cat. no. G4103; Wuhan Servicebio Technology Co., Ltd.). Cells were seeded in 96-well plates and treated with quercetin at gradient concentrations (0, 1, 2, 4, 8, 16, 32 and 64 µM) for 24, 48, 72 and 96 h at 37˚C in a 5% CO₂ humidified incubator, with an equal volume of dimethyl sulfoxide (DMSO) as the vehicle control. At each time point, 10 µl CCK-8 solution was added to each well and incubated for 2 h. The absorbance at 450 nm was measured using a microplate reader. The resulting optical density (OD) values were used to quantify cell proliferation and determine the effects of the treatments on cell proliferation.

The IC₅₀ of quercetin on endometriosis cell lines was determined using the CCK-8 cell viability assay. Cells in the logarithmic growth phase were digested, resuspended in complete medium and seeded into 96-well plates at a density of 5x10³ cells/well in 100 µl complete medium. After 24 h of incubation for cell adhesion, quercetin was added at gradient concentrations (0, 1, 2, 4, 8, 16, 32 and 64 µM), with six replicate wells for each concentration. A blank control with only culture medium was set simultaneously. After 24 h of drug treatment, the absorbance at 450 nm was measured to assess cell viability, and IC₅₀ values were calculated via dose-response curve fitting.

Cell proliferation was measured using the EdU assay kit (cat. no. G1601; Wuhan Servicebio Technology Co., Ltd.). Cells were divided into the DMSO solvent control group and quercetin treatment group: cells were treated with quercetin at the respective IC₅₀ concentration for each cell line, and incubated at 037˚C in a 5% CO₂ humidified incubator for 24 h; the control group was treated with an equal volume of DMSO matching the solvent concentration of the quercetin treatment group. After treatment, cells were incubated with 10 µM EdU for 2 h at 37˚C, followed by fixation with 4% paraformaldehyde for 15 min at room temperature. Cells were then permeabilized using 0.5% Triton and stained according to the manufacturer's instructions, with Hoechst 33342 used for nuclear counterstaining. The EdU-positive cells were visualized and quantified under a fluorescence microscope. The percentage of EdU-positive cells was calculated to assess the proliferative activity of the treated cells.

CETSA was performed to assess the stability of SERPINE1 in the presence of quercetin. Briefly, 12Z and VK2/E6E7 cells were treated with 20 µM quercetin for 4 h at 37˚C; an equal volume of DMSO was added to the DMSO control group. Following treatment, cells were harvested and resuspended in lysis buffer. The cell suspensions were subjected to gradient temperatures (37, 42, 47, 52, 57 and 62˚C) for 3 min, followed by rapid cooling on ice. After three freeze-thaw cycles, cell lysates were collected and analyzed by WB to measure the levels of SERPINE1. The degree of SERPINE1 stabilization was determined by comparing the protein levels across different temperature treatments.

SPR analysis was conducted to evaluate the interaction between quercetin and SERPINE1. Recombinant human SERPINE1 protein (cat. no. ab282405; Abcam) was immobilized on a CM5 sensor chip (cat. no. BR100530; Cytiva), and increasing concentrations of quercetin (0.3125, 0.625, 1.25, 2.5, 5 and 10 µM) were injected over the surface at 25˚C, with an association phase of 120 sec and a dissociation phase of 300 sec to allow binding. The changes in resonance units were monitored in real-time using a Biacore T200 system (Cytiva) to assess the binding affinity and kinetics. The dissociation constant (KD) was calculated from the binding curves to determine the strength of the interaction between quercetin and SERPINE1. Data analysis was performed using dedicated software (Biacore T200 Evaluation Software; version 3.2; Cytiva) to obtain precise binding parameters.

Statistical analysis was performed using GraphPad Prism 9.0 software (Dotmatics) and R software. Continuous data were first assessed for normality using the Shapiro-Wilk test. Normally distributed data are presented as the mean ± standard deviation, whereas non-normally distributed data are presented as median (interquartile range). For comparisons between two groups, unpaired Student's t-test was used for normally distributed data with equal variances, and unpaired Welch's t-test was applied when variances were unequal. P<0.05 was considered to indicate a statistically significant difference.

In the present study, a total of 32 bioactive ingredients were identified from BSHQ. Among these, 12 ingredients were sourced from Paeoniae Radix Alba, while 22 originated from Hedysarum multijugum Maxim. Notably, kaempferol (MOL000422) and mairin (MOL000211) were found to be common to both herbs (Table SI). After eliminating duplicates, the 32 active ingredients were linked to 231 unique target genes. To explore the interactions between these active ingredients and their corresponding targets, the present study constructed a detailed interaction network using Cytoscape 3.10.1. This network visually represents the connections between the bioactive ingredients and their respective target genes. In the diagram, the orange rectangular nodes at the top denote the two constituent herbs of BSHQ, ovals represent the active ingredients and the green rectangular nodes in the inner layer symbolize the target genes. This layout provides a clear and comprehensive overview of the interactions among these herbal ingredients and their targets (Fig. 1). Subsequent analysis within NetworkAnalyst (embedded in Cytoscape) was carried out to identify the top five active ingredients with the highest degree of connectivity in the network. The ingredients with the highest connectivity were identified as quercetin, kaempferol, 7-O-methylisomucronulatol, formononetin and β-sitosterol (Table I).

The GEO datasets GSE7305, GSE11691 and GSE58178 were retrieved from the GEO database and subsequently merged for further analysis. After normalization, boxplots of the processed data confirmed that the medians across all samples were comparable, indicating consistency in preprocessing steps (Fig. 2A-D). PCA was performed to visualize the global transcriptomic distribution of the normalized samples; a clear tendency toward separation was identified between the normal and endometriosis groups, supporting the presence of transcriptomic differences between the two tissue types (Fig. 2A-D). DEGs were identified using the criteria |logFC| >0.3 and P<0.05. The results of differential gene expression were visualized using a volcano plot (Fig. 2E).

In GSE7305, 2,977 genes were downregulated and 2,818 genes were upregulated. In GSE11691, there were 1,336 downregulated genes and 1,324 upregulated genes. In GSE58178, 1,295 genes were downregulated and 1,384 genes were upregulated. In the merged datasets, 1,389 genes were downregulated and 1,445 genes were upregulated (Fig. 2E). To further illustrate these findings, a heatmap was generated to display the top 10 most upregulated and downregulated genes from the merged dataset (Fig. 2F).

The gene expression profiles from 25 endometriosis samples and 25 normal endometrial samples (merged from GSE7305, GSE11691 and GSE58178 datasets) were used to construct co-expression modules using the WGCNA algorithm. Sample clustering was performed using the flashClust package. When the soft-thresholding power was set to 8, the scale-free topology fit index reached 0.85, optimizing the average connectivity of the network (Fig. 3A). Hierarchical clustering was performed to organize genes into co-expression modules, and the gene dendrogram with module color assignments is shown in Fig. 3E. Applying the dynamic tree cut algorithm, four distinct gene co-expression modules were identified. The eigengene adjacency heatmap was generated to visualize the correlation relationships among different modules (Fig. 3B). Subsequent analysis of the relationships between module eigengenes and clinical traits revealed that the turquoise module exhibited the strongest correlation with endometriosis, comprising 698 genes (Table SII; Fig. 3C). Furthermore, a strong positive correlation was observed between the module eigengenes and the genes within the turquoise module (Fig. 3D).

A total of 1,359 endometriosis-related genes were retrieved from the GeneCards database (with a relevance score >1). The effective ingredients of BSHQ, as analyzed through TCMSP, along with endometriosis-related genes from GeneCards, the turquoise module genes from WGCNA and DEGs from the merged datasets, were intersected, resulting in the identification of 16 BSHQ-related genes (Fig. 3F). These 16 genes were as follows: Cyclin dependent kinase inhibitor 1A (CDKN1A), claudin 4 (CLDN4), cytochrome P450 family 1 subfamily B member 1 (CYP1B1), estrogen receptor 1 (ESR1), Fos proto-oncogene (FOS), heme oxygenase 1 (HMOX1), heat shock protein family A member 5 (HSPA5), interleukin (IL)6, MET proto-oncogene (MET), matrix metallopeptidase 9 (MMP9), prostaglandin E receptor 3 (PTGER3), selectin E (SELE), SERPINE1, secretory leukocyte peptidase inhibitor (SLPI), tumor protein p53 (TP53) and vascular cell adhesion molecule 1 (VCAM1).

To explore the relationships between these 16 genes, the present study first performed chromosomal localization analysis (Fig. 4A). Following this, an enrichment analysis of the 16 genes was conducted. The top three most enriched items in each of the GO BP, CC and MF categories, and KEGG pathways were identified. For GO BP terms, the top three enriched items were ‘cellular response to oxidative stress’, ‘cellular response to chemical stress’ and ‘response to reactive oxygen species’. For GO CC terms, the top three were ‘caveola’, ‘plasma membrane raft’ and ‘collagen-containing extracellular matrix’. For GO MF terms, the top three were ‘RNA polymerase II general transcription initiation factor binding’, ‘oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen, reduced flavin or flavoprotein as one donor, and incorporation of one atom of oxygen’ and ‘core promoter sequence-specific DNA binding’. For KEGG pathways, the top three enriched items were ‘lipid and atherosclerosis’, ‘fluid shear stress and atherosclerosis’ and ‘endocrine resistance’ (Fig. 4B; Table SIII).

To further investigate the relationships between the expression levels of these 16 BSHQ-related genes, chord and correlation plots were constructed (Fig. 4C and D). In the correlation plot, the color intensities or color scale were used to indicate the strength of the correlation, with blue and red representing negative and positive correlations, respectively. Additionally, using the STRING database, a PPI network for the 16 BSHQ-related genes was constructed (Fig. 4E). The darker colors in the network indicate a higher influence on endometriosis. Fig. 4F displays the ranking of all 16 genes by degree. The top 10 hub genes (CDKN1A, ESR1, FOS, HMOX1, HSPA5, IL6, MMP9, SERPINE1, TP53 and VCAM) were identified using the MCC algorithm in the cytoHubba plug-in (Fig. 4G).

The present study analyzed scRNA-seq data from endometriosis (GSE213216), applying rigorous raw data preprocessing. After normalizing the data and performing PCA to reduce dimensionality, the present study used the UMAP method to visualize all cells, which were subsequently grouped into 20 distinct clusters. Finally, the expression of the 16 BSHQ-related genes across different cell types was visualized. The analysis revealed distinct clusters of immune cell types based on their gene expression patterns. UMAP clustering analysis identified separate clusters corresponding to different cell types. Notably, clear separation was observed for B cells, dendritic cells, endothelial cells, fibroblasts, lymphocytes, mast cells, myeloid cells, natural killer (NK) cells and T cells, with clear delineation between these cell populations (Fig. 5A). These clusters were defined by their unique transcriptional profiles, as indicated by the different colors representing the various cell types.

Next, the present study examined the expression of the 16 BSHQ-related genes, CDKN1A, CLDN4, CYP1B1, ESR1, FOS, HMOX1, HSPA5, IL6, MET, MMP9, PTGER3, SELE, SERPINE1, SLPI, TP53 and VCAM1. The feature plots (Fig. 5B) and violin plots (Fig. 5C) illustrated the differential expression of these selected marker genes across the various immune cell types. Several of these genes, such as CDKN1A, FOS and HSPA5, were highly expressed in nearly all cell types, whereas others, including CLDN4, MET, MMP9 and SELE, exhibited low expression levels across almost all cell types. This highlights the diversity in gene expression profiles within the immune cell populations. Notably, SERPINE1 was highly expressed specifically in fibroblasts.

Finally, a heatmap of gene expression across all cell types was generated, revealing distinct expression patterns for each gene (Fig. 5D). The heatmap confirmed the tissue- and cell-type-specific expression of the selected genes, offering valuable insights into the regulatory mechanisms that govern immune cell functions.

To identify key biomarkers for endometriosis, the present study applied a combination of LASSO regression and SVM-RFE algorithms for feature selection. The LASSO analysis was first conducted to identify the most predictive features, as illustrated in Fig. 6A and B, which shows the relationship between the fraction of deviance explained and the coefficients of each feature. As the model was regularized, the coefficients for each feature were progressively shrunk towards zero, with the optimal features emerging when the deviance explained by each feature plateaued. The SVM-RFE analysis, shown in Fig. 6C, was used to rank the features based on their importance. To refine the list of candidate biomarkers, the present study intersected the results from both the LASSO and SVM-RFE analyses. The Venn diagram (Fig. 6D) revealed three overlapping genes, which were selected as the most significant biomarkers for endometriosis: PTGER3, SERPINE1 and VCAM1.

To further assess the diagnostic potential of the identified hub genes (PTGER3, SERPINE1 and VCAM1), the present study developed a nomogram incorporating these biomarkers to predict the risk of endometriosis. The nomogram (Fig. 6E) visually represents the weighted contributions of each gene to the overall risk score. The individual scores for each biomarker were calculated based on their expression levels, with PTGER3, SERPINE1 and VCAM1 each contributing specific points corresponding to their respective expression levels. The total score, which is the sum of the individual gene scores, was then used to compute the linear predictor. The nomogram revealed a linear predictor ranging from -6 to 7, with higher scores associated with an increased predicted risk of endometriosis. The total points, ranging from 0 to 220, can be used to stratify patients based on their likelihood of developing endometriosis. A risk score between 0.1 and 0.9 was observed across the merged GEO dataset cohort, making this a reliable tool for risk stratification in clinical practice.

Moreover, the clinical utility of the nomogram was validated through DCA (Fig. 6F), which evaluated the net benefit of applying the nomogram at different risk thresholds. The DCA curve demonstrated that using the nomogram to predict endometriosis risk provided a notable net benefit over traditional strategies (such as treating all or no patients), particularly at higher risk thresholds. This indicated that the nomogram, incorporating PTGER3, SERPINE1 and VCAM1, may represent a valuable tool for clinical decision-making and personalized risk assessment in endometriosis.

To assess the diagnostic performance of the identified biomarkers, the present study generated ROC curves. In the merged datasets, PTGER3 achieved an AUC of 0.869, SERPINE1 had an AUC of 0.773 and VCAM1 showed an AUC of 0.939 (Fig. 6G), indicating that the three genes exhibited strong diagnostic potential for endometriosis. Boxplots showed that the expression levels of all three genes were significantly different between normal and endometriosis samples (P<0.05).

Finally, the present study validated the performance of these biomarkers in an independent validation cohort (GSE7846). The ROC curve analysis in the validation set (Fig. 6H) revealed that SERPINE1 maintained high diagnostic value with an AUC of 0.920, whereas PTGER3 and VCAM1 exhibited reduced AUC values of 0.560 and 0.480, respectively. Expression analysis further confirmed the significant differential expression of SERPINE1, whereas PTGER3 and VCAM1 did not show significant differences in the validation set. Therefore, the SERPINE1 gene was selected as the next research object.

Based on the aforementioned network pharmacology analysis of BSHQ active ingredients and their target genes, quercetin was identified as a common BSHQ-derived bioactive ligand that interacts with the proteins encoded by PTGER3, SERPINE1 and VCAM1. The molecular structure file of quercetin was retrieved from the TCMSP database and subjected to docking preparation. Molecular docking of the three proteins with the core active ingredient quercetin was conducted using AutoDock Vina, and the most optimal binding poses were selected. The resulting protein-ligand docking structures showed that the binding free energies of all the optimal poses were <-7 kcal/mol (Table II). Molecular interaction analysis of the docking poses was performed using PLIP and visualization was carried out with PyMOL.

In the binding pocket of PTGER3, the amino acids Q339, Q103, T107 and T106 form hydrophobic interactions with quercetin. Additionally, amino acids F140 and T136 form hydrogen bonds with quercetin (Fig. 7A). In the binding pocket of SERPINE1, the amino acids V32, V284, E283, T214, Y220 and R30 participate in hydrophobic interactions with quercetin, whereas N31, E283, D285 and G218 form hydrogen bonds with quercetin (Fig. 7B). In the binding pocket of VCAM1, the amino acids P120 and Y119 interact hydrophobically with quercetin and R123, D122, A144 and I88 form hydrogen bonds with quercetin (Fig. 7C).

Subsequently, based on the median expression level of the SERPINE1 gene, the present study classified the endometriosis samples from the merged dataset into two groups: The SERPINE1 low-expression group and the SERPINE1 high-expression group. DEGs (|logFC| >0.3 and P<0.05) between the two groups were identified, followed by GSEA to further investigate the biological pathways associated with these genes. GSEA was performed to identify the significantly enriched biological pathways within the dataset. The analysis revealed several prominent pathways with varying levels of enrichment, as indicated by their normalized enrichment scores and corresponding statistical parameters. The enriched biological pathways included: ‘Focal adhesion PI3K-AKT-mTOR signaling pathway’, ‘Signaling by receptor tyrosine kinases’, ‘Signaling by interleukins’, ‘Nuclear receptors meta-pathway’, ‘Pathways in cancer’ and ‘Adipogenesis’ (Fig. 7D).

To evaluate the effects of quercetin on the proliferative capacity of the endometriosis cell lines VK2/E6E7 and 12Z, a IC₅₀ assay was performed.

The dose-response curve was constructed by plotting the relationship between the drug concentration and the percentage of cell proliferation. Through curve fitting, the IC₅₀ values for each cell line were calculated. The results demonstrated a dose-dependent decrease in cell proliferation as the drug concentration increased. For the 12Z cell line, the calculated IC₅₀ value was found to be 41.49 µM, whereas for the VK2/E6E7 cell line, the IC₅₀ value was 9.518 µM (Fig. 8A).

To further investigate the effects of quercetin on the proliferative capacity of the endometriosis cell lines VK2/E6E7 and 12Z, the present study employed both CCK-8 and EdU assays. The CCK-8 assay results revealed a significant reduction in cell proliferation following quercetin treatment, as evidenced by the substantial decrease in OD450 values for both cell lines. Specifically, for the 12Z cell line, the cell proliferation rate was significantly diminished at various time points, with a notable reduction observed at 24 h (t=3.553; P<0.01), 48 h (t=34.59, P<0.001), 72 h (t=42.95, P<0.001) and 96 h (t=45.88, P<0.001). Similarly, the VK2/E6E7 cell line demonstrated a pronounced decrease in cell proliferation across the same time intervals, with statistical significance observed at 24 h (t=44.60; P<0.001), 48 h (t=58.21; P<0.001), 72 h (t=126.2; P<0.001) and 96 h (t=24.73; P<0.001) (Fig. 8B). These findings indicated that quercetin exerts a potent inhibitory effect on the proliferation of both VK2/E6E7 and 12Z cells, with the effect becoming more pronounced over time.

To further evaluate the effects of quercetin on cell proliferation, the present study conducted an EdU fluorescence staining assay using both the 12Z and VK2/E6E7 cell lines. The results demonstrated a significant reduction in the EdU-positive cell rate after quercetin treatment, indicating a marked inhibition of cell proliferation (Fig. 8C). In the 12Z cell line, the EdU-positive cell rate was 6.234±0.249% in the DMSO control group, whereas it decreased to 4.692±0.241% in the quercetin-treated group, with a significant difference between the two groups (t=20.40; P<0.001). Similarly, for the VK2/E6E7 cell line, the EdU-positive cell rate in the DMSO control group was 9.708±1.337% and after treatment with quercetin, it decreased to 5.791±0.504%, with statistical significance observed (t=13.15; P<0.001). These findings confirmed that quercetin significantly inhibits the proliferative capacity of both 12Z and VK2/E6E7 cells, as evidenced by the reduced EdU incorporation and the decreased number of proliferating cells.

To investigate the effect of quercetin on SERPINE1 expression in 12Z and VK2/E6E7 cells, both RT-qPCR and WB were conducted at 24 and 48 h after treatment. The results of the RT-qPCR experiments demonstrated a significant reduction in SERPINE1 mRNA expression levels following quercetin treatment in both cell lines (Fig. 9A). In the 12Z cell line, SERPINE1 mRNA expression decreased from 1.000±0.015 in the DMSO control group to 0.497±0.023 in the quercetin-treated group, with a highly significant difference (t=31.35; P<0.001). Similarly, in the VK2/E6E7 cell line, the mRNA expression of SERPINE1 decreased from 1.000±0.031 in the DMSO control group to 0.894±0.033 in the quercetin-treated group, with a statistically significant reduction (t=4.089; P=0.01).

WB further confirmed these findings at the protein level (Fig. 9B). In the 12Z cell line, quercetin treatment led to a substantial reduction in SERPINE1 protein expression, which decreased from 1.000±0.102 in the DMSO control group to 0.522±0.061 in the quercetin-treated group, with a significant difference (t=6.955; P=0.002). Similarly, in VK2/E6E7 cells, SERPINE1 protein expression decreased from 1.000±0.093 in the DMSO control group to 0.653±0.015 in the quercetin-treated group, with a statistically significant reduction (t=6.401; P<0.01). These results collectively suggested that quercetin treatment may significantly inhibit both the mRNA and protein levels of SERPINE1 in 12Z and VK2/E6E7 cell lines, indicating its potential role in modulating the expression of key genes involved in the pathogenesis of endometriosis.

To evaluate the effect of quercetin on the stability of cellular targets, the present study conducted a CETSA to assess the interaction between quercetin and SERPINE1 (Fig. 9C). The results showed that SERPINE1 protein levels decreased progressively with increasing temperature. However, this temperature-induced degradation was markedly mitigated in the quercetin-treated groups, suggesting that quercetin may effectively stabilize SERPINE1 by protecting it from temperature-dependent degradation.

Based on the findings from both molecular docking and CETSA, the present study preliminarily confirmed that quercetin can bind to SERPINE1. To further validate this interaction, the present study performed SPR experiments. The SPR sensorgram curves in Fig. 9D recorded real-time response units across serially diluted quercetin concentrations (0.3125-10 µM). Fitting these response curves yielded a KD of 2.03x10^-6 M for the interaction between quercetin and SERPINE1, indicating a strong affinity between the two (Fig. 9D). Consistent with this binding effect, quercetin treatment suppressed SERPINE1 mRNA and protein levels, thus reducing its expression. These results collectively supported the conclusion that quercetin interacts with SERPINE1, potentially influencing its stability and function in endometriosis cells.

To assess the overexpression efficiency of SERPINE1 in 12Z and VK2/E6E7 cells, the present study transfected both cell lines with a SERPINE1 overexpression plasmid, and performed RT-qPCR and WB at 24 and 48 h post-transfection, respectively. The RT-qPCR results demonstrated a significant increase in SERPINE1 mRNA expression in the overexpression groups for both cell lines (Fig. 10A). In the 12Z cell line, the SERPINE1 mRNA expression increased from 1.003±0.095 in the DMSO control group to 2,8658.739±779.729 in the SERPINE1 overexpression group, with a highly significant difference (t=63.66; P<0.001). Similarly, in the VK2/E6E7 cell line, SERPINE1 mRNA expression increased from 1.000±0.019 in the DMSO control group to 9,115.191±1479.679 in the overexpression group, also showing a statistically significant increase (t=10.67; P<0.001).

WB further confirmed the successful overexpression of SERPINE1 at the protein level (Fig. 10B). In the 12Z cell line, SERPINE1 protein expression increased from 1.000±0.139 in the DMSO control group to 2.399±0.489 in the SERPINE1 overexpression group, with a significant difference (t=4.769; P=0.009). Similarly, in VK2/E6E7 cells, the SERPINE1 protein expression increased from 1.000±0.119 in the DMSO control group to 1.880±0.118 in the overexpression group, with a highly significant difference (t=9.126; P<0.001). These results indicated that transfection with the SERPINE1 overexpression plasmid effectively enhanced both mRNA and protein levels of SERPINE1 in 12Z and VK2/E6E7 cell lines, confirming the successful overexpression of SERPINE1 in these models.

To evaluate the effect of SERPINE1 overexpression on cell proliferation in the presence of quercetin, the present study transfected 12Z and VK2/E6E7 cells with a SERPINE1 overexpression plasmid and treated them with quercetin. Cell proliferation was assessed using the CCK-8 assay at multiple time points: 0, 24, 48, 72 and 96 h. The results from the CCK-8 assay indicated a significant increase in the OD450 values in the SERPINE1 overexpression + quercetin treatment groups compared with in the quercetin-only treatment groups, suggesting enhanced cell proliferation (Fig. 10C). In the 12Z cell line, the OD450 values at 24 h (t=63.28; P<0.001), 48 h (t=26.05; P<0.001), 72 h (t=6.011; P<0.001) and 96 h (t=13.03; P<0.001) were all significantly higher in the SERPINE1 overexpression group treated with quercetin. Similarly, in the VK2/E6E7 cell line, significant increases in cell proliferation were observed at all time points, with OD450 values at 24 h (t=65.96; P<0.001), 48 h (t=70.81; P<0.001), 72 h (t=37.70; P<0.001) and 96 h (t=32.22; P<0.001) showing a marked elevation in the SERPINE1 overexpression + quercetin treatment group compared with the quercetin-only treatment.

To further confirm these findings, EdU fluorescence staining was performed to assess cell proliferation. The results revealed a significant increase in the percentage of EdU-positive cells in the SERPINE1 overexpression + quercetin treatment groups compared with the quercetin-only treatment groups (Fig. 10D). In the 12Z cell line, the percentage of EdU-positive cells was 10.315±1.458% in the quercetin-only treatment group and increased to 13.952±1.882% in the SERPINE1 overexpression + quercetin treatment group, with a significant difference (t=7.325; P<0.001). Similarly, in the VK2/E6E7 cell line, the percentage of EdU-positive cells was 11.010±1.111% in the quercetin-only group and increased to 13.002±1.466% in the SERPINE1 overexpression + quercetin treatment group, also showing a significant increase (t=5.080; P<0.001). These findings collectively suggested that quercetin suppresses cell proliferation, while SERPINE1 overexpression reverses this inhibitory effect and restores cellular proliferative capacity in both 12Z and VK2/E6E7 endometriosis cells.

The results suggested that quercetin alone may suppress the proliferation of endometriosis cells and reduce endogenous SERPINE1 expression, whereas SERPINE1 overexpression counteracts this inhibitory effect. Specifically, SERPINE1 overexpression markedly increased cell proliferation in both 12Z and VK2/E6E7 cell lines in the presence of quercetin, as evidenced by the elevated OD450 values and increased EdU-positive cell rates at all time points.