Closing the evidence loop-membrane-lipid homeostasis and vesicular transport link DEHP exposure to endometriosis

We first obtained the molecular structure of DEHP from PubChem and performed comprehensive target prediction by integrating three complementary databases—ChEMBL, PharmMapper, and SwissTargetPrediction. Using a “union for recall, intersection for confidence” strategy with de-duplication, we identified 1,364 putative protein targets, thereby establishing a sufficiently broad search space for subsequent convergence with disease-side expression evidence while minimizing single-algorithm bias (Fig.  2 A, B). On the disease side, we integrated GSE51981 , GSE6364 , and GSE7305 and applied batch correction. PCA showed pronounced platform-driven clustering before correction (PC1 = 62.4%, PC2 = 11.9%; Fig.  3 A), whereas after correction the three datasets were well mixed and case–control separation improved (PC1 = 21.9%, PC2 = 8.8%; Fig.  3 B), indicating substantial attenuation of batch effects. Batch correction substantially attenuated platform effects and improved case–control separation, yet mild dataset imprint persists along secondary PCs. Heatmap columns were clustered by sample (Euclidean distance, complete linkage), and all downstream DE/WGCNA used the ComBat‑adjusted matrix. Using thresholds of |log₂FC| ≥ 1 and FDR < 0.05, we identified differentially expressed genes and observed coherent cross-sample expression patterns (volcano/heatmap in Fig.  3 C–D). We then applied WGCNA to stratify the whole-genome matrix into co-expression modules (Fig.  3 E). Among them, MEbrown ( r  = 0.58, P  = 3.6 × 10 − ¹¹) and MEblue ( r  = − 0.57, P  = 6.5 × 10 − ¹¹) were the top disease-associated modules by magnitude (Fig.  3 F). The negative sign for MEblue indicates a disease-linked downregulated program (higher eigengene in controls). Considering association strength, module size, and biological interpretability, we designated MEblue (2,537 genes) as the principal candidate module and intersected it with the DEG set (with direction retained) to obtain 229 high-confidence Ems-related genes (Venn: 305 DE-only genes; 229 intersection; Fig.  3 G). This provides “expression + network” dual evidence for subsequent convergence with DEHP targets while preserving the up/down direction for interpretation. Interpretation in words. The functional cascade proceeds from membrane‑lipid remodeling, through vesicular transport, to detoxification and de‑esterification, providing a continuous route by which DEHP exposure could influence lesion microenvironments. To bridge chemical exposure with disease evidence, we intersected the 1,364 DEHP candidate targets with the 229 Ems-related genes identified above and obtained 17 overlapping key targets (Fig.  4 A). The protein–protein interaction network indicated that these genes form a compact functional subnetwork, with core nodes concentrated in membrane-lipid synthesis/remodeling (ELOVL6, HMGCR, LYPLA1), detoxification/de-esterification (EPHX1, NCEH1), transport and vesicular dynamics (SLC1A5, VAMP2), and genome-stability surveillance (MSH2) (Fig.  4 B, C). GO/KEGG enrichment further pointed to terms closely related to lipid metabolism and vesicle trafficking, including bile secretion (KEGG), serine/carboxylate ester hydrolase activity, lipase activity, and endocytic/phagocytic vesicles and their membranes (Fig.  4 D). The 17-gene overlap delineates a functional corridor in which membrane-lipid remodeling enzymes adjust acyl-chain composition and membrane fluidity; vesicular-trafficking components then package and route lipid- and protein cargo within the lesion microenvironment; finally, detoxification and de-esterification machinery clear xenobiotics and lipid-peroxidation products. In the reverse direction, increased vesicular flux facilitates the delivery of detox enzymes and substrates, and the ensuing removal of reactive species feeds back to stabilize membrane-lipid homeostasis. Together, these steps define a coherent, exposure-responsive circuit rather than isolated pathways. Nodes such as ELOVL6/HMGCR/LYPLA1 map to membrane-lipid homeostasis, VAMP2 to vesicle fusion, SLC1A5 to nutrient/glutamine input, and EPHX1/NCEH1 to detoxification/de-esterification, providing a mechanistic corridor through which DEHP exposure could influence lesion microenvironments. Based on these 17 overlapping targets, we constructed multiple machine-learning models and evaluated them in training and external validation cohorts, achieving overall stable discrimination (most AUCs > 0.75; Fig.  5 A). Single-gene ROC analysis showed the highest cross-cohort transferability for UGT8 (AUC = 0.869) and EPHX1 (0.853), followed by VAMP2 (0.797), ELOVL6 (0.784), HMGCR (0.776), and SLC1A5 (0.685) (Fig.  5 B). The directions of differential expression were consistent with function: EPHX1 up-regulation reflects activation of the detoxification/de-esterification axis, whereas down-regulation of UGT8, ELOVL6, and LYPLA1 suggests suppression of sphingolipid synthesis and very-long-chain fatty-acid elongation (Fig.  5 C). To achieve mechanistic interpretability, we quantified feature contributions using SHAP; ELOVL6, LYPLA1, UGT8, SLC1A5, HMGCR, EPHX1, and VAMP2 ranked highest (Fig.  5 D), and the waterfall plot for a representative sample illustrated the cumulative effects of multiple features on the predicted probability (Fig.  5 E). Pairwise dependence plots revealed ELOVL6–UGT8 synergy, HMGCR–LYPLA1 opposing effects, and EPHX1–SLC1A5 context dependence, delineating a non-linear interaction spectrum related to lipid reprogramming and nutrient input (Fig.  5 F, M). Integrating model robustness, single-gene diagnostic performance, and SHAP explanations, we defined these seven genes as the core Ems signature associated with DEHP. SHAP exposes non-linear feature relationships underpinning classification: ELOVL6–UGT8 synergy suggests coordinated control of acyl-chain elongation and sphingolipid conjugation, whereas HMGCR–LYPLA1 antagonism implies opposing contributions from cholesterol biosynthesis and depalmitoylation-driven membrane relocalization. EPHX1–SLC1A5 context dependence links detox flux to amino-acid uptake, aligning with an immunometabolic read-out of exposure effects. We docked DEHP against seven proteins implicated by the core signature (UGT8, HMGCR, ELOVL6, SLC1A5, EPHX1, LYPLA1, VAMP2) using AutoDock Vina. Docking grids were centered on annotated/structure-informed pockets and expanded to include membrane-proximal grooves (per-target centers and sizes below). All targets yielded plausible best-pose affinities (ΔG range − 6.7 to − 3.6 kcal/mol), with a recurring geometry of hydrophobic embedding plus one to two polar clamps: the aromatic/alkyl scaffold of DEHP packs against hydrophobic/π surfaces, while the di-ester carbonyl oxygens are engaged by polar sites. To mitigate potential noise from the deliberately large boxes, we examined local pose clustering within the top ten modes (number of neighbors with RMSD ≤ 3.5 Å is reported per target) and list representative atom-level contacts. HMGCR (−6.2; 6/10; center (−14.1, 7.2, −7.2); size 109.2 × 106.7 × 114.7 ų): Polar: Ile536@O 3.08 Å, Gln814@NE2 3.15 Å; Hydrophobic/π: Tyr533 (3.15 Å), Pro813 (3.43 Å), Val538 (3.47 Å), Ile762 (3.48 Å), Leu811 (3.59 Å). The large box captures membrane-side grooves; top modes form a coherent local cluster (Fig.  6 A). EPHX1 (−5.3; 1/10; center (−12.5, 4.0, 3.4); size 93.6 × 82.0 × 73.6 ų): Polar: Ser317@OG 2.97 Å, Arg86@N 3.04 Å, Arg66@O 3.34 Å, Lys69@O 3.39 Å; Hydrophobic/aliphatic: Arg71 (3.58 Å), Val319 (3.59 Å). Polar residues align around the carbonyl oxygen in a detox-relevant pocket (Fig.  6 B). ELOVL6 (−6.1; 1/10; center (−2.5, −0.4, 9.2); size 69.1 × 56.9 × 79.6 ų): Polar: Tyr181@OH 2.92 Å, His141@ND1 3.32 Å; Hydrophobic/π: Phe249 (3.48 Å), Tyr180 (3.57 Å). Consistent with the hydrophobic tunnel of membrane-embedded elongases (Fig.  6 C). SLC1A5 (−5.7; 1/10; center (27.9, 3.0, −8.7); size 106.6 × 89.4 × 108.9 ų): Polar: Glu291@OE2 3.05 Å, Arg465@NH2 3.08 Å, Lys288@NZ 3.47 Å; Hydrophobic/π: Trp461 (3.32 Å), Ala433 (3.52 Å). A “polar clamp + hydrophobic walls” configuration that accommodates the di-ester (Fig.  6 D). VAMP2 (− 3.6; 3/10; center (− 15.6, − 7.9, 29.5); size 111.4 × 75.6 × 109.2 ų): Polar: Arg86@NE 2.91 Å, Lys87@N 3.46 Å; Aromatic/hydrophobic: Trp90 (3.51 Å). A shallow surface-proximal site consistent with membrane/interface binding (Fig.  6 E). LYPLA1 (− 4.9; 6/10; center (− 0.8, 9.7, − 0.4); size 58.5 × 68.9 × 53.1 ų): Polar: Thr8@OG1 2.90 Å, Asn95@O 3.31 Å, Leu99@N 3.31 Å; Hydrophobic: Leu99 (3.52 Å), Ile96 (3.59 Å). A typical interleaved polar–hydrophobic lipase pocket (Fig.  6 F). UGT8 (best − 6.7 kcal/mol; 7/10 neighbors ≤ 3.5 Å; grid center (2.8, 37.8, − 11.9); size 73.4 × 103.3 × 107.1 ų): Polar candidates: Tyr200@OH 2.71 Å, Tyr472@OH 2.89 Å, Gly169@O 3.30 Å; Hydrophobic/π: Phe473 (3.41 Å), Tyr200 (3.55 Å). An aromatic–hydrophobic groove stabilizes the ring system and carbonyl oxygens (Fig.  6 G). Taken together, DEHP consistently fits viable pockets across targets via hydrophobic packing and one to two polar anchors, dovetailing with the expression/network and interpretable-ML evidence along the lipid-metabolism–vesicular–detox axis. Large boxes were used to ensure coverage of membrane-proximal grooves; follow-ups (higher exhaustiveness, re-docking, short MD, and rescoring) can further test robustness without altering the present conclusion of structural feasibility. For the core proteins HMGCR, EPHX1, ELOVL6, SLC1A5, VAMP2, LYPLA1, and UGT8, we performed standardized molecular docking (receptor structures from PDB/AlphaFold; AutoDock Vina scoring; PyMOL visualization). All complexes adopted lowest-energy, plausible poses within accessible pockets, with DEHP forming hydrogen bonds of approximately ~ 2.2–3.3 Å alongside extensive hydrophobic/π interactions, exhibiting a “pocket fit + multi-point anchoring” pattern indicative of stable binding (Fig.  6 A, G). These results align with the expression/network/machine-learning evidence and suggest the structural feasibility of direct DEHP modulation of lipid-metabolism and membrane-associated proteins. To assess the temporal robustness of docking poses, we conducted 100-ns GROMACS simulations for the UGT8–DEHP, ELOVL6–DEHP, and HMGCR–DEHP complexes. RMSD trajectories indicated overall stabilization, with ELOVL6–DEHP showing the smallest fluctuations and earliest equilibration (Fig.  7 A), UGT8–DEHP exhibiting moderate fluctuations (Fig.  7 A), and HMGCR–DEHP displaying larger early deviations followed by a plateau (Fig.  7 A). The radius of gyration (Rg) remained essentially constant across simulations, indicating stable compactness of the complexes (Fig.  7 B). SASA changed slowly with small amplitudes, with the HMGCR system showing a higher surface area and slightly larger fluctuation, reflecting local solvent rearrangement in the binding pocket (Fig.  7 C). Protein–ligand hydrogen bonds fluctuated between 0 and 4, most often 1–2, providing sustained interaction constraints (Fig.  7 D). Residue-level RMSF was generally low (most sites < 4 Å), with higher flexibility confined to loop/terminal regions (Fig.  7 E, G). Collectively, these metrics indicate stable, compact, and persistent binding on the nano- to sub-microsecond scale, further reinforcing the mechanistic plausibility that DEHP engages key nodes of lipid metabolism and membrane dynamics in Ems. Fig. 1 Flow-chart of datasets analysis in this paper. Inputs: GEO training matrices ( GSE51981 / GSE6364 / GSE7305 ) and external validation ( GSE11691 / GSE23339 / GSE25628 ); DEHP targets from ChEMBL/Pharm Mapper/Swiss Target Prediction. Preprocessing: background correction, quantile normalization, probe‑to‑gene aggregation, ComBat batch adjustment; quality control (QC) by PCA. Disease‑side discovery: DEGs and WGCNA. Exposure–disease convergence: intersection, PPI, and enrichment context. Modeling: nested‑cross-validation (CV) training, SHAP interpretability, and independent validation. Structure/dynamics: Vina docking and 100‑ns GROMACS MD on representative complexes. Outputs: 17‑gene overlap, seven‑gene core signature, and translational hypotheses Flow-chart of datasets analysis in this paper. Inputs: GEO training matrices ( GSE51981 / GSE6364 / GSE7305 ) and external validation ( GSE11691 / GSE23339 / GSE25628 ); DEHP targets from ChEMBL/Pharm Mapper/Swiss Target Prediction. Preprocessing: background correction, quantile normalization, probe‑to‑gene aggregation, ComBat batch adjustment; quality control (QC) by PCA. Disease‑side discovery: DEGs and WGCNA. Exposure–disease convergence: intersection, PPI, and enrichment context. Modeling: nested‑cross-validation (CV) training, SHAP interpretability, and independent validation. Structure/dynamics: Vina docking and 100‑ns GROMACS MD on representative complexes. Outputs: 17‑gene overlap, seven‑gene core signature, and translational hypotheses Fig. 2 Systematic identification of DEHP targets. A Chemical structure of DEHP (PubChem). B Integrated in silico prediction from ChEMBL, PharmMapper, and SwissTargetPrediction; after de-duplication, the union comprises 1364 candidate protein targets Systematic identification of DEHP targets. A Chemical structure of DEHP (PubChem). B Integrated in silico prediction from ChEMBL, PharmMapper, and SwissTargetPrediction; after de-duplication, the union comprises 1364 candidate protein targets Fig. 3 Transcriptomic integration and module discovery in endometriosis (Ems). A PCA before batch correction (PC1 = 62.4%, PC2 = 11.9%). B PCA after batch correction (PC1 = 21.9%, PC2 = 8.8%). C Volcano plot of differentially expressed genes (|log2FC| ≥ 1, FDR < 0.05). D Heatmap of DEGs across samples (controls, n  = 47; cases, n  = 64). E WGCNA dendrogram with module assignment. F Module–trait correlations (MEbrown: r  = 0.58, P  = 3.6 × 10⁻¹¹; MEblue: r  = 0.57, P  = 6.5 × 10⁻¹¹). G Venn diagram showing the 229-gene intersection between DEGs and the MEblue module (2,537 genes) Transcriptomic integration and module discovery in endometriosis (Ems). A PCA before batch correction (PC1 = 62.4%, PC2 = 11.9%). B PCA after batch correction (PC1 = 21.9%, PC2 = 8.8%). C Volcano plot of differentially expressed genes (|log2FC| ≥ 1, FDR < 0.05). D Heatmap of DEGs across samples (controls, n  = 47; cases, n  = 64). E WGCNA dendrogram with module assignment. F Module–trait correlations (MEbrown: r  = 0.58, P  = 3.6 × 10⁻¹¹; MEblue: r  = 0.57, P  = 6.5 × 10⁻¹¹). G Venn diagram showing the 229-gene intersection between DEGs and the MEblue module (2,537 genes) Fig. 4 Convergence of DEHP targets with endometriosis genes and functional context. A Venn diagram intersecting 1,364 predicted DEHP targets with 229 disease genes (DEGs ∩ WGCNA MEblue), yielding 17 overlaps. B STRING PPI network (confidence ≥ 0.4), highlighting connectivity among lipid-metabolic, vesicular, and detoxification nodes. C Curated core subnetwork emphasizing HMGCR, ELOVL6, UGT8, SLC1A5, VAMP2, LYPLA1, EPHX1. D GO/KEGG enrichment for the 17 genes (dot size = gene count; color = adjusted P), with themes in bile secretion, serine/carboxylate ester hydrolase activity, lipase activity, and endocytic/phagocytic vesicle membranes, consistent with a lipid–vesicle–detox triad Convergence of DEHP targets with endometriosis genes and functional context. A Venn diagram intersecting 1,364 predicted DEHP targets with 229 disease genes (DEGs ∩ WGCNA MEblue), yielding 17 overlaps. B STRING PPI network (confidence ≥ 0.4), highlighting connectivity among lipid-metabolic, vesicular, and detoxification nodes. C Curated core subnetwork emphasizing HMGCR, ELOVL6, UGT8, SLC1A5, VAMP2, LYPLA1, EPHX1. D GO/KEGG enrichment for the 17 genes (dot size = gene count; color = adjusted P), with themes in bile secretion, serine/carboxylate ester hydrolase activity, lipase activity, and endocytic/phagocytic vesicle membranes, consistent with a lipid–vesicle–detox triad Fig. 5 Interpretable machine learning identifies a DEHP-related core signature in endometriosis. A Cross-cohort AUC heatmap across multiple algorithms (most AUCs > 0.75). B Single-gene ROC curves with top transferability for UGT8 (AUC 0.869) and EPHX1 (0.853), followed by VAMP2 (0.797), ELOVL6 (0.784), HMGCR (0.776), SLC1A5 (0.685). C Volcano plot with directionality (e.g., EPHX1 up, UGT8/ELOVL6/LYPLA1 down). D – E Global SHAP importance and a representative waterfall plot. F – M SHAP interaction plots illustrating ELOVL6–UGT8 synergy, HMGCR–LYPLA1 antagonism, and EPHX1–SLC1A5 context dependence—mechanistic patterns consistent with lipid reprogramming and nutrient influx Interpretable machine learning identifies a DEHP-related core signature in endometriosis. A Cross-cohort AUC heatmap across multiple algorithms (most AUCs > 0.75). B Single-gene ROC curves with top transferability for UGT8 (AUC 0.869) and EPHX1 (0.853), followed by VAMP2 (0.797), ELOVL6 (0.784), HMGCR (0.776), SLC1A5 (0.685). C Volcano plot with directionality (e.g., EPHX1 up, UGT8/ELOVL6/LYPLA1 down). D – E Global SHAP importance and a representative waterfall plot. F – M SHAP interaction plots illustrating ELOVL6–UGT8 synergy, HMGCR–LYPLA1 antagonism, and EPHX1–SLC1A5 context dependence—mechanistic patterns consistent with lipid reprogramming and nutrient influx Fig. 6 Molecular docking of DEHP with core proteins. A HMGCR. B EPHX1. C ELOVL6. D SLC1A5. E VAMP2. F LYPLA1. G UGT8. Yellow dashed lines indicate hydrogen bonds (~ 2.2–3.3 Å); additional hydrophobic and π–π contacts reflect pocket complementarity Molecular docking of DEHP with core proteins. A HMGCR. B EPHX1. C ELOVL6. D SLC1A5. E VAMP2. F LYPLA1. G UGT8. Yellow dashed lines indicate hydrogen bonds (~ 2.2–3.3 Å); additional hydrophobic and π–π contacts reflect pocket complementarity Fig. 7 Molecular-dynamics validation of DEHP–protein complexes. A RMSD. B Radius of gyration (Rg). C Solvent-accessible surface area (SASA). D Protein–ligand hydrogen-bond counts. (E) RMSF profile for ELOVL6. (F) RMSF profile for HMGCR. (G) RMSF profile for UGT8. Across metrics, all three complexes show stable, compact, and persistent binding over 100 ns Molecular-dynamics validation of DEHP–protein complexes. A RMSD. B Radius of gyration (Rg). C Solvent-accessible surface area (SASA). D Protein–ligand hydrogen-bond counts. (E) RMSF profile for ELOVL6. (F) RMSF profile for HMGCR. (G) RMSF profile for UGT8. Across metrics, all three complexes show stable, compact, and persistent binding over 100 ns

