EMs transcriptomic microarray data were obtained from the NCBI’s Gene Expression Omnibus (GEO; https://ncbi.nlm.mih.gov/geo/), and data on the clinical features of patients examined in the microarray analyses were obtained from EMBL-EBI’s ArrayExpress (https://www.ebi.ac.uk/arrayexpress). The expression of genes of interest was analyzed in endometriotic tissue using the ENDOMET Turku Endometriosis Database (9). A total of 198 ectopic EMs lesion samples (including 6 biological replicates) from patients with EMs were analyzed from the GSE141549 dataset, which contains the expression matrix and clinical information of 198 EMs lesions. Additionally, three other independent transcriptomic profiles (GSE25628, E-MTAB-694, and GSE23339) were analyzed for further validation.

Expression values from both datasets were log2-transformed before cross-platform normalization. The ComBat function from the SVA package (in the R environment, version 3.6.2) (10) was used for the meta-analysis and data cleaning to remove batch effects. The principal component analysis (PCA) was used to evaluate whether the known batch effects were removed. Samples from the same dataset were considered to have no obvious batch effect if they did not cluster together.

Normalized expression values of 198 EMs lesions from the GSE141549 dataset were used to identify molecular subtypes by applying the consensus clustering method, which was implemented using the ConsensusClusterPlus package (in the R environment, version 1.58.0) (11). Consensus clustering was implemented with the following settings: maximum cluster number (maxK) = 10, number of repeats (reps) = 10,000, proportion of items to sample (pItem) = 0.8, proportion of features to sample (pFeature) = 1, cluster algorithm (clusterAlg) = “km” (K-means), and distance = “Euclidean”. The optimal cluster number was determined based on the consensus matrix and the cluster consensus score. The consensus score for k = 2 was larger than that of other clusters. Finally, the EMs lesions were clustered into two molecular subtypes: stroma-enriched subtype (S1) and immune subtype (S2).

The biological function of each subgroup was investigated by identifying clusters of co-expressed characteristic genes using the WGCNA package (in the R environment, version 4.1.3) (12), which was also used to construct the co-expression network. The process also involved calculating the Pearson’s correlation matrices among all of the genes. Using the pickSoftThreshold function, 8 was determined to be the appropriate soft threshold to create the weighted adjacency matrix. Subsequently, hierarchical clustering and the dynamic tree cut function were used to assign genes to modules. The blockwiseModules function was applied using the following parameters: maxBlockSize = 30,000, minModuleSize = 30, mergeCutHeight = 0.25, reassignThreshold = 0.

Pearson’s correlation coefficients with P values estimated using the corPvalueStudent function were used to assess the correlations between molecular subtypes and module eigengenes.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Gene Ontology (GO) analysis were performed using the clusterProfiler package (in the R environment, version 4.2.2) to explore the module-related pathways.

Cell type composition scores were estimated for the S1 and S2 subtypes using the xCellAnalysis function from the xCell package (in the R environment, version 1.1.0) (13). Correlations between module eigengenes and computed cell type scores were visualized for all cell types that were scored in more than 25% of the samples, which allowed us to infer the representative cell type of each different module.

The infiltrating cells in lesions were also estimated using CIBERSORT (14), which yielded similar results to xCell.

The EaSIeR package (in the R environment, version 1.0.0) was used to predict the outcome of immune therapy for the different subtypes (15). EaSIeR provided the estimation of immune responses for each lesion based on RNA-seq data combined with prior knowledge, and lesions with a higher relative score had a stronger positive correlation with the immunotherapy response.

To explore the subtype-related biological processes, the GSEA (in GSEA software version 4.2.3) was conducted for each subtype (16). The KEGG pathway enrichment results were considered statistically significant based on the net enrichment score (NES), gene ratio, and P value. Gene sets with a |NES| > 1 and a NOM p < 0.05 were considered to be significantly enriched.

The NetworkAnalyst (https://www.networkanalyst.ca/) online tool was used to construct and visualize the PPI network of gene models related to both subtypes (17). Genes with a connectivity degree of greater than 10 were identified as hub genes.

Differentially expressed genes (DEGs) between the two different subtypes were identified using the limma package (in the R environment, version 3.50.1). Statistically significant DEGs were identified with the P value cutoff of <0.01 and a fold change of ≥2 or ≤−2.

Subtype-related predictive genes were assessed in the GSE1141549 dataset using the glmnet package (in the R environment, version 4.1-4) based on the least absolute shrinkage and selection operator (LASSO) method and the randomForest package (in the R environment, version 4.1.3) based on Breiman’s random forest algorithm. The identified genes were used to build a linear regression model to examine the associations between the predictor genes and the molecular subtypes. The reproducibility was validated using the enhanced bootstrap method (100 bootstrap rounds). Subsequently, the model was tested on the three other datasets (GSE23339, GSE25628, and E-MTAB-694). C-statistics and Brier scores were calculated to assess the prediction performance. The C-index values were between 0.5 and 1.0, with 0.5 and 1.0 representing random opportunity and an excellent ability of the model to predict the subtype, respectively. The Brier score values were between 0 and 1, with 1 representing an entirely inaccurate forecast. A lower Brier score of a set of predictions indicated that the predictions were better calibrated, with the best possible Brier score of 0 representing total accuracy.

All tissue samples were obtained with informed consent from patients and approved by the ethics service under the research ethics committee numbers (kyy-2019-105). Premenopausal women diagnosed with EMs were included from the Obstetrics and Gynecology Hospital of Fudan University between January 2020 and September 2022. Samples were obtained from patients undergoing surgical resection of EMs lesions (conservative or radical surgery) and confirmed by pathological examination. All participants underwent hormone therapy before surgery. Detailed clinical information about previous EMs surgery and hormone therapy was collected for all recruited participants. Patients who underwent hormone therapy for less than 1 month were considered unexposed (18). The criteria for patients considered as failed/intolerant to hormone therapy were as follows: 1) patients who underwent hormone therapy for more than 1 month without any symptom relief and 2) patients who underwent hormone therapy for more than 1 month and who terminated treatment due to intolerable side effects.

