Microarray data were obtained from a public dataset (GSE51981), which was first processed and batch-corrected, followed by differential expression and pathway enrichment analyses to ensure agreement of data processing with previous literature ( Figure 1 ). Table 1 describes the study population consisting of 105 samples across various disease stages (34 control, 24 stage I–II, 47 stage III–IV) and cycle phases (47 proliferative endometrium (PE), 24 early secretory endometrium (ESE), 34 mid-secretory endometrium (MSE)).

Batch correction successfully mitigated laboratory-associated variations ( Figure S1 ). Based on PCA dimensionality reduction plots, samples tended to cluster by disease/stage (case) versus control and by phase ( Figure 1B ). Heatmap clustering, focusing on genes highlighted as differentially expressed (FDR < 0.05, and log2 fold change >1) in any sample stratification, revealed strong clustering based on case versus control status ( Figure S2A ). Phase-stratified analysis revealed overall concordance in the results ( Figures S2B – D ). Several differentially expressed genes were identified across multiple phases, but also some phase-specific associations with endometriosis ( Figures S2B, C ). Among significantly upregulated genes, 79 were common to all menstrual phases such as FOSB, FOS, JUNB, and EGR1, and 182 were shared across at least two phases. In addition, there were 27 genes specific to ESE, 106 to MSE, and 428 to PE. Among the significantly downregulated genes, 246 genes were common across all menstrual phases, including CTSZ, SNTN, AGR3, and OLFM4, and 693 genes were shared across at least two phases. In addition, there were 64 genes specific to ESE, 201 specific to MSE, and 962 specific to PE. Stratified analysis allowed the identification of phase and disease stage where dysregulation occurred ( Supplementary Table S1 ). For example, CTSW, with functional roles in natural killer cells and cytotoxic T-cells and identified as upregulated in disease in previous unstratified analysis (4), was upregulated only in MSE in our analysis. The endothelial cell-related gene HSP90B1 was only elevated in PE, MSE, and stage I–II endometriosis. Immune-associated genes GNLY and C1QA were upregulated only in stage I–II and MSE, respectively. Discordance across phases is consistent with previous reports (7, 28); however, it is important to note that arbitrary differentially expressed genes (DEG) cutoffs may have amplified some of these differences as fold–fold plots revealed a high degree of correlation between phase-stratified disease versus control DEG results ( Figure S2D ). Similarly, concordant DEG signatures were also observed when data were stratified by disease stage ( Figure S3 ).

Pathway analysis, conducted with GSEA targeting MSigDB’s Hallmark Pathways, confirmed the general concordance noted with disease versus control DEG fold changes across stratifications and also recapitulated known biology in that many immune pathways were among those associated with disease. Interferon alpha and gamma responses, TGF-beta and IL-2 STAT5 signaling, and complement pathways were upregulated in endometriosis consistently across most menstrual phase and disease stage stratifications ( Figure S4 ; FDR < 0.05). TNF alpha signaling and allograft rejection, on the contrary, were significantly downregulated in disease across all cycle phases (FDR < 0.05).

Cell-type deconvolution provides a powerful opportunity to computationally disentangle bulk transcriptomic data into individual cell types. After confirming the agreement between our data processing with previous literature, we turned toward applying this technique by adapting a gene set enrichment-based deconvolution method (26) for use in the human endometrium ( Figure 1 ). The original method provided a comprehensive signature compendium for 64 classical human cell types derived from multiple organ types based on 6 human tissue consortia. To allow for discovery of new cell types and relationships between cell types, rather than using prior knowledge to select a subset of the signatures for analysis, we opted to utilize all signatures and rely on statistical metrics to infer likelihood that a cell type, of the 64, may be present in the endometrium.

