Gene expression analysis in endometriosis: Immunopathology insights, transcription factors and therapeutic targets

Abstract

Background: Endometriosis is recognized as an estrogen-dependent inflammation disorder, estimated to affect 8%-15% of women of childbearing age. Currently, the etiology and pathogenesis of endometriosis are not completely clear. Underlying mechanism for endometriosis is still under debate and needs further exploration. The involvement of transcription factors and immune mediations may be involved in the pathophysiological process of endometriosis, but the specific mechanism remains to be explored. This study aims to investigate the underlying molecular mechanisms in endometriosis.

Methods: The gene expression profile of endometriosis was obtained from the gene expression omnibus (GEO) database. Gene set variation analysis (GSVA) and gene set enrichment analysis (GSEA) were applied to the endometriosis GSE7305 datasets. Cibersort and MCP-counter were used to explore the immune response gene sets, immune response pathway, and immune environment. Differentially expressed genes (DEGs) were identified and screened. Common biological pathways were being investigated using the kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis. Transcription factors were from The Human Transcription Factors. The least absolute shrinkage and selection operator (Lasso) model identified four differential expressions of transcription factors (AEBP1, HOXB6, KLF2, and RORB). Their diagnostic value was calculated by receiver operating characteristic (ROC) curve analysis and validated in the validation cohort (GSE11691, GSE23339). By constructing the interaction network of crucial transcription factors, weighted gene coexpression network analysis (WGCNA) was used to search for key module genes. Metascape was used for enrichment analysis of essential module genes and obtained HOXB6, KLF2. The HOXB6 and KLF2 were further verified as the only two intersection genes according to Support Vector Machine Recursive Feature Elimination (SVM-RFE) and random forest models. We constructed ceRNA (lncRNA-miRNA-mRNA) networks with four potential transcription factors. Finally, we performed molecular docking for goserelin and dienogest with four transcription factors (AEBP1, HOXB6, KLF2, and RORB) to screen potential drug targets.

Results: Immune and metabolic pathways were enriched in GSVA and GSEA. In single sample gene set enrichment analysis (ssGSEA), most immune infiltrating cells, immune response gene sets, and immune response pathways are differentially expressed between endometriosis and non-endometriosis. Twenty-seven transcription factors were screened from differentially expressed genes. Most of the twenty-seven transcription factors were correlated with immune infiltrating cells, immune response gene sets and immune response pathways. Furthermore, Adipocyte enhancer binding protein 1 (AEBP1), Homeobox B6 (HOXB6), Kruppel Like Factor 2 (KLF2) and RAR Related Orphan Receptor B (RORB) were selected out from twenty-seven transcription factors. ROC analysis showed that the four genes had a high diagnostic value for endometriosis. In addition, KLF2 and HOXB6 were found to play particularly important roles in multiple modules (String, WGCNA, SVM-RFE, random forest) on the gene interaction network. Using the ceRNA network, we found that NEAT1 may regulate the expressions of AEBP1, HOXB6 and RORB, while X Inactive Specific Transcript (XIST) may control the expressions of HOXB6, RORB and KLF2. Finally, we found that goserelin and dienogest may be potential drugs to regulate AEBP1, HOXB6, KLF2 and RORB through molecular docking.

Conclusions: AEBP1, HOXB6, KLF2, and RORB may be potential biomarkers for endometriosis. Two of them, KLF2 and HOXB6, are critical molecules in the gene interaction network of endometriosis. Discovered by molecular docking, AEBP1, HOXB6, KLF2, and RORB are targets for goserelin and dienogest.

Keywords: ceRNA network; endometriosis; goserelin and dienogest; immune disorder; transcription factors.

Figure 1

GSVA and GSEA analysis in exploration GSE7305 dataset for endometriosis. (A) The pathway enrichment by GSVA. (B) Selected top five pathways from GSVA enrichment. (C) Main pathways analysis by GSEA. C0: normal endometrium. C1: endometriosis.

Figure 2

Relationships between immune cells from Cibersort analysis based on GSE7305. (A) The abundance values of sixteen immune cells in each sample. (B) Eighteen types of immune infiltrates with correlation. The squares in the heat map represent correlation strength, orange represents positive correlation, and blue represents negative correlation. The darker the color, the more meaningful it is. The P value is shown. The darker the color, the smaller the P value, and the lighter the color, the larger the P value (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001).

Figure 3

Relative cellular fraction of 16 immune infiltrating cell types in endometriosis assessed by Cibersort (A). Eight immune cell types, endothelial cells, and fibroblasts in endometriosis evaluate by MCP-counter (B). Differential expression of immune response gene (C) and immune response signaling pathways (D) in endometriosis. *P < 0.05; **P < 0.01; ***P < 0.001; ns, not significant with P > 0.05.

Figure 4

GO and KEGG analysis for identified DEGs. DEGs between endometriosis and non-endometriosis were displayed in volcano plot (A) and the top 30 DEGs in the heat map (B). (C) The chord diagram and circle diagram showed the enrichment analysis of DEGs.

Figure 5

TFs in DEGs and TF-TF interactions. (A) Overlap between DEGs and TFs. (B) The enrichment analysis of GO and KEGG for twenty-seven TFs regrouped from the intersection. (C) The correlation of twenty-seven TFs between each other.

Figure 6

TFs-immune interaction. (A) The relationship between the target twenty-seven TFs and immune-related genes, associated immune pathways (B), and immune infiltrating cells (C).

Figure 7

Lassol model screened TFs diagnostic values both in the training and validation cohorts. (A) Results of the Lasso multivariate model. (B) ROC curves of AEBP1, DLX65, HOXB6, KLF2, and RORB in the GSE7305 dataset. (C) ROC curves of AEBP1, DLX65, HOXB6, KLF2, and RORB in the GSE11691 dataset. *P < 0.05; **P < 0.01; ***P < 0.001; ns, not significant with P > 0.05.

Figure 8

Correlation network generated for AEBP1, HOXB6, KLF2, RORB using String. The interaction network diagram for AEBP1, HOXB6, KLF2 and RORB respectively (A-D) and combined comprehensive network diagram (E). (F-I) The associated pathway from GSEA analyses for AEBP1, HOXB6, KLF2, and RORB.

Figure 9

WGCNA analysis of gene networks for closely related four genes (AEBP1, HOXB6, RORB, and KLF2). (A) The soft thresholds. (B) A tree of modules, each color represents a specific coexpression module and the upper branches represent genes. Genes that do not belong to any module will be marked grey. (C) The heat maps of different co-expression modules. (D) The turquoise module (uptrend) and the blue module (downtrend), respectively.

Figure 10

Enrichment analysis of interaction network by Metascape for sixty-two genes extracted from WGCNA, including GO (A) and KEGG (B).

Figure 11

Sankey diagram of the lncRNAs-miRNAs-mRNAs network. Each rectangle represents a gene, and the size of the rectangle indicates the degree of connectivity of each gene.

Figure 12

Results of molecular docking simulations. Protein-ligand docking pocket for each transcription factor (AEBP1, HOXB6, KLF2 and RORB). The pocket with the lowest binding energy when the transcription factor binds Dienogest or Goserelin. PHE, Phenylalanine residue; MET, methionine; HIS, Histidine residue; ARG, Arginine; ALA, alpha Linolenic acid; LEU, leucine; ASN, ASP, Aspartic acid.