Comparative Analysis of Shapley Values Enhances Transcriptomics Insights across Some Common Uterine Pathologies

1. Introduction

Disorders affecting the uterus represent significant burdens on women’s health worldwide. These conditions, characterized by aberrant cellular proliferation and tissue growth within the uterine environment, manifest with diverse clinical presentations and pose substantial challenges in diagnosis and management [1]. Endometriosis is defined as a debilitating gynecological disorder characterized by the presence of endometrial-like tissue outside the uterine cavity. This chronic widespread condition affecting up to 10% of women of reproductive age commonly causes pelvic pain, infertility and menstrual irregularities, among other symptoms [2]. Uterine leiomyomas, often known as fibroids, are benign smooth muscle tumors that arise within the uterine wall and are prevalent in up to 70% of women by the age of 50, often causing symptoms such as abnormal uterine bleeding, pelvic pressure, and reproductive dysfunction [3]. Endometrial cancer, originating from the malignant transformation of the endometrial lining, ranks as the most common gynecologic malignancy in developed countries, with an increasing incidence observed globally [4]. Despite the high prevalence of these common gynecological conditions and the ongoing debate about the existence of a genetic overlap and comorbidity among them, the molecular basis of these pathologies has yet to be determined [1,5]. Thus, understanding the molecular mechanisms underlying uterine pathologies is crucial for the development of targeted therapeutic interventions and improved patient outcomes [6].

Advances in omics technologies, particularly in microarray analyses, have paved the way to the comprehensive exploration of gene expression patterns associated with uterine conditions [7,8,9,10,11,12,13]. Microarray technologies provide the measurement of the expression levels for thousands of genes at a glance, which allows to obtain a deeper insight into the dysregulated molecular pathways implicated in the pathogenesis of several diseases [14,15,16,17,18]. The identification of differentially expressed genes (DEGs) represents a keystone in microarray data analysis. Classical approaches for microarray data analyses usually apply Welch’s t-test as a statistical method for the identification of DEGs by comparing expression levels between two experimental groups or conditions [19,20]. However, traditional methods may overlook significant changes at gene expression level, especially in complex diseases such as those affecting the uterus, which possess heterogeneous molecular profiles [21,22]. To overcome this limitation, alternative methods can be applied to help untangle the complexity underlying transcriptomics datasets [23,24,25].

In the present study, we aim to investigate the gene expression profiles associated with three of the most common uterine pathologies through the application of two different methods for microarray data analysis: i) a conventional method using Welch’s t-test and ii) a complementary analysis based on the Comparative Analysis of Shapley value (CASh) method derived from Game Theory, a method we have previously demonstrated that significantly increases the power to identify DEGs [25].

2. Materials and Methods

2.1. Microarray Expression Data Acquisition and Processing

Microarray data were obtained from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). For the selection of datasets of interest, raw data from Affymetrix commercial microarrays Affymetrix Human Genome U133A Array (HG-U133A), Affymetrix Human Genome U133A 2.0 Array (HG-U133A_2), Affymetrix Human Genome U133 Plus 2.0 Array (HG-U133_Plus_2), and Affymetrix Human Gene 1.0 ST Array [transcript (gene) version] (HuGene-1_0-st), were accessed preferentially, when possible.

CEL files from two datasets of endometriosis (GSE7846, GSE17504), two datasets of uterine leiomyomas (GSE12814, GSE23112), and two datasets of endometrial cancer (GSE36389, GSE63678) were retrieved from GEO repository. Raw data were downloaded for each dataset and preprocessing, quality control and normalization based on relative log expression (RLE), normalized unscaled standard error (NUSE), and Robust Multi-Array Average expression measure (RMA) methods were computed using ‘affy’ and ‘affyPLM’ packages in RStudio (version 2021.09.0) [26,27,28]. Finally, expression matrices were generated and samples were classified into experimental and control groups for further analyses (Supplementary Table S1).

