RNA-seq data for EM were downloaded from the National Center for Biotechnology Information Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) (15, 16). Data were collected and analyzed using R 3.6.3 software. The GEOquery package was used to download the Matrix file from the GEO database (17).

The TRGs were obtained from TelNet (http://www.cancertelsys.org/telnet/) (18), a database of genes involved in telomere maintenance, including a systematic assessment of the telomere maintenance machinery (TMM) signature and its inclusion in TMM pathway models or genome-wide mutational analyses of the cancer genome.

The R package sva was used to remove batch effects from different datasets (19). Inter-sample correction effects were calculated using principal component analysis (PCA) clustering. The R package limma (20) was used to screening for differentially expressed genes (DEGs) and telomerase-related genes. The package ggplot2 was used to plot volcanoes of DEGs to visualize differential expression patterns. Genes with adjusted P (Padj) <0.05 and |log₂FC| > 1 were considered statistically significant.

Two machine learning methods, random forest model (RFM) and support vector machine (SVM)-recursive feature elimination (RFE), were used to screen for TRGs in EM. RFM is a classifier consisting of multiple decision trees, which ranks the importance of genes in regulatory relations. The SVM-RFE algorithm identifies the optimal variables by eliminating the feature vectors generated by the SVM (21). Receiver operating characteristic (ROC) curves and area under the curve (AUC) were used to evaluate the models' screening powers, resulting in a final set of EM telomere-associated signature genes.

A nomogram model was developed using the rms and rmda software packages to assess the predictive power of EM signature genes. Model accuracy was assessed using calibration curves, the Hosmer–Lemeshow test, clinical impact curve (CIC), and decision curve analysis (DCA) (22).

Unsupervised cluster analysis was applied to identify the different telomere gene subtypes and to classify samples for further analysis (23). The consensus clustering algorithm determined the number of clusters and their stability. The R package ConsensuClusterPlus was used to perform these steps, with 1,000 replications to ensure classification stability (24).

C57BL/6 mice (specific pathogen-free, female) were obtained from the Department of Laboratory Animals, Capital Medical University. The animal experiment protocol used herein was approved by the Animal Ethical Use Committee of Capital Medical University (Ethical Number: AEEI-2021-219). Five C57BL/6 mice were housed in a well-controlled, pathogen-free environment in a barrier unit with a regulated light/dark cycle (12/12 h, 23–25°C). One mouse was used as the endometrial tissue donor after euthanization. The donor uterus was evenly divided into 8 endometrial fragments of 2 mm diameter using a disposable biopsy punch (2 mm) (Integra Miltex, Shanghai, China). These were injected into the peritoneal cavities of anesthetized mice (~100 mg tissue/0.5 mL PBS per mouse) via a 2 ml syringe. All mice survived the duration of the experiment. There were no significant between-mouse growth rate differences. Mice were euthanized 4 weeks later and both EM and uterine tissues were collected for immunohistochemical analysis.

Mouse EM and uterine specimens were fixed in 4% neutral paraformaldehyde solution, followed by paraffin embedding. Serial sections (5 μm) were then dewaxed. The fixed endometrial and EM lesion sections were incubated overnight at 4°C with anti-microtubule-associated protein (MAP) 7 (1:400, Servicebio, Wuhan, China), anti-grainyhead-like 2 (GRHL2) (1:100, Affinity, Changzhou, China), anti-retinoic acid receptor response protein 2 (RARRES2) (1:50, ABclonal, Wuhan, China), anti-ubiquitin carboxyl-terminal hydrolase-L1 (UCHL1) (1:1,000, Servicebio), anti-estrogen receptor alpha (ESR1) (1:1,000, Servicebio), and anti-protein reversionless 3-like (REV3L) (1:100, Affinity). Horseradish peroxidase-labeled goat anti-mouse IgG (immunoglobulin) (1:200, Servicebio) was then added and incubated for 50 min at room temperature. Sections were washed with distilled water, counterstained with hematoxylin, dehydrated, and mounted for microscopy. Colored areas were quantified using SlideViewer and Image Pro Plus software. The whole tissue sections were scored for staining intensity and percentage. The scoring scale was: 0 (no staining), 1 (light brown staining), 2 (brown staining), and 3 (dark brown staining). The percentage of positive cells was graded into one of 4 levels: 1 (<5%), 2 (5–30%), 3 (31–60%), 4 (61–100%). Immunohistochemistry staining score was calculated as follows: intensity score × percentage score.

An ordered list of genes was generated based on correlations between all genes and telomerase-related gene expressions using Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA). The GO knowledgebase is the world's largest gene function information source (http://geneontology.org/) (25–27). KEGG is a collection of databases dealing with genomes, biological pathways, diseases, drugs, and chemical substances (www.kegg.jp/kegg/kegg1.html) (28–30). DEGs were defined by an absolute fold change >1.5 and Padj <0.05.