To define dataset eligibility and preprocessing to produce a harmonized training/validation matrix with minimized batch effects, we retrieved six endometriosis (Ems) transcriptomic datasets from Gene Expression Omnibus (GEO) (raw CEL/TXT files): (i)Training set: GSE51981 , GSE6364 , and GSE7305 (controls n  = 47; cases n  = 64); (ii)Independent validation set: GSE11691 , GSE23339 , and GSE25628 . Using R (v4.2.2), background correction and quantile normalization were performed with affy and limma12; probe IDs were mapped to gene symbols via platform annotations (multiple probes per gene aggregated by median). To mitigate platform heterogeneity, batch effects were adjusted with ComBat (sva), and principal component analysis (principal component analysis (PCA) via prcomp) was used to assess pre/post-correction structure. The training set samples were merged into a unified matrix for model development, while the validation set was processed separately using the same pipeline for external testing. Group sizes were inspected for balance and handled by weighting or stratification where appropriate. The complete analytical workflow is schematized in Fig.  1 . Dataset eligibility and selection (training cohorts: GSE51981 , GSE6364 , GSE7305 ). We screened GEO for human endometriosis case–control transcriptomes and pre-specified the following eligibility: (i) human tissue (ectopic endometriotic lesion and/or eutopic endometrium) from reproductive-age participants; (ii) no in-vitro drug/chemical stimulation before RNA profiling; (iii) raw data (CEL/TXT) and full sample annotations available; (iv) microarray or RNA-seq platforms with broad gene coverage (≥ 10k genes); (v) ≥ 10 total samples/study; and (vi) clear case–control labels. We excluded cell-line and in-vitro exposure studies, non-human datasets, and cohorts lacking raw files or phenotype labels. GSE51981 , GSE6364 , and GSE7305 met all criteria and were merged as the training set (controls n  = 47; cases n  = 64). Raw matrices underwent background correction, quantile normalization, probe-to-gene mapping (median across probes), and ComBat batch adjustment; PCA confirmed attenuation of platform-driven structure and improved case–control separation. The identical pipeline was applied to GSE11691 , GSE23339 , and GSE25628 , which we retained solely for external validation (no cross-study leakage). To assemble a high-recall DEHP target space via multi-algorithm prediction and to contextualize it with PPI and enrichment analyses, we identified differentially expressed genes (DEGs) between Ems and controls using limma (|log₂FC| ≥ 1; FDR < 0.05, Benjamini–Hochberg) [ 12 ]. Volcano and heatmaps were generated to confirm cross-sample coherence. A weighted gene co-expression network was constructed with WGCNA: soft-threshold power (β) was chosen by the scale-free topology criterion; adjacency and topological overlap matrices (TOM) were computed; dynamic tree cutting produced modules. Module–trait correlations (Pearson) were calculated between module eigengenes (MEs) and the Ems phenotype. The disease-associated module with optimal correlation, size, and biological interpretability was selected for integration with DEGs to obtain a high-confidence Ems candidate set. To evaluate the structure-level feasibility of DEHP engagement with core proteins, we performed molecular docking as an orthogonal validation of chemical–protein complementarity to support actionability. The DEHP structure (CID: 8343) was retrieved from PubChem.Putative targets were predicted using ChEMBL, PharmMapper, and SwissTargetPrediction; results were merged (union) and de-duplicated to yield unique candidate proteins. Protein–protein interactions were queried in Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) v11.5 (Homo sapiens; minimum confidence 0.4; disconnected nodes hidden for the main graph). clusterProfiler was used for GO/KEGG enrichment (p.adjust < 0.05), with visualization via ggplot2/enrichplot [ 13 , 14 ]. When experimental structures were unavailable downstream, AlphaFold Protein Structure Database entries were also consulted alongside Protein Data Bank (PDB) records [ 15 ]. We implemented a modeling panel in caret/mlr [ 16 ] that included logistic regression with L1/L2 penalties, lasso, ridge, naïve Bayes, random forest, support vector machine (SVM), gradient boosting/Extreme Gradient Boosting (XGBoost), and glmBoost [ 17 ]. All model development and parameter tuning were performed exclusively on the training set ( GSE51981 + GSE6364 + GSE7305 ). Performance was assessed by stratified 5-fold cross-validation (10 repeats) with nested CV for hyperparameter tuning to prevent leakage; class imbalance was addressed via class weights and repeated undersampling. The final model validation was conducted on the independent validation set ( GSE11691 + GSE23339 + GSE25628 ). The primary metric was AUC (with accuracy, sensitivity, specificity, and DeLong 95% CIs). Single-gene ROC/AUC estimation with 2,000 bootstrap replicates used pROC [ 18 ]. Feature contributions were quantified using SHAP (packages iml/shapley), and we generated SHAP importance rankings, global summary plots, representative waterfall plots, and pairwise dependence/interaction plots. The three-dimensional ligand (DEHP) was downloaded from PubChem; receptor structures were taken from the PDB or AlphaFold models when experimental structures were unavailable. Structures were prepared in AutoDock Tools 1.5.7 (water removal, hydrogen addition, Gasteiger charges, rotatable bonds). Docking was performed with AutoDock Vina 1.2.0 using grid boxes centered on the active pocket (covering key residues); default exhaustiveness = 8 was increased (up to 32) when needed for convergence. Lowest binding-energy poses were retained; binding energies < − 5.0 kcal/mol were considered indicative of strong interactions. PyMOL 2.5.0 was used for 3D visualization and LigPlot + 2.2 for interaction schematics (hydrogen bonds, hydrophobic and π–π/π–cation contacts). We selected UGT8, ELOVL6, and HMGCR based on (i) top docking scores and pose convergence, (ii) centrality within the lipid‑metabolism/vesicular axis, and (iii) complementary biology (sphingolipid synthesis, acyl‑chain elongation, cholesterol biosynthesis). Additional complexes will be pursued in follow‑up to balance cost with coverage. Definitions. PME = particle‑mesh Ewald for long‑range electrostatics; RMSD = root‑mean‑square deviation; RMSF = root‑mean‑square fluctuation (Cα unless specified). Representative complexes—UGT8–DEHP, ELOVL6–DEHP, and HMGCR–DEHP—underwent all-atom MD with GROMACS 2022 and the AMBER99SB-ILDN force field [ 19 ]; ligand parameters were generated using General Amber Force Field 2 (GAFF2)/Antechamber (Austin Model 1–Bond Charge Corrections (AM1-BCC) charges). Systems were solvated with Transferable Intermolecular Potential with 3 Points (TIP3P) water and neutralized with 0.15 M NaCl. After energy minimization, equilibration comprised 1 ns constant number, volume, and temperature (NVT) and 1 ns constant number, pressure, and temperature (NPT); production runs were 100 ns with a 2 fs timestep at 300 K (V-rescale thermostat) and 1 bar (Parrinello–Rahman barostat). Long-range electrostatics used PME; covalent bonds were constrained with Linear Constraint Solver (LINCS). Trajectories were analyzed for RMSD, radius of gyration (Rg), solvent-accessible surface area (SASA), hydrogen bonds, and RMSF (Cα-based), sampled every 10 ps. Unless otherwise stated, we applied the two‑sided Wilcoxon rank‑sum test and adjusted p‑values by Benjamini–Hochberg FDR. All analyses were conducted in R 4.2. Group comparisons used the Wilcoxon rank-sum test; p values were adjusted via Benjamini–Hochberg FDR (default threshold 0.05). Random seeds were fixed ( set.seed  = 202405). Software/database versions and access dates were recorded. All data were obtained from public repositories; ethics approval was not applicable. Analysis scripts and parameter files are available upon reasonable request.