As a summary, we first evaluated the statistical significance of xCell results using a permutation test. We then evaluated the specificity of xCell’s signatures toward endometrial cell types using a published scRNAseq dataset on healthy human endometrium. For in-depth evaluation on xCell’s signatures for immune subtypes, we performed immune cell only scRNAseq analysis on the scRNAseq dataset and examined the relationship between xCell’s immune subtype signatures and identified endometrial immune subtypes. Lastly we applied xCell enrichment analysis to artificial mixtures constructed from sorted endometrial cells to further validate our approach.

With metrics built for evaluating the applicability of xCell’s 64 signatures to the human endometrium, we deconvoluted 105 human endometrium bulk transcriptomic profiles into the 64 cell types using the adapted xCell approach ( Figure 1A ). Signatures that passed the permutation test were retained for downstream analysis. The data were obtained in proliferative, early secretory, and mid-secretory phases of the menstrual cycle from women with stage I–II or more severe (stage III-IV) endometriosis, as well as those without disease (control) ( Table 1 ). The overall clustering of the deconvoluted dataset was largely explained by disease versus control, followed by the menstrual cycle phase ( Figure 1B , right). Cellular compositions that contributed to the changes in phase and disease were then assessed via stratified differential enrichment analysis (FDR < 0.05, no FC cutoff, Figure 4 , Table 2 ) and interpreted alongside signature specificity metrics.

Our signature evaluation confirmed the specificity of xCell’s signatures to most major non-immune endometrial cell types, including epithelial cells and endothelial cells, but not fibroblasts. We observed increased enrichment scores in PE with endometriosis for epithelial cells and endothelial cells. This is consistent with observations of increased endothelial proliferation in women with endometriosis and menorrhagia versus controls (36–38).

eMSC is an endometrial cell type that exhibits mesenchymal stem cell characteristics in vivo (39) and in vitro (9, 30). Based on different characterization metrics, this cell type has been referred to as mesenchymal stem cells (9, 30, 40), pericytes (35), perivascular cells (41), or smooth muscle cells (27), each of which is represented by a different xCell signature. Our evaluation using both single-cell data and artificial mixtures discovered the lack of specificity of xCell’s signatures (MSC, pericyte, and smooth muscle cell) to eMSC, especially with low eMSC abundance ( Figure 2C ), due to concurrent expression of these signatures in eSF. This observation confirms the close relationship between these two endometrial cell types and their common association with progenitor MSC and pericytes. eMSCs are implicated in endometriosis (9, 30), and future studies should use unique markers identified for this cell type, such as RGS5, GUCY1A2, and NOTCH3 (9, 27, 35).

xCell’s fibroblast signature did not pass the permutation test despite receiving a moderate ratioNext score and “onTarget” classification ( Figure 2B ). This discrepancy may be due to many factors. Firstly, expression levels of fibroblast signature genes that passed the thresholds for ratioNext calculations (i.e., adjusted p-value and log2(fold change), Materials and Methods) often showed only the borderline fold changes, and as Subramanian et al. explain, signatures of this nature can be expected to score poorly in a gene set enrichment-based method (42). Furthermore, seemingly unrelated signatures such as chondrocytes, astrocytes, smooth muscle, and MSC showed the highest enrichment for the eSF cluster of the single-cell data ( Figure 2B ), and artificial mixture analysis confirmed that one of these cognate cell types, eMSC, could contribute to fibroblast-signature enrichment scores ( Figures 2C,D ). Our signature evaluation method ( Figure 1A ) thus considered output of this signature with low confidence. It is known that eSF have their own unique phenotype that is distinct from other fibroblasts of the body in many ways, and that expression profile as well with hormonal changes in the endometrium (9, 27, 35). Future studies should use a published dataset (9, 27, 35) to identify signatures that are specific to human eSF.