The GSEA computational method allows determination of overrepresented classes in large sets of genes or proteins that may be significantly association with disease phenotypes (31). The predefined gene set is from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) (32). Herein, an ordered list of genes was generated based on the correlation between all genes and telomerase-related gene expressions using GSEA. Enriched pathways were determined based on P-values and normalized enrichment scores.

The Alizadeh Lab CIBERSORT analytical tool was developed by Newman et al. to estimate the abundances of member cell types in a mixed cell population using gene expression data (33). It was used herein to assess the relative proportions of the 24 immune infiltrating cells in the telomerase and gene clusters of EM samples (33, 34). Correlation analyses were then performed between these 24 immune cells and critical genes.

To identify potential EM therapeutic agents, DEGs were uploaded to the connective map (CMAP) database (https://clue.io/) (35). Enrichment analysis was used to screen for relevant drugs with therapeutic potential (36). A negative enrichment value generally indicates that the drug is more likely to treat the disease, and a larger absolute value means that it is more disease-specific. Based on this, we screened for drugs with enrichments of <−0.5.

The study flowchart is shown in Figure 1. Using key words “endometriosis, Homo sapiens,” the GEO database expression profile was searched and the following datasets were selected for inclusion:

All three are mRNA-seq data GEO datasets, with large sample sizes and containing human EM and normal endometrial control samples. The SVA algorithm was used to batch-correct the three datasets and combine them into a single dataset of 35 normal and 28 EM samples. PCA was performed before and after removing the batch effect removal (Figures 2A, B).

We screened for DEGs between the disease and control groups, and created a volcano plot of the results (Figure 2C). A total of 2,093 TRGs were obtained from TelNet, including 165 validated genes, 923 predicted genes, and 1,005 screened genes. Further screening of genes differing between the disease and control groups was performed based on TRGs. The top 15 genes with the most significant differences are shown in the heatmap (Figure 2D). The positions of these DEGs were then labeled on the chromosomes. The top three genes with the most significant differences, ESR1, RARRES2, and MAP7, were on chromosomes 6 and 7 (Figure 2E). Further analyses revealed correlations among these DEGs (Figure 2F).

Two machine learning algorithms, RF and SVM, were used to screen telomere signature EM genes, plot box plots, and reverse cumulative distributions of residuals (Figures 3A, B). ROC curves show that RF (AUC = 1.0) had better predictive power compared with SVM (AUC = 0.969) (Figure 3C).

Based on the RFM, an importance score >2.0 was used as the threshold value. A total of 6 genes were included in the characteristic gene model: MAP7, GRHL2, RARRES2, UCHL1, ESR1, and REV3L (Figures 3D, E). A nomogram model was constructed based on the six included signature genes to predict disease risk (Figure 3F). Next, we validated the accuracy and predictive power of the model using CIC and DCA, with the Hosmer–Lemeshow test (P = 0.828) (Figure 3G). The CIC (Figure 3H) and DCA (Figure 3I) show good model predictive power at risk thresholds from 0 to 0.92.

To validate the role of the six screened genes in EM, an EM mouse model was developed (see flowchart in Figure 4A). Fourteen days after intraperitoneal implantation, mice injected with endometrial fragments developed endometrioid lesions in the intestine, mesentery, and peritoneum. Adhesions and vascular formations around the endometriotic implants were also detected (Figure 4B). Endometrial and endometriotic lesion tissues were collected for immunohistochemical analysis. Immunohistochemistry results further confirmed our database-based analysis (Figure 4C).

To investigate the modification patterns of telomerase genes in EM, we performed an unsupervised consensus clustering analysis of 28 EM samples based on the expressions of 15 telomere gene regulators. Two telomere gene modification subtypes in EM were identified by setting K value ranges at 2–9 and selecting the optimal K = 2 (Figures 5A, B). Among them, clusters 1 and 2 contained 21 and 7 samples, respectively. PCA showed that these two subtypes could clearly distinguish the samples (Figure 5E). In addition, we identified 10 telomerase gene regulators, which were significantly differentially expressed in the two isoforms (Figures 5C, D).

The heat map of the immune cell correlation analysis is in Figure 5F. HMOX1 was positively correlated with macrophages, with the most significant correlation coefficient of 0.87. KLF2 was negatively correlated with immature dendritic cells, with the smallest negative correlation coefficient of −0.65. Further analysis of immune cell differences indicated that five different immune cells differed significantly between clusters A and B: CD56 bright natural killer cell, CD56dim natural killer cell, macrophage, natural killer cell, and regulatory T cell (Figure 5G).

GO, KEGG pathway analysis, and GSEA were used to determine telomere gene roles and their potential mechanisms in EM. DEGs analysis was performed for both cluster A and cluster B samples (logFC >1, Padj <0.05) (Figure 6A). For the 138 DEGs obtained, we performed GO and KEGG analyses. Critical genes were mainly associated with extracellular matrix organization, extracellular matrix structural constituent, proteoglycans in cancer, and other pathways (Figures 6B–F).

Further GSEA analysis showed that critical genes were enriched in Reactome cytokine signaling in immune system, Naba core matrisome, Reactome signaling by interleukins, Reactome cell cycle checkpoints, Reactome cell cycle mitotic, and other pathways (Figures 6G, H).

We further subtyped the samples according to DEG expressions. An unsupervised consensus clustering analysis was performed on 28 EM samples, based on cluster 1 and 2 DEG expressions. Two distinct isoforms of differential gene modifications were identified by setting K values in the range from 2–9 and selecting the best K = 2 (Figures 7A, B). Clusters 1 and 2 contained 18 and 10 samples, respectively. Figure 7C shows the expression heat map of the 138 DEGs in clusters 1 and 2. Of the 15 telomere signature genes, eight were significantly differentially expressed in clusters 1 and 2: ESR1 (P < 0.001), HMOX1 (P < 0.001), KLF2 (P < 0.001), RARRES2 (P < 0.001), UCHL1 (P < 0.01), GRHL2 (P < 0.01), CCNB1 (P < 0.01), and FOS (P < 0.05) (Figure 7D). Further immune infiltration analysis revealed significant differences in the expressions of five immune cells in DEG analysis: macrophage (P < 0.001), natural killer T cell (P < 0.05), natural killer cell (P < 0.01), regulatory T cell (P < 0.001), and T folic helper cell (P < 0.05) (Figure 7E).

Finally, PCA was performed on the telomere gene scores for each sample. Telomere gene scores differed significantly between telomere genetic subtypes and differential genetic subtypes (P < 0.001). In both subtypes, telomere gene scores were significantly upregulated in subtype A/subtype 1. These results further demonstrate the accuracy and reliability of our model (Figures 7F, G).

By screening the DEGs identified in the CMAP database, we identified potential therapeutic agents for EM treatment. Eleven potential small-molecule compounds were identified based on the screening criteria of normalized connective score <2 and log₁₀q >15. Drug mechanism analysis suggested that crizotinib, AZD-4547, C-646, docetaxel, rociletinib, ENMD-2076, AMG-232, dioscin, mibefradil, NNC-05-2090, and nutlin-3 have potential as novel EM treatments (Supplementary Table 1).