This study constructs a closed evidentiary loop that proceeds from chemical exposure to protein targets, through disease networks and interpretable machine learning, and onward to structural docking and molecular dynamics, using DEHP as an index phthalate to illuminate how membrane-lipid homeostasis, vesicular trafficking, and detoxification metabolism may converge in Ems. By integrating multi-cohort transcriptomes and WGCNA, we derived a high-confidence disease gene set; intersecting it with DEHP target predictions yielded 17 overlap genes, from which SHAP prioritized a seven-gene core signature (ELOVL6, LYPLA1, UGT8, SLC1A5, HMGCR, EPHX1, VAMP2). Molecular docking and 100-ns simulations supported the structural feasibility and temporal stability of key complexes. Altogether, the results align with established hallmarks of metabolic remodeling and membrane-lipid abnormalities in Ems and offer actionable molecular anchors for translating environmental exposure into microenvironmental change [ 5 , 20 , 21 ]. From a pathophysiological standpoint, Ems exhibits broad perturbations in energy and lipid programs, increasingly recognized as central to disease persistence and immune remodeling. Recent reviews, lipidomics, and genetic epidemiology converge on a “lipid metabolism–inflammation–cellular plasticity” axis. Our enrichment terms (bile secretion, lipase/ester-hydrolase activity, endocytic/phagocytic vesicle membranes) and the PPI backbone (lipid synthesis—vesicle dynamics—genome surveillance) are consistent with these trends. Notably, a 2025 multi-omics Mendelian randomization analysis reported putative causal links between specific lipid metabolites and Ems risk, providing orthogonal support for our lipid-centric hypothesis [ 5 , 22 – 24 ]. Functionally, the seven-gene signature maps to coherent nodes: ELOVL6 extends long-chain fatty acids, shaping membrane composition; HMGCR controls cholesterol biosynthesis; UGT8 couples to sphingolipid synthesis; LYPLA1 depalmitoylates proteins, regulating membrane localization; SLC1A5 governs glutamine uptake and metabolic rewiring; VAMP2 drives vesicle fusion; and EPHX1 detoxifies epoxide intermediates. This “lipid–vesicle–detox/nutrient” triad mirrors well-known DEHP/MEHP effects on PPAR and LXR/SREBP axes, unifying our differential directions (suppressed lipid remodeling; activated detoxification) with plausible upstream endocrine-disrupting mechanisms [ 5 , 9 , 21 , 25 – 27 ]. The appearance of vesicle-related terms and VAMP2 further implicates extracellular vesicles (EVs) as couriers of exposure-to-phenotype signals. Contemporary reviews and primary studies document EV involvement in immune crosstalk, fibrogenesis, neuroangiogenesis, and ectopic implantation in Ems, with diagnostic promise from menstrual-blood-derived or peripheral EV cargos. Rigorous validation should follow Minimal Information for Studies of Extracellular Vesicles (MISEV 2023) guidelines and leverage multi-omics profiling of EV lipids/proteins/miRNAs to test our signature in clinically accessible biospecimens [ 28 ]. Regarding exposure plausibility, population biomonitoring indicates that although several DEHP metabolites (e.g., MEHHP, MEOHP) have declined over the last decade, exposures persist with substantial geographic and demographic variation; urinary biomarkers remain the preferred readout for short half-life phthalates. Future prospective designs should collect repeated urine samples to control within-person variability and incorporate individualized mixture profiles (phthalate substitutes, BPA, etc.) and time-window effects to resolve the “exposure–signature–phenotype” triangle [ 29 , 30 ]. Translationally, the signature suggests two directions. First, as a diagnostic/stratification panel, it can be validated across lesion tissue, eutopic endometrium, peritoneal fluid, and menstrual-blood EVs, and linked to pain, recurrence, and infertility outcomes. Second, as therapeutic entry points, both the HMGCR axis and glutamine transport (SLC1A5) have pharmacological precedents elsewhere, yet their efficacy and safety in Ems remain unsettled. Preclinical and early clinical signals around statins are mixed, warranting cautious, fertility-aware trials; SLC1A5 inhibitors (e.g., V-9302 and successors) show activity and combinational potential in immunometabolic oncology but face specificity/toxicity hurdles that must be addressed before gynecologic application [ 31 – 35 ]. Clinical implications: diagnostics and therapeutics. The seven-gene signature (ELOVL6, LYPLA1, UGT8, SLC1A5, HMGCR, EPHX1, VAMP2) provides a translation-ready scaffold. Diagnostics: a minimal qPCR/reverse-transcription droplet digital PCR (RT‑droplet digital PCR (ddPCR)) panel in lesion/eutopic tissue, peritoneal fluid, or menstrual‑blood–derived extracellular vesicles can be optimized for rule‑in/rule‑out performance and linked to pain, recurrence, and fertility outcomes [ 36 – 38 ]. Feature directions (e.g., EPHX1 up, UGT8/ELOVL6/LYPLA1 down) offer immediate interpretability and could be integrated with clinical variables to yield nomograms. Therapeutics: the HMGCR axis suggests hypothesis‑driven evaluation of statin‑based strategies with fertility‑aware monitoring [ 39 – 41 ]; SLC1A5‑mediated glutamine uptake nominates immunometabolic interventions to be tested in preclinical endometriosis models before gynecologic translation [ 42 , 43 ]. Given DEHP’s network engagement, combination strategies (lipid‑modulating agents with hormonal or anti‑inflammatory regimens) merit systematic, exposure‑stratified exploration. Next steps: under MISEV‑compliant workflows, we recommend paired exposure biomarker collection (urinary DEHP metabolites45), longitudinal sampling, and perturbational experiments (CRISPR/siRNA of HMGCR/ELOVL6/UGT8/SLC1A5) to test causality along the lipid–vesicle–detox corridor. Limitations include: (i) clinical heterogeneity and sample imbalance in public transcriptomes; (ii) reliance on in silico DEHP targeting, despite multi-database integration and convergence with disease evidence; (iii) docking/MD constraints (force fields, pocket definitions, finite timescales) that cannot fully capture slow conformational transitions and membrane effects—even with 100-ns stability; and (iv) absence of paired exposure data precluding causal-chain testing. Next steps should, under MISEV2023 compliance, couple primary cells, organoids, and animal models with longitudinal clinical cohorts: CRISPR/perturbational tests of core genes; concurrent exposure plus lipidome/EV profiling; and synergy/antagonism trials with hormonal or immunomodulatory regimens [ 44 , 45 ].