Patients who met the following criteria were excluded: 1) patients with a history of any inflammatory condition, autoimmune disease, or malignant tumor; 2) patients who stopped treatment due to drug allergy or severe side effects; 3) patients who underwent hormone therapy for other conditions.

All tissues used to construct TMAs were evaluated by an experienced gynecological pathologist in advance of the presence of EMs lesions. When typical endometrial glands were found in the lesions, the pathological diagnosis of EMs was made. CD10 immunohistochemical staining of lesions was performed when the pathologists failed to find typical endometrial glands (e.g., due to heat injury). In total, 97 formalin-fixed, paraffin-embedded EMs lesions (22 intestinal lesions; 21 deep rectovaginal lesions; 11 other deep lesions, including ureteral or vaginal lesions; 36 ovarian lesions; and 7 peritoneal lesions) were used to construct the TMAs. The intestinal lesions were arrayed in two different paraffin blocks, whereas the deep rectovaginal lesions were arrayed in three different paraffin blocks. Ureteral lesions and vaginal lesions were arrayed in one paraffin block, ovarian lesions were arrayed in three different paraffin blocks, and peritoneal lesions were arrayed in one paraffin block. Tissue cylinders with a diameter of 5 mm were used to extract tissue from the targeted area confirmed by the pathologist of each donor tissue block, which were then deposited into recipient blocks. After the array blocks were constructed, they were cut into multiple 4-μm sections until all 97 tissue samples were represented on a single section. Each section was placed on a microscopic slide and histologically investigated after hematoxylin and eosin (H&E) staining to determine the adequacy of the arrayed tissues. These sections were separately placed on charged polylysine-coated slides for immunohistochemistry (IHC).

The paraffin-coated microarray sections were placed on a heating block at 60°C for 2 hours and continuously washed with xylene to remove the paraffin. The slides were then rehydrated in varying concentrations of alcohol. Subsequently, 3% hydrogen peroxide was applied for 30 min at room temperature to block endogenous peroxidase activity. Then, the slides were incubated in the retrieval buffer (sodium citrate solution, pH = 6.0) and heated in the microwave to restore the antigen. The slides were pre-incubated with serum for 1 hour to reduce non-specific background. The slides were incubated overnight with primary antibodies diluted in phosphate-buffered saline (PBS) at 4°C. Four monoclonal antibodies were used to detecting PR (ab32085, 1:150), PI16 (ab127014, 1:500), FHL1 (ab133661, 1:100), and SORBS1 (ab224129, 1:200). The next day, the slides were washed three times with PBS (1×) for 5 min each. Subsequently, the slides were incubated with the secondary antibody (abs20040, 1:500) for 2 h at room temperature. The slides were developed in diaminobenzidine solution and stained with hematoxylin. Finally, images of the representative fields in each case were collected using the Stream Software, and expression was quantified using Image J software. The FHL1 and SORBS1 staining results were assessed as negative or positive by a pathologist who was unfamiliar with the clinical pathological data. Correlations between FHL1 or SORBS1 expression and the effects of hormone therapy were estimated using the chi-square test. Correlations between the expression of other genes and the effects of hormone therapy were estimated using the t-test, with or without Welch’s correction. All analyses were performed using GraphPad Prism v.9.0 software. P values of <0.05 were considered statistically significant.