Our scRNAseq analysis identified 13 transcriptomically distinct immune types/subtypes in healthy endometrium, which were previously broadly categorized into lymphocytes and macrophages (27). The use of a signature-based method (31) and classical cell-type markers allowed us to confidently annotate pDC and monocytes, which have not been confidently identified at single-cell resolution or functionally examined in the endometrium, but have been characterized collectively with cDC1 by marker IRF-8 (43). With confirmed applicability of xCell’s signature for both cell types in the endometrium, our deconvolution results revealed relative increases in pDC and monocytes during MSE in women with endometriosis ( Figures 4 , 5 ), suggesting likely involvement in inflammation. pDC have known involvement in the inflammatory response normally and in pathologic settings through interaction with vasculature and T cells (44). Increased monocytes are congruent with the increased expression of monocyte chemoattractant protein-1 (MCP-1) in the endometrium of women with endometriosis (45, 46). Moreover, we observed a greater increase in monocytes in women with stage III–IV endometriosis, which may contribute to a lower implantation rate and live birth rates compared to women with stage I–II disease (47).

We have previously shown that endometrial macrophages (M1 and M2) in endometriosis are predominantly pro-inflammatory (10). Phase- or disease-stratified abundance quantification of endometrial myeloid cell types, including macrophages, monocytes, and dendritic cells, however, is limited. Here we report the unique markers that discriminate diverse endometrial myeloid cell types/subtypes for future studies.

Our scRNAseq analysis of immune cells in healthy endometrium has identified immune cell types/subtypes that are beyond the definition of xCell’s 64 signatures, such as four NK cell subtypes, one CD8+ T-cell subtype, one CD4 expressing T-cell subtype, and a KIT+ lymphocyte cell type. Intriguingly, xCell’s Th1, Th2, Tgd, and CD8+ Tem signatures were more enriched in NK cell subtypes rather than the T-cell subtypes in the single-cell dataset ( Figures 5 , S7 ). Therefore, the xCell outputs of these signatures across the cycle and in endometriosis likely reflect the changes of endometrial NK cell abundance or phenotypes more than the changes in the nominal cell types. However, these results do not conclusively suggest that these cell types are not present in the endometrium with or without disease, especially considering Th1 involvement in cytokine secretion and Tgd intraepithelial presence. Rather they are likely low in abundance, and their transcriptomic signals may be interfered by those of the more abundant NK cells. Lastly, although Tregs and cDC were identified in the scRNAseq dataset of healthy endometrium and demonstrated moderate enrichment of associated xCell signatures ( Figures 3 , S7A, B ), their xCell signatures did not pass our permutation test. Previous studies have identified both cell types in eutopic endometrium (43, 48–51) and have shown that Tregs increased in abundance during the secretory phase in women with endometriosis compared to controls in both the eutopic endometrium (52, 53) and peritoneal fluid (54) with potential interplay with eSF (55). Enriching for aforementioned cell types/subtypes with classical markers and single-cell level identification and analysis is warranted in future studies.

Notable for some xCell signatures that passed our permutation analysis are the elevated abundance scores of eosinophils, neutrophils, basophils, NKT, and immune progenitors in the endometrium of women with endometriosis and their absence in women without disease and in the annotations of the scRNAseq dataset of healthy endometrium. Eosinophils, initiators of inflammatory responses, were enriched in all phases of the cycle, compared to the control endometrium wherein they appear mainly during menses, confirmed herein and by others (56). Thus, eosinophils likely contribute to the pro-inflammatory phenotype observed in bulk-tissue analysis of the endometrium from women with endometriosis. Our finding of neutrophils, key participants in the innate immune response to foreign pathogens and enriched in the endometrium of women with endometriosis and independent of the cycle phase, compared to controls, is consistent with other reports, although others found cycle dependence of this cell population in women with versus without disease (57, 58). Basophils also initiate inflammatory responses and were found herein to be enriched in the endometrium of women with disease in PE and were significantly increased throughout the cycle. We are unaware of other reports on this cell type in the endometrium of women with endometriosis, and this finding warrants further study.