In summary, we propose and substantiate a DEHP–Ems model scaffolded by lipid metabolism, vesicular transport, and detoxification, nominating ELOVL6, LYPLA1, UGT8, SLC1A5, HMGCR, EPHX1, and VAMP2 as actionable hubs for diagnostics and intervention. This framework offers testable hypotheses and tractable sampling/therapeutic routes for mechanistic dissection and clinical translation, while exemplifying a generalizable pipeline for connecting environmental exposures to chronic gynecologic disease.

Endometriosis is a chronic, estrogen-dependent condition characterized by diagnostic delay and frequent recurrence. Multiple human studies—including meta-analyses—now link di(2-ethylhexyl) phthalate (DEHP) and its metabolites to endometriosis [ 1 – 3 ]. A 2019 meta-analysis of 30 epidemiologic studies reported an elevated risk for phthalates overall and a significant association for DEHP (odds ratio ≈ 1.42, 95% CI 1.19–1.70) [ 4 ], and a 2021 PRISMA meta-analysis found higher urinary MEOHP/MEHHP and higher blood DEHP/MEHP among cases versus controls [ 5 ]. Case–control investigations further corroborate higher DEHP/MEHP burdens in plasma, urine, and even peritoneal fluid of affected women [ 6 – 8 ]. Mechanistically, independent in vivo and in vitro work shows that DEHP perturbs lipid and endocrine programs relevant to endometriosis by activating LXR/SREBP-1c and PPAR-α/γ signaling, amplifying inflammatory pathways, and broadly reshaping lipid metabolism [ 9 – 11 ]. Despite this convergence, a cohesive exposure-to-biology chain that pinpoints actionable molecular nodes has remained incomplete. Here, we develop an integrated framework centered on DEHP, spanning statistical, network, structural, and dynamical analyses. We aggregate predicted targets (ChEMBL, PharmMapper, SwissTargetPrediction), integrate multi-cohort endometriosis transcriptomes with differential expression and weighted gene co-expression network analysis, contextualize the exposure–disease overlap through PPI and GO/KEGG, and use interpretable machine learning to prioritize features. Molecular docking (AutoDock Vina 1.2.0 with AlphaFold structures) and 100-ns MD simulations (GROMACS 2022) provide orthogonal structural/dynamical validation. The result is a compact, exposure-anchored seven-gene signature (ELOVL6, LYPLA1, UGT8, SLC1A5, HMGCR, EPHX1, VAMP2) that links membrane-lipid homeostasis, vesicular transport, and detoxification, supporting translation toward diagnostics and therapeutics (e.g., the HMGCR axis; SLC1A5-mediated glutamine uptake).