Microarray EMs transcriptomic data and matched clinical data were obtained from NCBI’s GEO (https://ncbi.nlm.mih.gov/geo/) and EMBL-EBI’s ArrayExpress (https://www.ebi.ac.uk/arrayexpress). In total, the data of 198 EMs samples (91 deep lesions, 79 peritoneal lesions, and 28 ovarian lesions) from the GSE1412549 dataset were used as the training dataset. The other datasets, including GSE25628 (7 deep lesions), E-MTAB-694 (18 peritoneal lesions), and GSE23339 (10 ovarian lesions), were used for further external validation. A detailed flowchart explaining the process is illustrated in Figure 1 . To remove any batch effects caused by the data being from different platforms and different batches, the ComBat method was applied between the datasets. The PCA was also conducted on these datasets to establish the relationships among the four validation datasets. As shown in Figure 2A , samples from the four independent datasets formed different clusters before removing the batch effects; however, they clustered together after batch effect removal ( Figure 2B ). This clustering indicated that cross-platform normalization was successful in removing the batch effects.

To understand the heterogeneity of EMs and identify potential subtypes, the GEO dataset GSE141549 was downloaded along with clinical data ( Supplementary Table 1 ) to screen EMs gene expression. EMs subtypes were identified from the high-dimensional dataset by applying the unsupervised clustering algorithm to classify the ectopic lesions based on the heterogeneity of the gene expression profiles. Normalized expression values of the 198 EMs lesions from GSE141549 were used to identify the different molecular subtypes. As shown in Figure 3 , the EMs lesions were finally clustered into two molecular subtypes: S1 and S2. We observed that the S1 and S2 subtypes contained 109 and 88 lesions, respectively. The consensus matrix suggested a high degree of intra-group homogeneity and distinct heterogeneity between the two subtypes ( Figure 3A ). The cluster consensus value for each subgroup suggested that 2 was the optimal number of clusters ( Figure 3B ). The PCA analysis suggested that the two subtypes were separated at the transcriptional level ( Figure 3C ).

The detailed clinical characteristics of the 198 EMs patients along with clinical data from the GSE141549 dataset are shown in the additional file ( Supplementary Table 1 , Supplementary Figure 1 ). Given that the clinical characteristics might have an impact on the gene expression profile of lesions, the clinical characteristics of the two subtypes were investigated by examining the age, cycle, medication, revised American Fertility Society(rAFS) stage, and lesion locations of the EMs cases from the GSE141549 dataset. The proportion of lesion distribution in S1 was different from that in S2 ( Supplementary Table 1 ). However, there were no significant differences between S1 and S2 with regard to age, menstrual cycle, and rAFS stage ( Supplementary Table 1 , Supplementary Figure 2 ). These results indicated that the transcriptome classification might represent certain intrinsic biological characteristics.

The WGCNA was performed to identify clusters of co-expressed genes that were characteristic of the biological function of each subgroup. This process revealed 14 modules of highly co-expressed genes ( Supplementary Figure 3 ). After calculating the correlations between each module and clinical traits, it was evident that several gene modules were strongly positively correlated with a subtype (MEturquoise for S1, MEblue for S2) ( Figure 4A ). The functional enrichment analysis based on the KEGG database indicated that genes in MEturquoise were mainly enriched in pathways related to fibrosis, whereas genes in Meblue were mainly enriched in pathways related to the immune response and inflammation ( Figure 4B ). The GO analysis showed similar results ( Supplementary Figure 4 ).

The GSEA was performed for functional enrichment of S1 and S2 ( Supplementary Table 2 ). The KEGG enrichment terms revealed that S1 was associated with fibrosis-related terms, including vascular smooth muscle contraction and dilated cardiomyopathy, whereas S2 was associated with immune response-related terms, including primary immunodeficiency and allograft rejection. Notably, steroid biosynthesis pathways were also enriched in S2 ( Figure 4C ).

Next, hub genes were identified from the representative genes of each subtype, and PPI networks were developed for gene co-expression modules associated with molecular subgroups; the gene names of the top 30 nodes with the highest degree of connectivity are labeled in Figure 4D . The hub genes identified for S1 included PTPN11, PPKCA, CAV1, PPP1CB, PAP1A, CD44, CCND2, and PI16, which might have the largest impact on fibroblast activation and integrin signal transduction in milieu. The hub genes identified for S2 included LYN, STAT1, PTPN6, FGR, CXCR4, MYD88, CEBPA, and HCK, which have an important function in regulating the innate and adaptive immune responses, responses to growth factors and cytokines, and migration of immune cells.

We postulated that the observed gene co-expression patterns might also be indicative of cell type features associated with EMs. The differences in immune-related signatures between the EMs lesions of S1 and S2 were further explored. In this study, the xCell algorithm was used to analyze gene expression data (GSE141549) and calculate the module eigengene and immune-stroma scores of ectopic endometrial samples between S1 and S2. Strong positive correlations were observed between MEblue and obvious immune cell infiltration in S2, including that of macrophages and T lymphocytes ( Figure 5A ). The results revealed that the ectopic EMs lesions of S1 had significantly lower immune scores. In contrast, the stroma scores for ectopic EMs lesions of S1 were significantly higher than those of S2 (P < 0.0001). There was a significant difference between the microenvironment scores of the two subtypes (P > 0.05, Figure 5B ). To further explore the association between the S2 subtype and the response to immunotherapy, we used the EaSIeR method to predict the immunotherapy effect of different lesions according to the gene expression profile. The results suggested that the immune signature score of the S2 subtype was higher than that of the S1 subtype, which might represent a better response to immunotherapy ( Figure 5C ). Similar results were obtained using the CIBERSORT algorithm ( Figure 5D ). Compared with the neutrophil enrichment in the S1 subtype, higher infiltrating M0 macrophage, M2 macrophage, Treg cell, and activated dendritic cell ratios were observed in the S2 subtype of the endometrial ectopic lesions, indicating the characteristic of immune tolerance.

The gene labels identifying the two subtypes were explored by performing differential analysis, which identified 159 DEGs. Of these, 120 DEGs were significantly upregulated in S1, whereas 39 DEGs were downregulated ( Figure 6A , Supplementary Table 3 ). Subtype-related genes were identified using LASSO and the random forest algorithm based on the computed LASSO coefficient and the mean decrease in the Gini coefficient ( Figures 6B–D ). FHL1 and SORBS1 were selected to construct a predictive model using LASSO regression. These two genes were also differentially expressed between ectopic lesions and endometrium samples from patients ( Supplementary Table 4 ). Internal validation computed by 100 rounds of the bootstrap method indicated excellent discrimination (area under the receiver operating characteristic curve [AUC]: 0.97) and calibration (0.068) ( Figures 6E, F ).

For external validation, three independent EMs datasets (GSE25628, GSE23339, E-MTAB-694) were analyzed to verify the predictive performance of this model. After removing batch effects among the datasets, consensus clustering of endometriotic lesion samples from the training dataset and the three independent validation datasets is shown in Figure 7A . The final results revealed that the AUC was 0.86 and the Brier score was 0.16 ( Figures 7B, C ).

Immune infiltration estimates for the validation datasets determined by xCell yielded similar results to those of the training set ( Figures 7D–F ).

Hormone therapies are usually employed to treat EMs. However, our understanding of the molecular EMs subtypes and their associations with the clinical response to hormones remains incomplete. Furthermore, the significance of the EMs subtypes and their biomarkers requires validation across cohorts of well-annotated clinical samples and clinical features.

The unsupervised classification of the tissue transcriptome and gene signatures in EMs patients was validated by recruiting a validation cohort of 83 EMs patients. All patients participated in this study anonymously because of privacy and security concerns. The detailed baseline demographic information of the cohort is shown in Table 1 . A total of 97 ectopic EMs samples were collected during surgery. On the basis of the above, we measured the protein expression of the subtype markers and functional molecules identified in Figures 6A–D using IHC and distinguished hormone-sensitive and hormone-resistant populations. In brief, with the help of IHC analysis of the TMAs, FHL1 and SOBRS1-positive and high PI16 expression patients were classified into the S1 subtype, while FHL1 and SORBS1-negative and low PI16 expression patients were classified into the S2 subtype.

The subtype-specific gene signatures were identified in all patients. The associations between established EMs subtypes and clinical characteristics are reported in Table 1 . Overall, the subtypes were significantly associated with the lesion locations (P < 0.001). Next, the correlation between FHL1 positivity and the response to hormone therapy was examined for ectopic EMs lesions. A higher number of hormone-resistant EMs patients were found in FHL1-positive subgroup (S2) ( Figure 8 , top), which is consistent with a previous report showing that chronic inflammation might induce a progesterone-resistant state. However, we did not observe any correlation between SORBS1 expression and the response to hormone therapy. As recent research has reported in human tissues, a new universal fibroblast transcriptional subtype was identified across tissues. The PI16-positive fibroblast serves as a reservoir that can yield specialized fibroblasts across a broad range of steady-state tissues and activated fibroblasts under pathological conditions. In this study, IHC was applied to verify the expression of the universal fibroblast marker PI16, the same as PR, which was significantly decreased in the ectopic EMs tissues from resistant patients ( Figures 8 , middle and bottom). Overall, the highest numbers of endometrial stromal cells and fibroblasts were observed in the FHL1-positive subtype, indicating a subtype-specific response, which was likely involved in the distinct biological behavior of the different EMs subtypes.