Each dataset was processed independently in order to identify DEGs. To conduct differential expression analyses between patients and controls, two approaches were performed: i) a conventional approach based on the utilization of Welch’s t-test, and ii) an alternative method rooted on CASh technique.

2.1.1. Classical Approach

Conventional analyses for the detection of DEGs were performed using the Welch’s t-test implemented in the ‘multtest’ package in RStudio (version 2021.09.0) [29].

The significant DEGs were detected after multiple testing correction using the Benjamini & Hochberg method to control for False Discovery Rates (FDR) [30]. A significance threshold of an adjusted p-value (FDR) < 0.05 was applied.

2.1.2. Alternative Approach

We applied the Comparative Analysis of Shapley value (CASh) method to identify DEGs based on their cooperative contribution to overall gene expression changes [31]. The Shapley value, a concept derived from Game Theory, quantifies the marginal contribution of each gene to the collective expression change observed in the dataset [32]. CASh is a statistical technique that combines the Microarray Game algorithm (applied to transcriptomic values obtained from microarray experiments) with the Bootstrap technique, that applies random resampling of certain values, aiming to compensate for potential outliers in the data matrix [31,33,34,35]. Therefore, CASh considers gene expression as a cooperative game, where each gene contributes to the observed expression changes in a collaborative manner.

In our study, CASh method was applied to the detection of DEGs using two levels of restriction by stablishing 0.01 (more restrictive) and 0.05 (less restrictive) as cutoff p-values. Bootstrap resampling with 1000 iterations was computed in each analysis. In addition, Fold Changes (FC) were evaluated. Genes with p-values below 0.01 and 0.05 and |FC|>1.5 were considered as statistically significant.

2.1.3. Gene Set Enrichment Analysis and Functional Annotation

