Gene expression profiles of RIF were screened from the GEO (http://www.ncbi.nlm.nih.gov/geo) database. Inclusion criteria were as follows: 1) Homo sapiens expression profiling by the array; 2) samples were endometrium of RIF patients or control individuals (CON) during the window of implantation; 3) datasets contained ten or more samples with at least five patients in each group, and 4) RIF patients and fertile controls were included in one experiment. This study ultimately included four datasets sed on the above selection criteria, including GSE111974, GSE92324, GSE26787, and GSE71835. Details of all data are shown in Table 1.

We merged GSE111974 and GSE92324 microarray data as test datasets. Specifically, the first step is to convert the series matrix file from gene probe IDs to gene symbol codes, averaged for the case of one gene corresponding to multiple probes. The second step is to remove the batch effect, we first used limma’s removeBatchEffect function (Ritchie et al., 2015), yet it failed to eliminate the batch effect between GSE111974 and GSE92324 (Supplementary Figure S1). We then used sva’s combat function (Leek et al., 2012) to eliminate the batch effect between the datasets, and Supplementary Figure S2 showed that the batch effect was successfully eliminated for the merged data. This may be explained by the fact that the removeBatchEffect function removes known batch effects from the data (Ritchie et al., 2015), while sva package not only removes known batch effects but also adjusts for other potentially unwanted sources of variation in the data for subsequent analysis (Leek et al., 2012).

The final step is to normalize the expression values through the limma package to have a similar distribution in a set of arrays. Here, we have drawn up a flow chart of the analysis process ( Supplementary Figure S3 ).

The differentially expressed genes (DEGs) between RIF patients and CON were selected by using the limma package with the |Log2FC (fold change) | > 1 and adjusted p-value < 0.05. Analyses of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) for DEGs were performed by the clusterprofiler package (Yu et al., 2012), the significantly different GO terms were determined by thresholds adjusted p-value < 0.05, and KEGG pathways with a p-value <0.05 was selected. In addition, to discover candidate genes sets or pathways that likely contribute to RIF, gene set enrichment analysis (GSEA) (Subramanian et al., 2005) was performed by the clusterprofiler package to scrutinize the gene expression profile at an entire level, and C5 (ontology gene sets) was chosen for functional enrichment analyses. The normalized enrichment score (|NES| >1), p-value <0.05, and adjusted p-value <0.05 were set as threshold criteria.

In this study, 28 immune cell types and associated 782 marker gene signatures were first obtained from two previously reported studies (Supplementary Table S1) (Barbie et al., 2009; Charoentong et al., 2017). Then, the abundance of these immune cell types in endometrial samples was calculated by the single sample gene set enrichment analysis (ssGSEA) method based on the genomic variance analysis (GSVA) algorithm. Recently, the application of ssGSEA in deconvolution of bulk gene expression data has been widely performed (Zhao et al., 2021; Liu et al., 2022).

CIBERSORT is another method to calculate cell composition based on expression profiles. In the present study, we used CIBERSORT to assess immune cell infiltration in endometrial tissue between RIF and CON. The leukocyte signature matrix (LM22) was used as a reference expression signature with 1,000 permutations (Zhou et al., 2021). LM22 signature matrix contains 22 infiltrating immune cell components and the corresponding 547 signature genes (Supplementary Table S2) (Newman et al., 2015). Then the Wilcoxon test was conducted to determine significant differences in immune cell types between RIF and CON.

The weighted gene co-expression network analysis (WGCNA) is an algorithm that can find modules of a co-expressed gene with high biological significance (Langfelder and Horvath, 2008). In this study, to reduce the whole network’s computation size but maintain a scale-free topological network, we selected the genes in the top 75% based on the magnitude of the variance. Then we entered them into the WGCNA package in R to identify the gene modules associated with significantly altered immune cell types. Briefly, genes with similar expression patterns were assigned to co-expression modules by weighted correlated adjacency matrices and clustering analysis. Firstly, the weighted adjacency matrix is constructed by calculating an appropriate soft threshold β that satisfies the criteria for a scale-free network. Afterward, the weighted adjacency matrix was converted to a topological overlap matrix (TOM), and the corresponding dissimilarity degree (1-TOM) was generated. Then, module identification was performed using the dynamic tree-cutting method, and modules with differences less than 0.25 were merged. In addition, the relationship between module eigengene values and immune cell types was assessed by Pearson correlation.

Additionally, we performed the local maximal quasi-clique merger (lmQCM) to network mining (Zhang and Huang, 2014) based on the merged matrix by the lmQCM package in R. The parameters for lmQCM were set as follows: gamma = 0.55, t = 1, lambda = 1, beta = 0.4, and minimum cluster size = 10. The lmQCM is a weighted network mining algorithm that detects weak quasilinear modules in a weighted graph and applies it to the discovery of functional gene clusters. The algorithm is characterized by a greedy approach using hierarchical clustering that does not allow overlap between modules but allows genes to be shared between multiple modules. This is in accordance with the fact that genes are often involved in multiple biological processes (Bichindaritz et al., 2021).