Our results generally agree with a prior deconvolution study on the endometrium from women with and without endometriosis (59), although fewer xCell signatures with disease-associated changes (9 in total) were identified compared to our study. Differences may be due to our adapted usage of the deconvolution method. Additionally, this study (59) did not design or apply metrics for statistical significance evaluation and result interpretation or develop de novo identification of normal endometrial immune cell signatures to enrich the xCell data interpretation in the endometrial context.

There are several limitations in this study. One is the limited sample size, especially in the ESE phase. Another limitation of our current approach is its limited capacity in inferring cell-type-specific phenotypic state. Although such insights can still be inferred for phenotype-specifying signatures, such as several immune cell subtypes mentioned above, for signatures without tissue-matching phenotype specifications, such as fibroblasts, such insights cannot be obtained directly from the deconvolution results. This is remarkable as the eSF changes transcriptomically across the cycle and displays marked abnormalities in endometriosis (9, 60) and is a key regulator of successful embryo implantation. Other cell-type deconvolution tools such as Cibersortx (15) provide the possibility to infer cell-type-specific gene expression profiles through additive combinations of input cell-type signatures. Successful application of this approach requires that highly specific cell-type signatures be used and that all potential cell-type signatures be included (61).

Further studies leveraging single-cell technologies as well as integrating different types of omics measurements including proteomics, epigenetics, and others will enable further corroboration of our findings and linking transcriptional phenotypes with endometriosis-associated cell types, especially considering the decreasing cost of single-cell analysis through strategies such as multiplexing. Functional studies will help elucidate the roles these cell types play in disease.

Through integrated whole-tissue deconvolution and single-cell analysis, we identified endometrial cellular compositions that are dynamic across menstrual cycle phases and altered in women with endometriosis. Guided by our signature evaluation metrics, we report cell-type candidates—immune cell types/subtypes of myeloid lineage, as well as non-immune cells, including epithelial and endothelial cell types—that most likely contribute to the pro-inflammatory endometrial phenotype previously observed in women with endometriosis (4, 7). Our results can help guide the selection of cell types for functional evaluation of cellular mechanisms that contribute to or result from endometriosis. Moreover, our analytical framework can be used in studies of other tissue types.

Background correction and quantile normalization were performed with the justRMA function of the affy package (62). Then batch correction was performed with ComBat (63) to reduce signals coming from the lab of origin, while protecting signals associated with the disease stage and the cycle phase. Direct principal component analysis (PCA) and the pvca package (64), which combines PCA with variance component analysis to estimate the proportion of variation in data that are associated with a set of potential sources, were used to assess the success of batch correction ( Figure S1 ).

Differential expression (DE) analysis was carried out via linear modeling using the Limma package (65) and the log-transformed and batch-corrected expression matrix as input. Simple, ~ single variable, formulas were used for linear model designs. When cutoffs for significant differential expression were used, they were FDR < 0.05, and abs log2 fold change > 1. Such analysis was run on various stratifications of the data ( Figures 1 , S2 , S3 ). For unstratified (all samples) analysis, DE was run 1) between disease versus control, 2) between stages I–II or stages III–IV versus control, and 3–5) between all combinations of pairwise cycle phase comparisons (PE vs ESE, PE vs MSE, ESE vs MSE). DE was also run between stages I–II versus stages III–IV, but zero genes met DE cutoffs. For stage-stratified analysis, DE was run separately on control samples only, stages I–II samples only, and stages III–IV samples only between all combinations of pairwise cycle phase comparisons. Lastly, for phase-stratified analysis, DE was run separately for PE samples only, ESE samples only, or MSE samples only, 1) between disease versus control and 2) between stages I–II or stages III–IV versus control.