g:Profiler functional profiling tool g:GOSt (https://biit.cs.ut.ee/gprofiler/gost) was used to determine the enriched biological processes and pathways modulated by the DEGs based on Gene Ontology (GO) terms [36,37].

When needed, transcripts IDs were annotated and converted into official gene symbols using g:Convert tool (https://biit.cs.ut.ee/gprofiler/convert) from g:Profiler webserver [36,37]. For ambiguous transcripts names, IDs with the most GO annotations were utilized.

Three main GO categories were evaluated: biological processes (BP), cellular components (CC), and molecular functions (MF). GO terms with a Benjamini-Hochberg FDR <0.05 were considered to be significantly enriched. Top-ten significantly enriched GO terms in each category were plotted for CASh 0.05 comparisons using ‘ggplot2’ RStudio package [38].

3. Results

3.1. Datasets and Samples Analyzed

Gene expression data from six datasets covering a total of 68 samples was accessed. Table 1 describes the main characteristics of the datasets included in our study.

Datasets were analyzed for the detection of DEGs using two different strategies. First, the use of conventional methods based on Welch’s t-test was applied. Then, an alternative analysis based on CASh method was performed. The use of Welch’s t-test did not allow us to identify any DEGs, while several transcripts were revealed when using CASh method with both 0.01 and 0.05 cutoff p-values (Table 2). Total lists of DEGs for each dataset detected after each comparison are shown in Supplementary Table 2. Our analyses revealed that the application of CASh method allows a better detection of differentially expressed genes in the six datasets analyzed.

3.2. Functional Enrichment Analysis of the Differentially Expressed Genes

Given the restrictive criteria applied when running CASh 0.01 method, the number of DEGs detected did not allow to obtain a number of significantly enriched pathways associated to some genesets (data not shown). However, functional enrichment analysis of the differentially expressed genes obtained after the application of CASh 0.05 method revealed relevant significantly enriched processes in the analyzed datasets. In endometrial cancer datasets (GSE36389 and GSE63678), DEGs were mainly related to biological processes (BP) such as development and morphogenesis, whereas cellular components (CC) and molecular functions (MF) were mainly associated with extracellular locations and binding, respectively (Figure 1).

Regarding datasets of endometriosis (GSE7846 and GSE17504), the top significantly enriched BP were related to development, regulation of several cellular processes and morphogenesis. CC results revealed extracellular regions, cell surface and vesicles to be significantly relevant, and MF analysis detected functions mainly associated to protein activity (Figure 2).

Gene set enrichment analysis of the differentially expressed genes obtained after the application of CASh 0.05 method in uterine leiomyomas datasets (GSE12814 and GSE23112) revealed the regulation of several biological processes as significantly enriched BP, while membrane and binding processes were detected as significantly enriched CC and MF, respectively (Figure 3).

4. Discussion

Uterine pathologies impact women’s health and quality of life considerably. In recent years, the advent of omics technologies has facilitated a comprehensive exploration of molecular patterns associated with some of the most common gynecological conditions [39,40,41,42]. Microarray technology emerged about three decades ago with the aim of studying whole gene expression profiles, and the analysis of the amount of data derived from the application of this powerful tool has provided unprecedented insights into the discovery of dysregulated molecular pathways implicated in disease pathogenesis [43,44]. In the present study, we analyzed data from six datasets generated from the application of Affymetrix microarray devices: two datasets from endometrial cancer, two datasets from endometriosis and two datasets from uterine leiomyomas.

Raw data were downloaded from GEO public repository, and gene expression files were pre-processed, quality controlled and normalized. For the detection of DEGs, two strategies were adopted: i) a traditional approach based on the use of classical statistical t-tests and ii) an alternative approach using CASh method [31]. We were not able to detect any DEGs using traditional approaches, while the use of CASh method revealed a number of statistically significant genes in the six datasets analyzed. The t-test selects genes according to their differential expression between the two study conditions at an individual level. Thus, genes are considered significant when its p-value is below a stablished threshold (0.05 adj. p-value in our study). On the other hand, CASh method considers not only the expression of each gene under two conditions but the contribution of those genes that consistently interact with other genes, using the Shapley value to measure this contribution. Therefore, CASh highlights the most relevant genes as those that not only explain the differences between two conditions, but also play an important role over all possible permutations of genes [31,33,34,35]. In brief, CASh offers a more nuanced understanding of gene interactions and their collective impact on disease pathophysiology.

Interestingly, functional enrichment analysis of the DEGs detected using CASh method confirmed previous findings on the molecular bases of the uterine pathologies analyzed in our study. Some processes related to cell cycle and proliferation events have been shown to be significantly dysregulated in our sets of DEGs. Given the nature of endometrial cancer and endometriosis, it is plausible to believe that alterations at gene expression levels in some genes involved in these proliferative pathways may contribute to the phenotype of these diseases, as it has been previously proposed [45,46]. Further, a possible role of the degradation and remodeling of the extracellular matrix in endometriosis datasets has been revealed in our study. Endometriotic tissues have been shown to be significantly associated to extracellular matrix reorganization in some studies, which may explain some of the molecular mechanisms underlying the progression of the disease [47,48,49,50]. Regarding uterine leiomyomas, we were able to detect some significantly enriched biological processes that have been previously reported in association with the disease such as hormone secretion and cell signaling [51].

Our preliminary results underscore the potential of CASh as a valuable tool for analyzing microarray data. Further extensive research, including validation studies on larger cohorts and functional assays, is warranted to confirm the robustness and clinical relevance of the identified molecular signatures.

5. Conclusions

The Comparative Analysis of Shapley value seems to improve the detection of differentially expressed genes from microarray datasets compared to classical statistical techniques based on the use of t-tests in the context of complex diseases such as those affecting the uterus. By identifying key DEGs and elucidating their functional roles and interactions, we seek to advance our understanding of the underlying mechanisms driving disease origin and progression, thus paving the way for the development of targeted therapeutic strategies tailored to the management of women with uterine pathologies.