Among the weighted network modules constructed from WGCNA and lmQCM, we selected the module with the highest or lowest correlation coefficient as the specific module associated with the differentially expressed immune cell types. This research defined the modules most relevant to γδ T cells and activated memory CD4 T cells as crucial modules. The genes contained in the crucial modules were defined as differentially expressed immune cell types-related genes (DE ICTRGs).

The overlapping genes between DEGs and DE ICTRGs were identified with the Venn online platform (http://bioinformatics.psb.ugent.be/webtools/Venn/) and defined as potential hub genes. If too many genes were overlapping, they would be further filtered according to the protein-protein interaction information from STRING (https://string-db.org) with confidence scores ≥0.4. The interaction file (string_interactions.tsv) was downloaded. Subsequently, ten algorithms of cytoHubba (Chin et al., 2014) in Cytoscape 3.9.0 (Shannon et al., 2003) were conducted to score each node gene, namely, MCC (Maximal Clique Centrality), MNC (Maximum Neighborhood Component), Degree, EPC (Edge Percolated Component), BottleNeck, EcCentricity, Closeness, Radiality, Betweenness, and Stress. Lastly, the ten node genes with the highest scores for each algorithm were examined for hub genes using the UpSet package in R. Furthermore, GeneMANIA (http://genemania.org) analyses were performed to examine protein and gene interactions, pathways, co-expression, co-localization, and protein domain similarities (Franz et al., 2018).

To further validate the accuracy and reliability of the hub genes selected from test datasets, two external datasets, GSE26787 and GSE71835 microarray data, were downloaded from the GEO database and combined using the approach mentioned above. Firstly, the expression of hub genes was extracted from the test sets and validation sets and analyzed by the Wilcoxon test, with a p-value of <0.05 defined as statistical significance. Then, we constructed a prediction model using the differentially expressed hub genes by the generalized multivariate regression with the test sets. Finally, we calculated the sensitivity rate and specificity rate of the model, Receiver Operating Characteristic (ROC) analysis was also performed to detect the Area Under the Curve (AUC).

The multiMiR package in R is a comprehensive collection of predicted and validated miRNA–target interactions and their associations with diseases and drugs (Ru et al., 2014), including 14 databases. In the present study, we selected three databases (DIANA-microT, miRanda, and TargetScan) to predict the targeted miRNAs of the hub genes. The top 35% of miRNAs in the prediction scores of these three databases were intersected. The intersection was used for subsequent analysis. For all long non-coding RNAs (lncRNAs)-miRNA interaction data were acquired in the starbase database (https://starbase.sysu.edu.cn/) (Li et al., 2014), and the target lncRNAs were filtered according to clipExpNum >4. Eventually, the established network was visualized by Cytoscape software.

In this study, all data analysis and visualization were performed using R software (version 4.0.5; https://www.r-project.org/) with appropriate packages. p-value <0.05 was considered significant.

We downloaded GSE111974 and GSE26787 from the GEO database as test datasets, including 34 RIF and 32 normal endometrial tissues. Box plots and principal component analysis shows the data before batch correction (A, C, and E) and after batch correction (B, D, and F) (Supplementary Figure S2), which indicates that the batch effect was successfully eliminated from the combined data.

A total of 324 DEGs were identified (223 significantly up-regulated genes and 101 significantly down-regulated genes) (Figures 1A,B; Supplementary Table S3). Then, the DEGs were subjected to GO and KEGG pathway enrichment analyses. The GO enrichment revealed that these DEGs were mainly associated with carboxylic acid transport, organic acid transport, and detoxification in the biological process. Regarding cellular components, the genes were primarily enriched in the apical plasma membrane, apical part of the cell, and collagen-containing extracellular matrix. As for molecular function, the genes were enriched primarily in anion transmembrane transporter activity, organic anion transmembrane transporter activity, and active transmembrane transporter activity (Figure 1C; Supplementary Table S4). Likewise, the KEGG analysis demonstrated that these DEGs are relevant to immune pathways such as TNF signaling pathway, Leukocyte transendothelial migration, and NF-kappa B signaling pathway (Figure 1D; Supplementary Table S4). In addition, the GSEA results showed that 385 gene sets were significant at an adjusted p-value <0.05, and most of the enriched gene sets were related to various immune responses (Supplementary Table S5). Figure 1E shows the five most enriched immune-related gene sets based on the adjusted p-value. These were related to leukocyte cell-cell adhesion, leukocyte migration, T cell activation, positive regulation of lymphocyte activation, and antigen receptor-mediated signaling pathway.

Next, we explored immune cell changes in the test set. First, ssGSEA identified 16 immune cell subtypes, including activated CD8 T cells, activated dendritic cells, CD56 dim natural killer cells, central memory CD4 T cells, central memory CD8 T cells, and γδ T cells. Their cell-specific marker genes were lower expressed in the RIF group (Figure 2A). Conversely, the specific marker genes of three immune cells (CD56 bright natural killer cells, effector memory CD4 T cells, and eosinophils) were expressed at higher levels in the RIF group. Furthermore, compared with the CON group, CIBERSORT analysis demonstrated that activated memory CD4 T cells had statistically higher abundance. In comparison, γδT cells had statistically lower abundance in the RIF group (Figure 2B), which is consistent with the results of ssGSEA. The above results indicated that γδT cells and activated memory CD4 T cells were the significantly altered cell types in the RIF group. In addition, the constituency of the 22 immune cell types in each sample was plotted as a histogram by performing CIBERSORT (Figure 2C). Meanwhile, the correlation between these 22 immune cell types in endometrial tissue from the RIF group was calculated (Figure 2D). Figure 2D shows a significant positive correlation between memory resting CD4 T cells and monocytes (R = 0.61). Also, plasma cells were positively correlated with resting mast cells (R = 0.57). In contrast, CD8 T cells were negatively correlated with macrophage M0 (R = -0.71). Likewise, γδT cells were negatively associated with activated memory CD4 T cells (R = -0.48). Nevertheless, the association mentioned above of immune cells was attenuated to null in CON (Supplementary Figure S4).

In WGCNA analysis, the 13,194 genes in the top 75% based on the magnitude of the variance were included in the WGCNA analysis, and the soft power of β = 10 (scale-free R ² > 0.85) was determined as soft-thresholding to acquire co-expressed gene modules (Supplementary Figures S5A, B). Then, dynamic hybrid cuts were conducted to construct hierarchical clustering trees by dividing the dendrogram at relevant transition points (Supplementary Figure S5C). Of which, single genes were represented as tree leaves, multiple genes with analogous expression data were presented as branches of the dendrogram tree, and branches containing similarly expressed genes were considered gene modules. Similarly, another differentially expressed immune cell types co-expression network was also constructed by lmQCM analysis. Ultimately, we got 14 WGCNA modules (Figure 3A) and 15 lmQCM modules (Figure 3B).

Among these, the dark olive-green module (correlation = 0.31, p-value = 0.01) was the most relevant module identified by WGCNA for activated memory CD4 cells, and the salmon module (correlation = -0.29, p-value = 0.02) was the most relevant module identified by lmQCM for activated memory CD4 cells. Thus the dark olive-green module was selected as the activated memory CD4 cell-associated key module for further analysis. Similarly, the dark turquoise module (correlation = 0.35, p-value = 0.004) was the most relevant module identified by WGCNA for γδT cells, and the green module (correlation = 0.3, p-value = 0.01) was the most relevant module identified by lmQCM for γδT cells. Thus the dark turquoise module was selected as the key module associated with the γδT cells module for further analysis. In this study, we defined genes in these two modules most relevant to γδ T cells and activated memory CD4 T cells as DE ICTRGs.

We collected shared genes from DE ICTRGs and DEGs using a Venn diagram. It turned out that two overlapping genes between γδ T cells-associated genes and DEGs (ALOX5AP, SLC7A7) (Figure 3C), and 213 overlapping genes between activated memory CD4 T cells-associated genes and DEGs (Figure 3D) (Supplementary Table S6). Considering there were too many overlapping genes in Figure 3D, the PPI network was constructed for 213 genes and filtered by cytoHubba in Cytoscape. The results of the CytoHubba were listed in Supplementary Table S7, and PTGS2 was determined as the hub gene (Figure 3E). Finally, we explored three hub genes (ALOX5AP, SLC7A7, and PTGS2) and their 20 interacting genes using the GeneMANIA database (Figure 4). The network illustrated that these genes were relevant to immune processes such as the leukotriene metabolic process, antigen binding, and regulation of inflammatory response.

For validating the identified hub genes, another two datasets, GSE26787 and GSE71835, were merged after removing the batch effect (Supplementary Figure S6), including 11 RIF and 11 normal endometrial tissues. The expression levels of ALOX5AP, SLC7A7, and PTGS2 were presented in the heatmap (Figure 5A). As shown in Figures 5B,C, in both the test set and validation set, the expression levels of ALOX5AP and SLC7A7 were significantly decreased, and PTGS2 was increased dramatically in the RIF group (p < 0.05). In addition, we used the datasets GSE26787 and GSE71835 as validation sets to investigate the predictive effect of hub genes for RIF. Encouragingly, in the validation set, the prediction model showed a specificity of 90.8% and sensitivity of 90.8%, and the ROC analysis showed that the AUC was 0.908 (Figure 5D).

To explore possible interactions between lncRNAs, miRNAs, and mRNA in RIF, we structured a ceRNA regulatory network. In the present study, we collected 53 miRNAs, including hsa-miR-3180-3p, hsa-miR-548p, and hsa-miR-1297 (Figures 6A–C). Next, we mapped the abovementioned 53 miRNAs into the starbase database and searched for the target lncRNAs. As a result, 35 lncRNAs that interacted with 19 of the 53 miRNAs in the starbase database were selected (Supplementary Table S8). Eventually, the ceRNA regulatory network was structured, and the visualization was carried out in the CytoScape software (Figure 6D).