Pathway enrichment analysis was performed on the Broad’s hallmark gene sets (obtained via https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp#H) by gene set enrichment analysis. This was carried out with the fgsea function of the fgsea package (66) on log2 fold changes of all genes (both significant and nonsignificant), for all stratifications and differential expression comparisons, with additional parameters: minSize = 15, eps = 0, and maxSize = 1500. Pathways with FDR corrected p-values below 0.05 were considered differentially enriched ( Figure S4 ).

Of the numerous deconvolution and enrichment methods, those that attempt to deconvolve a sample into additive mixtures of reference cell-type signatures have a strong reliance on both concordance between reference signatures and cell-type profiles of the target tissue, as well as on the presence of reference signatures for all cell types that might exist in the target sample. Given that we could not be certain that we would include signatures for all cell types that might exist in the endometrium, and that many immune cells profiled from the endometrium have shown noncanonical transcriptional profiles, we chose to use xCell’s enrichment-based approach, which is more robust to signature absence and inconsistencies, and includes signatures for 64 human immune and stromal cell types (including adaptive and innate immune cells, hematopoietic progenitors, epithelial cells, and extracellular matrix cells derived from thousands of expression profiles) (26). xCell was run on the log-transformed and batch-corrected expression data of human endometrial tissue described above. Due to uncertainty in the applicability of xCell signatures to endometrial tissue, only the rawEnrichmentAnalysis and transformScores steps were utilized for the calculation of enrichment scores. The spillOver adjustment step was not utilized due to notable deviations between xCell signatures versus nominal endometrial cell profiles, which would have been carried over into the compensation matrix derived from in silico mixtures of reference cells.

As discussed by its authors, xCell often produces nonzero scores, which may result in false-positive interpretation for nonexistent cell types or unsuitable signatures. xCell signatures that might not apply well to endometrial samples were first identified based on comparison to a permuted background distribution ( Figures S5 , S6 ). A background distribution of enrichment scores was generated for every cell-type signature, and for each cycle phase, by running xCell with 1,000 permutations of our expression matrix where rownames (gene symbols) were shuffled. Significance testing was then performed for each cycle phase among control or disease samples individually. Median enrichment scores of all iterations, among the current stratification samples, formed the background distribution for each cell type. For a given stratification, a cell-type signature was then considered as expressed if the “true” median enrichment score, from the nonpermuted data, was greater than the 90th quantile of its background distribution. To ensure that cells present in only certain conditions might still be accurately assessed, this filtering procedure was run on a per-disease status and per-phase basis, and cell-type signatures were retained for future analyses as long as the median enrichment score was above the background cutoff for at least one stratification.

Meanwhile, even when an xCell score is statistically significant, its contributing xCell cell-type signature may not be specific to its nominal cell-type target in the tissue of interest due to intertissue variability of the same cell type or ambiguity in cell-type naming.

To test sensitivity, for each of the 64 xCell signatures, a signature score was calculated with respect to each endometrial cell type identified in the scRNAseq dataset (27). To identify differentially expressed genes, Wilcoxon’s rank sum test (two-sided) was performed, and fold change (FC, dummy variable = 10^-2) was calculated between cells from an endometrial cell type and the remaining cells. P-values obtained from Wilcoxon’s rank-sum test were adjusted for multiple comparison by the Benjamini– Hochberg’s procedure to obtain p.adj. A signature score was quantified as the percentage of genes in the given xCell signature that were differentially expressed between cells in an endometrial cell type compared to the remaining cells (p.adj < 0.05, log2(FC) > 1). For each of the xCell signatures, the resulting score was normalized by the median of scores of all eight endometrial cell types identified in the single-cell dataset (normalized s_{i, j} = s_{i, j}/Med(s_i,₁, s_i,₂, … s_i,₈), where i is an xCell signature and j an endometrial cell type identified in the single-cell dataset).

Each xCell signature was categorized as either “reference” or “no-reference,” based on whether there is an endometrial cell type or subtype in the single-cell dataset that the signature is potentially targeting. A map between each xCell signature and each endometrial cell type was constructed to describe this relationship ( Figure 2B , boxes). As shown in Figure 1 , we kept this relationship relatively broad such that a signature is considered targeting a single cell type/subtype if it targets directly the identified endometrial cell type, or a subcategory of the identified cell type, or a related category of the identified cell type, to account for ambiguity in naming cell types and for potential existence of subtypes within the annotated cell populations.

Two specificity score metrics were then established. Given the target map, for the first specificity metric, “onTarget,” an xCell signature was tagged as “on-target” if the highest-ranking endometrial cell type from the single-cell expression data matches the cell type targeted by the xCell signature and “off-target” otherwise. Signatures without a clear reference cell type within the single-cell dataset were given an “NA” label ( Figure 2B ).

Separately, to evaluate how specific an xCell signature is to the endometrial cell type it represents, we calculated a “ratioNext” score representing the ratio between the highest and the second highest-ranking signature scores. Importantly, to avoid overpenalizing, if subtypes exist for the highest-ranking cell type (e.g., epithelial cells), scores in the subtypes were ignored in determining the second highest signature score ( Figure 2B ).

Dimensional reduction was performed on cells from the two clusters annotated as “Lymphocytes” and “Macrophages” in the original analysis (27) using Seurat’s (v3.2.0) (67) implementation of uniform manifold approximation and projection (UMAP). Specifically, top 2,000 variable genes among the immune cells were identified via FindVariableFeatures(). Principal component analysis was performed via RunPCA() on the top variable genes. Dimension reduction was performed on the top 20 principal components (PCs) via RunUMAP() based on the distribution of variances explained by the top PCs. Cell types/subtypes were identified using Seurat’s FindNeighbors(dims = 1:20) and FindClusters(resolution = 0.6). For each identified cell type/subtype, FindNeighbors() and FindClusters() were iterated one additional round to identify further heterogeneity. A cluster is classified as a candidate immune cell type/subtype if it can be defined by statistically significant uniquely expressing markers. pDC and B cells were present in 4 samples, macrophages were present in 7 women, and the rest of identified immune types/subtypes were in all 10 women.

For each identified cell type/subtype, uniquely expressing genes were found via FindAllMarkers(only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25, test.use = “wilcox”, slot = “data”) and ordered based on log2FC.

As elaborated in the text, annotation of each immune cell type/subtype was performed through iterative evaluation of classical marker expression, signature level scoring of xCell’s immune cell signatures, and RNA expression pattern of uniquely expressing genes identified above reported by the Human Protein Atlas (32). For signature scoring, we used the method reported in (31). Briefly, for each xCell’s immune signature, the score was quantified as the ratio between transcripts (UMI) that encode genes in the signature to all transcripts (UMI) detected in each single cell. We further examined the distribution of each signature in each identified immune cell type/subtype.

Microarray expression data from four cell types of FACS-purified endometrial cells (from participants with and without endometriosis) were used to generate 20 different artificial mixtures with varying proportions of each cell type. The microarray expression data from endothelial cells (n=11 samples), epithelial cells (n=7), mesenchymal stem cells (n=28), and stromal fibroblasts (n=31) were first summarized by their median expression for all probes. These median cell profiles were then additively combined into 20 different mixtures in which one cell type made up 10, 30, 50, 70, or 90% of the mixture, and the remaining cell types made up the remaining 90, 70, 50, 30, or 10%, respectively. xCell was then run on these mixtures both with and without the spillOver step ( Figures 2C, D ).

Log2 enrichment ratios (log2ER), between groups, were calculated for each cell type signature. P-values were generated by performing two-sided Mann–Whitney U tests between enrichment scores of all samples, between groups. These were then corrected for multiple hypothesis testing via the FDR method based on the number of signatures assessed. FDR < 0.05 was the sole cutoff used for differential cell-type enrichment. Such analysis was run on the same stratifications of the samples and for the same comparisons for each of those stratifications, as for differential gene expression analysis, described previously ( Figure 1 ).