Enhancing Transcriptomics Insights into Neurological Disorders Through the Comparative Analysis of Shapley Values

1. Introduction

Complex disorders affecting neurological processes are responsible of great health, social and economic costs worldwide. Despite the heterogeneity of these complex disorders, they all pose a significant global burden, since the misunderstanding of their causes and associated factors that intensify the importance of these phenotypes is the main cause of the insufficiency of diagnosis, and also the lack of effectiveness in medical treatment for patients negatively impacts the well-being of those affected.

Autistic Spectrum Disorder (ASD) (Hirota & King, 2023; Lord et al., 2018) is a phenotype that contains from the most severe autism, when social and communicative functions are very limited, until Asperger syndrome, characterized by mild symptoms. Anyway, all diagnostic features show a rigid behavior, a pathological selection for some issues, and the capacity for attention and communication is affected (Sharma et al., 2018). Some body systems are also affected, as digestive (Martin et al., 2018; Srikantha & Mohajeri, 2019), immune (Heidari et al., 2021; Meltzer & Van de Water, 2017; Ormstad et al., 2018; Robinson-Agramonte et al., 2022), circulatory (Kealy et al., 2020; Mouridsen et al., 2016) and nervous (Doroszkiewicz et al., 2021; Matta et al., 2019; Sharon et al., 2016). Microbiota (Kang et al., 2017) and genetic causes (Qiu et al., 2022) have been proposed in the early development of this disorder, and studies support that risk of suffering ASD rises when relatives are affected (Arberas & Ruggieri, 2019). These symptoms are harmful for patients’ autonomy and the welfare of caregivers (Dissanayake et al., 2019; Ruzzo et al., 2019). The World Health Organization communicated in 2023 that one out of one hundred children are suffering this disorder and prevalence is rising in the previous few years. Due to that, and given the fact that origins and development of this condition are not agreed by specialists and researchers in this field of medicine (Silverman et al., 2010; Tseng et al., 2022), plenty of research teams are thinking of strategies in order to discover the etiology and main factors for understanding this disorder.

Schizophrenia is a neurological disorder that is characterized by positive (hallucinations, lack of social skills and cognitive distortions) and negative (general apathy, in social and job issues) (Khan et al., 2013) symptoms. This disease is linked to increased vulnerability to cardiovascular (Barnett et al., 2007; Hagi et al., 2021), metabolic (Trubetskoy et al., 2022) and infectious (Fuglewicz et al., 2017) diseases, which rise the risk of an early death. Furthermore, it has a direct link with suicide index growth (Balhara & Verma, 2012; Carlborg et al., 2010; Sher & Kahn, 2019), and also caregivers and nearly people are negatively affected in social terms, since patients suffer from diminished autonomy (Gulayín, 2022; Ribé et al., 2018). Its prevalence worldwide is 24 million people (McGrath et al., 2008; Saha et al., 2005), with a percentage of 0.32%, 0.45% in adults (World Health Organization, 2022), and it tends to appear in teenagers at an advanced age, probably with bonds to neural pruning (Germann et al., 2021). Up to date, the origin of this disorder remains unknown (Pino et al., 2014). In this line, there is some consensus in the relevance of some gene factors implied in its onset (Trubetskoy et al., 2022), but that are not determinant to its origin (Vilain et al., 2013), considering that other factors as social environment, drug abuse (including alcohol) (Häfner & an der Heiden, 1997; Janoutová et al., 2016; Stilo & Murray, 2019), and neural pruning, usual in adolescence, can influence in a decisive way (Germann et al., 2021). Myelin sheaths (Valdés-Tovar et al., 2022) and central nervous system architecture (Bobilev et al., 2020; Cheng et al., 2022; Heckers & Konradi, 2002) are also bonded with this disorder. The most extended belief nowadays is that schizophrenia is a multifactorial disorder (Kahn & Sommer, 2015; Morera-Fumero & Abreu-Gonzalez, 2013). In order to clarify the disease causes, omics techniques as transcriptomics have been applied (Gandal et al., 2018). Nevertheless, this is a pathological situation that harms in a severe way the life quality of patients, what generates a medical and social interest that concerns to pharma industry, which intends to alleviate this suffering with drugs that minimize the secondary effects associated to available treatments (Krause et al., 2018), usually adverse for daily life of patients (Perlick et al., 2010). Thus, efficient research is crucial to solve the social and economic problems attached to this disease (Evans-Lacko et al., 2014).

Bipolar disorder (BD) (Smith et al., 2012; Tondo et al., 2017) is a neurological condition characterized by the alternance of maniac episodes (euphoria, excessive joy, uncontrolled enthusiasm, etc.) with depressive ones (anhedonia, sadness, lack of interest for living, etc.) (Fagiolini et al., 2015). Genetic causes have been studied (Gandal et al., 2018), and some environmental factors as alcoholism and other types of drug abuse has been proposed as disease cause (Aldinger & Schulze, 2017). The development of genomics and transcriptomics may help to understand the disorder and treat it efficiently. Its prevalence was 40 millions of people in 2019 (Institute of Health Metrics and Evaluation, 2022), and the suffering rises suicide index for these patients (Beyer & Weisler, 2016; Clemente et al., 2015; Miller & Black, 2020). There is still not much understanding about this disorder, but some drugs, including lithium, have been reported to alleviate its symptoms (Katz et al., 2022; Malhi et al., 2017).

Major depressive disorder (MDD) (Filatova et al., 2021) is a neurologic disease of unknown origin (Gómez Maquet et al., 2020), with more severe symptoms that common depression (Kennedy, 2008). Among these are anhedonia, sadness and lack of desire to live (Bauer et al., 2019). Genetic causes are considered, which has led to the development of transcriptomics and epigenetic studies (Kendall et al., 2021), but also a physiological and hormonal origin are reported, as well as environmental factors like stress, psychological and social aspects (Suda & Matsuda, 2022). Due to the fact that its origin remains unknown, it is classified as a complex disorder (Harder et al., 2022), which causes a great social and economic burden for the community environment of the affected people (Greenberg et al., 2021; Keshavarz et al., 2022). World prevalence is about 350 million of people (Gutiérrez-Rojas et al., 2020; Smith et al., 2019), but there is not much consensus. In fact, this prevalence differs between regions (3% in Japan and 16.3% in USA) (Gutiérrez-Rojas et al., 2020; Smith et al., 2019). Every year even 850,000 suicides have been registered because of major depressive disorder (Li et al., 2022; Serra et al., 2022). Different techniques, such as omics ones and neuroimage, and several biomarkers, as certain fatty acids and miRNA have been used, but there is no consensus (Figueroa-Hall et al., 2020; Gadad et al., 2018; Zhou et al., 2019). Nowadays, there are lots of medicines that treat this disease, taking advantage of the limited knowledge we have about encephalon.

Despite high prevalence worldwide, the origins of these disorders are still unknown. Because of that, it is necessary to apply techniques able to detect key factors for prevention and treatment, pointing towards its main causes and improving as much as possible the health and quality of life of these patients.

Advances in omics technologies, particularly microarray analysis, have revolutionized the comprehensive exploration of gene expression patterns associated with neurological conditions (Bettencourt et al., 2023; Legati et al., 2021; Xu et al., 2022). Microarray technology enables the simultaneous measurement of thousands of genes, providing deep insights into the dysregulated molecular pathways implicated in the pathogenesis of various diseases (Bryant et al., 2004; Copland et al., 2003; Krokidis & Vlamos, 2018; Rai et al., 2016; Ward, 2006). A critical aspect of microarray data analysis is the identification of differentially expressed genes (DEGs), which serve as key indicators in understanding disease mechanisms. Traditionally, these analyses have relied on ranking genes based on individual p-values; however, this approach does not always correlate with biological significance. In some cases, very small p-values, indicative of high statistical significance, may not correspond to biologically relevant signals, while larger p-values, often disregarded, could be linked to genes crucial for specific biological processes (Esteban & Wall, 2011). Classical microarray analysis methods typically utilize Welch’s t-test and linear models such as Empirical Bayes to identify DEGs by comparing gene expression levels between experimental groups or conditions (Jeffery et al., 2006; Selvaraj & Natarajan, 2011). However, these traditional approaches may miss significant gene expression changes, particularly in complex diseases like those affecting the brain, which are characterized by heterogeneous molecular profiles (Ganapathy et al., 2019; Villani & Marzetti, 2023).

To address the limitations of p-value-based approaches, which often result in the excessive suppression of biologically relevant signals due to multiple testing correction methods, more robust methodologies have been developed (Breitling & Herzyk, 2005; Cordero et al., 2007; Esteban & Wall, 2011). Notably, one such method integrates Game Theory, utilizing a computational index known as the Shapley value (Esteban & Wall, 2011). This approach provides a more nuanced assessment of gene significance by evaluating the cumulative contribution of each gene within the context of the entire gene set analyzed. The Shapley value quantifies the importance of each gene by considering its contribution in conjunction with the contributions of all other genes in the same experiment (Moretti & Patrone, 2008). By combining Game Theory with traditional statistical analyses, this methodology offers a powerful tool for enhancing the detection and interpretation of meaningful gene expression differences (Esteban & Wall, 2011).

We applied the microarray games methodology in this study, specifically harnessing Shapley values, to analyze gene expression data related to various neurological pathologies. This approach integrates Game Theory to improve the detection and functional analysis of genes involved in complex biological conditions, such as ASD, schizophrenia, bipolar disorder, and major depressive disorder (Esteban & Wall, 2011). By evaluating the average marginal contribution of each gene across all possible coalitions, this technique reveals critical insights into the genetic underpinnings of these diseases, potentially leading to innovative diagnostic and therapeutic strategies. The game-theoretic approach not only enhances the identification of key genetic players but also enriches our understanding of their biological roles within complex, multi-genic pathologies.

To achieve a comprehensive understanding of gene expression profiles associated with four prevalent neurological pathologies, we employed two distinct methods for microarray data analysis: (i) a conventional approach utilizing Welch’s t-test and Empirical Bayes methods, and (ii) a complementary analysis based on the Comparative Analysis of Shapley value (CASh) method, derived from Game Theory. Previous research (Castro-Martínez et al., 2024; Esteban & Wall, 2011) has demonstrated that the CASh method significantly increases the power to detect differentially expressed genes (DEGs), providing a more robust framework for analyzing complex biological data.

2. Materials and Methods

2.1. Microarray Expression Data Acquisition, Processing and Exploratory Analysis

Microarray data were sourced from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). For dataset selection, raw data from Affymetrix commercial microarrays—specifically, the Affymetrix Human Genome U133A Array (HG-U133A), Affymetrix Human Genome U133 Plus 2.0 Array (HG-U133_Plus_2), and Affymetrix Human Gene 1.1 ST Array [transcript (gene) version] (HuGene-1_1-st)—were prioritized whenever available.

CEL files from three autism datasets—GSE6575, GSE18123, and GSE25507 (Alter et al., 2011; Gregg et al., 2008; Kong et al., 2012); two schizophrenia datasets—GSE17612 and GSE62333 (Cattane et al., 2015; Maycox et al., 2009); two bipolar disorder datasets—GSE5389 and GSE7036 (Matigian et al., 2007; Ryan et al., 2006); and two datasets encompassing schizophrenia, bipolar disorder, and major depressive disorder samples—GSE12654 and GSE53987 (Iwamoto et al., 2004; Lanz et al., 2019) were retrieved from the GEO repository. Raw data were downloaded for each dataset, followed by preprocessing, quality control, and normalization using relative log expression (RLE), normalized unscaled standard error (NUSE), and Robust Multi-Array Average expression measure (RMA) methods. These processes were carried out using the ‘affy’ (version 1.82.0) and ‘affyPLM’ (version 1.80.0) packages in RStudio (version 2021.09.20) (Bolstad et al., 2003; Bolstad et al., 2005; Irizarry et al., 2003). Finally, expression matrices were generated, and samples were classified into experimental and control groups for further analysis (Supplementary Table 1).

Each dataset was processed independently to identify differentially expressed genes (DEGs). Two distinct approaches were employed for differential expression analysis between patients and controls: (i) a conventional approach utilizing Welch’s t-test and Empirical Bayes methods, and (ii) an alternative method based on the Comparative Analysis of Shapley value (CASh) technique.

To gain a comprehensive understanding of gene expression patterns, we applied various exploratory techniques commonly used in microarray data analysis. Principal Component Analysis (PCA), heatmaps, and volcano plots were generated to evaluate the distribution of gene expression patterns. PCA was used to illustrate gene expression distribution at two levels: the entire gene set in each dataset versus the differentially expressed genes identified by CASh analysis (p-value < 0.01). Heatmaps were created to visualize DEGs detected through both Empirical Bayes (raw p-value < 0.05) and CASh (p-value < 0.01) analyses, highlighting the clustering of samples according to disease status. Additionally, volcano plots were employed to compare the p-values obtained from Empirical Bayes and CASh analyses, providing a visual representation of the statistical relationships between the methods used for detecting DEGs (see Supplementary Figures 1, 2, and 3 for further details).

2.1.1. Classical Approaches

Conventional analyses for detecting DEGs were conducted using Welch’s t-test, implemented through the ‘multtest’ (2.60.0) package in RStudio (version 2021.09.0) (Pollard et al., 2005). In microarray experiments, the small number of replicates and the large number of genes typically analyzed pose significant challenges, leading to the issue of low statistical power with ordinary t-tests. This limitation makes t-tests less effective for filtering out regulated genes (Åstrand et al., 2008; Gottardo et al., 2003). Moreover, most multiple testing adjustments tend to be quite conservative, particularly when the number of replicates is small (Gottardo et al., 2003). To address this issue, we employed Bayesian-based methods, specifically the Empirical Bayes approach, as implemented in the Bioconductor ‘limma’ R package (bioconductor.org/packages/release/bioc/html/limma.html).

Significant DEGs were identified following multiple testing correction using the Benjamini & Hochberg method to control for False Discovery Rates (FDR) (Benjamini & Hochberg, 1995).

2.1.2. Comparative Analysis of Shapley Value (CASh) Approach

We utilized the Comparative Analysis of Shapley value (CASh) method to identify DEGs by assessing their cooperative contribution to overall gene expression changes (Moretti et al., 2008). The Shapley value, a concept rooted in Game Theory, quantifies the marginal contribution of each gene to the collective expression changes observed within the dataset (Moretti, 2010). CASh combines the Microarray Game algorithm, which is applied to transcriptomic data from microarray experiments, with the Bootstrap technique that involves random resampling of certain values to mitigate the impact of potential outliers in the data matrix (Cesari et al., 2018; Moretti, 2010; Moretti et al., 2008, 2010). In this approach, gene expression is treated as a cooperative game, where each gene collaboratively contributes to the observed expression changes, providing a more nuanced understanding of gene interactions within the dataset.

In this context, a cooperative game is defined by a set $N$ of players (genes) and a characteristic function $v$ that assigns a value to each subset of genes, representing their combined contribution. The Shapley value ${\varphi }_i$ for a gene $i$ is calculated as its average marginal contribution across all possible subsets of genes, providing a robust measure of each gene's importance in the study. In the following lines of this section, we reproduce the mathematical formulation as described elsewhere (Castro-Martínez et al., 2024; Esteban & Wall, 2011).

Formally, given a coalitional game $(N,v)$ for each player $i\in N$ the Shapley value $\varphi \left(v\right)$ is defined by:

$${\varphi }_i={\displaystyle \frac{1}{n!}}{\sum }_{\pi }\left[v\left(P\left(\pi ,i\right)\cup \left\{i\right\}\right)-v\left(P\left(\pi ,i\right)\right)\right],$$

where $\pi$ is a permutation of players, $P(\pi ,i)$ is the set of players that precede player $i$ in the permutation $\pi$ and $n$ is the cardinality of $N$.

We refer to a Boolean matrix (see below) $B\in {\left\{\mathrm{0,1}\right\}}^{n\times k}$, where $k\ge 1$ is the number of arrays, and the Boolean values 0–1 represent two complementary expression properties, for example, normal expression (coded by 0) and over-expression (coded by 1). Let $B.j$ be the $j$-th column of $B$, we define the support of $B.j$, denoted by $sp(B.j)$, as the set $sp\left(B.j\right)=\left\{i\in \left\{1,\dots ,n\right\}:B_{ij}=1\right\}.$The microarray game corresponding to $B$ is defined as the coalitional game $(N,w)$, where $w:2^N\to R^{+}$ is such that $w\left(T\right)$ is the rate of occurrences of coalition $T$ as a winning coalition, i.e., as a superset of the supports in the Boolean matrix $B$; in formula, $w\left(T\right)$, for each $T\in 2^N\{\varnothing \}$, is defined as the value

$$w\left(T\right)={\displaystyle \frac{c\left(\Theta \left(T\right)\right)}{k}},$$

(1)

where $c\left(\theta \left(T\right)\right)$ is the cardinality of the set $\theta \left(T\right)=\left\{j\in K:sp\left(B.j\right)\subseteq T,sp\left(B.j\right)\ne \varnothing \right\},$ with the set of arrays $K=\{1,\dots ,k\}$ and $v(\varnothing )\ne 0.$ Since it is computationally too expensive to calculate the Shapley value $\varphi \left(w\right)$ of game $(N,w)$ according to relation ${\varphi }_i$, Moretti et al. (2008) introduced an easy way to calculate $\varphi \left(w\right)$ for whatever microarray game $\left(N,w\right).$ We have adapted the scripts (available under request) from these authors run under R (r-project.org/).

In our study, the CASh method was applied to detect DEGs using two levels of stringency by setting cutoff p-values at 0.01 (more restrictive) and 0.05 (less restrictive). These genes were then analyzed to distinguish between over-regulated and under-regulated expression levels, based on standard deviations from the control group. Boolean matrices were constructed to represent these expression states, which were subsequently used to define microarray games and calculate the corresponding Shapley values.

A final matrix was generated from the original data, incorporating the expression levels of a selected number of genes and samples. This matrix included genes with raw p-values below 0.01 or 0.05 and categorized the samples into distinct groups (e.g., patients with specific conditions and healthy controls). To identify over-regulated gene expression levels relative to controls, each continuous value in the gene expression vector was coded as 1 if it was equal to or greater than the mean plus the standard deviation of the control group expressions, and as 0 otherwise. This process resulted in a Boolean matrix with values {0, 1} reflecting these conditions.

A similar approach was employed to identify under-regulated gene expressions, where each value less than the mean minus the standard deviation of control expressions was coded as 1, with all other values coded as 0. This produced another Boolean matrix with rows corresponding to genes and columns to samples. These Boolean matrices were then divided according to sample group distinctions, forming separate matrices for each group. From these matrices, microarray games were defined for each condition, and Shapley values were calculated to assess the significance of each gene’s contribution to the conditions under study.

To mitigate the influence of random high Shapley values, a resampling procedure was applied, similar to that described by Moretti et al. (2008). Bootstrap resampling with 1,000 iterations was performed for each analysis. This method, known as Comparative Analysis of Shapley value (CASh), refines the selection of genes significantly associated with the conditions under investigation.

To further reduce the likelihood of false positives, corrections for multiple testing were applied, and Shapley values were compared against statistically significant thresholds. Additionally, Fold Changes (FC) were evaluated, with genes exhibiting p-values below 0.01 and 0.05 and |FC| > 2 considered statistically significant.

2.2. Gene Set Enrichment Analysis and Functional Annotation

In this study, we utilized the g:Profiler functional profiling tool (version e111_eg58_p18_30541362), specifically the g:GOSt module (https://biit.cs.ut.ee/gprofiler/gost), to perform an extensive analysis of the biological processes and pathways influenced by differentially expressed genes (DEGs). This tool leverages Gene Ontology (GO) terms to create a comprehensive, annotated landscape of gene functions and interactions (Kolberg et al., 2023; Raudvere et al., 2019). Gene Ontology provides a structured vocabulary that classifies and integrates biological data across species into three main categories: biological processes, cellular components, and molecular functions. By inputting the list of DEGs into g:GOSt, the tool maps these genes to known GO terms, allowing us to identify which biological pathways and processes are enriched in these genes. This enrichment analysis is instrumental in understanding the roles these DEGs may play in the specific conditions under study. The g:GOSt module performs its analysis by comparing the input gene list against databases of known gene and protein functions, identifying statistically significant over-representations of specific functions or pathways. This is achieved using various statistical methods, including Fisher’s Exact Test, to calculate enrichment p-values, which help determine which processes or pathways are more involved with the DEGs than expected by chance. The results from g:GOSt not only highlight the predominant biological themes associated with the DEGs but also provide insights into the potential molecular mechanisms driving the disease or condition. For example, if a significant number of DEGs are involved in inflammatory response pathways, it may indicate that inflammation plays a crucial role in the pathology of the condition being studied. Moreover, the outcomes from such analyses can guide experimental design by identifying key pathways that could be targeted for further experimental validation or therapeutic intervention. Tools like g:Profiler are indispensable in the genomic era, enabling researchers to translate large datasets of gene expression information into actionable biological insights (Kolberg et al., 2023; Raudvere et al., 2019).

During the analysis of gene expression data, it is crucial to ensure that transcript identifiers (IDs) are accurately annotated and standardized to official gene symbols. This step is fundamental for consolidating data from different sources and facilitating meaningful biological interpretation. To achieve this, we used the g:Convert tool available on the g:Profiler web server (https://biit.cs.ut.ee/gprofiler/convert). This tool converts various biological identifiers into recognized gene symbols, enhancing the consistency and reliability of genomic data analysis. The g:Convert module supports a wide range of biological identifiers, including Ensembl IDs, UniProt IDs, RefSeq, and others, allowing researchers to input data from various experimental outputs and databases. Upon inputting transcript IDs into g:Convert, the tool maps these IDs to the official gene symbols based on the most up-to-date and comprehensive databases. This ensures that subsequent analyses, such as gene expression profiling or functional enrichment, are conducted on verified and universally accepted nomenclature, thereby reducing the risk of errors and inconsistencies. In cases where transcript names are ambiguous, which can occur due to multiple identifiers for a single gene or updates in genomic databases, we prioritized IDs with the most extensive Gene Ontology (GO) annotations. This approach is based on the rationale that identifiers with more comprehensive annotations are likely more researched and documented, offering a higher degree of reliability. By selecting IDs with the most GO annotations, we aimed to enhance the robustness of our dataset, ensuring that the functional analysis reflects well-supported gene functions and interactions. The use of the g:Convert tool in this manner not only streamlines the process of gene annotation but also significantly enhances the quality of the data being analyzed. Accurate annotation is critical, as it directly impacts the interpretation of biological data and the conclusions drawn from research studies. As gene databases are continually updated and refined, tools like g:Convert are invaluable for maintaining the accuracy and relevance of genomic research, providing researchers with confidence in their analytical outputs (Kolberg et al., 2023; Raudvere et al., 2019).

In addition, we conducted an in-depth analysis of Gene Ontology (GO) categories, which categorize gene products based on their involvement in biological processes (BP), cellular components (CC), and molecular functions (MF). To determine the significance of these GO categories, we applied a rigorous statistical criterion, the Benjamini-Hochberg False Discovery Rate (FDR), requiring that GO terms exhibit an FDR value below 0.05 to be considered significantly enriched. This stringent threshold ensures that only the most robust associations are identified, minimizing the likelihood of false positives. We then identified the top ten significantly enriched GO terms within each category for further investigation. To visually represent the findings, the top ten significantly enriched GO terms in each category were plotted for CASh 0.05 comparisons using the ‘ggplot2’ (version 3.5.1) package in RStudio (Wilkinson, 2011).

3. Results

3.1. Datasets and Samples Analyzed

Gene expression data from nine datasets, encompassing a total of 506 samples, were accessed for this study. Table 1 provides an overview of the main characteristics of these datasets.

To detect differentially expressed genes (DEGs), we employed two distinct strategies. First, conventional methods based on Welch’s t-test and Empirical Bayes were applied. Following this, we conducted an alternative analysis using the CASh method. The conventional approach, utilizing Welch’s t-test and Empirical Bayes, failed to identify any DEGs. In contrast, the CASh method successfully revealed several transcripts when using both 0.01 and 0.05 cutoff p-values for the preselection of DEGs (Table 2). Complete lists of DEGs detected for each dataset after these comparisons are provided in Supplementary Table 2. Our analyses demonstrate that the application of the CASh method significantly improves the detection of differentially expressed genes across the nine datasets analyzed.

3.2. Functional Enrichment Analysis of the Differentially Expressed Genes

Due to the restrictive criteria applied when using the CASh 0.01 method and the FDR correction of p-values, the number of DEGs detected was insufficient to identify significantly enriched pathways associated with some gene sets (data not shown). However, functional enrichment analysis using the DEGs identified with the CASh 0.05 method revealed several significantly enriched processes in the analyzed datasets.

In the Autism Spectrum Disorder datasets (GSE6575, GSE18123, and GSE25507), DEGs were primarily associated with biological processes (BP) such as structure development and cardiac muscle cell development, while cellular components (CC) were mainly related to intracellular organelles (Figure 1).

For the schizophrenia datasets (GSE17612, GSE62333, GSE12654, and GSE53987), the most significantly enriched BPs were related to the regulation of programmed cell death, regulation of primary metabolic processes, and multicellular organism development. The CC results highlighted cytoplasm, nucleoplasm, and extracellular space, while molecular function (MF) analysis identified activities mainly associated with glycine-tRNA ligase activity (Figure 2).

In the bipolar disorder datasets (GSE5389, GSE7036, GSE12654, and GSE53987), gene set enrichment analysis revealed telomeric and metabolic processes as a significantly enriched BP. Additionally, synapse, phosphatidylcholine-sterol O-acyltransferase activator activity and chromosome, and protein binding were identified as significantly enriched CC and MF, respectively (Figure 3).

Finally, in the major depressive disorder datasets (GSE12654 and GSE53987), gene set enrichment analysis highlighted epiboly and wound healing spreading of cells as significantly enriched BPs, while nucleoplasm, site of polarized growth and growth cone were identified as significantly enriched CCs and protein binding as significantly enriched MF(Figure 4).

4. Discussion

Neurological pathologies inflict significant suffering and pose substantial burdens on millions of people worldwide. In recent years, the advent of omics technologies has enabled a comprehensive exploration of the molecular patterns associated with many of the most common neurological conditions. Microarray technology, which emerged nearly three decades ago with the goal of studying whole gene expression profiles, has since provided unprecedented insights into the dysregulated molecular pathways involved in disease pathogenesis (Ducray et al., 2007; Shai, 2006). In the present study, we analyzed data from nine datasets generated using Affymetrix microarray devices, including three datasets from Autism Spectrum Disorder, two from Schizophrenia, two from Bipolar Disorder, and two datasets encompassing samples from Schizophrenia, Bipolar Disorder, and Major Depressive Disorder.

Raw data were downloaded from the GEO public repository, and gene expression files were pre-processed, quality controlled, and normalized. To detect differentially expressed genes (DEGs), we employed two strategies: (i) a traditional approach using classical statistical t-tests and (ii) an alternative approach utilizing the CASh method (Moretti et al., 2008). The traditional t-test approach identified few DEGs, whereas the CASh method revealed a significant number of statistically relevant genes across the nine datasets analyzed. The t-test identifies genes based on their differential expression between two conditions, considering a gene significant when its p-value falls below a pre-established threshold (0.05 adjusted p-value in our study). In contrast, the CASh method not only considers the expression of each gene under two conditions but also evaluates the contribution of each gene across all possible permutations using the Shapley value as a measure. This holistic approach mitigates the impact of confounding variables by considering the overall gene network rather than isolated gene expressions. However, a current limitation of the CASh method is that it does not explicitly account for potential confounding effects, which should be addressed in future applications (Cesari et al., 2018; Moretti et al., 2008, 2010; Sun et al., 2020). In summary, CASh offers a more nuanced understanding of gene interactions and their collective impact on disease pathophysiology.

Interestingly, the functional enrichment analysis of the DEGs detected using the CASh method, confirmed previous findings on the molecular bases of the neurological pathologies studied. For instance, processes related to cardiac muscle cell development in ASD samples are directly linked to vascular abnormalities observed in patients with this phenotype (Yao et al., 2006). In the case of Schizophrenia, the regulation of primary metabolic processes and glycine-tRNA ligase activity emerged as significant processes, which are particularly relevant given the metabolic issues associated with Schizophrenia (Von Hausswolff-Juhlin et al., 2009). Similarly, Bipolar Disorder was linked to several key findings in our study, including the positive regulation of lipoprotein lipase activity and synapse and phosphatidylcholine-sterol O-acyltransferase activator activity, which align with the known association of this disorder with altered fatty acids (Saunders et al., 2016). For Major Depressive Disorder, characterized by inflammation and neurological damage, we identified processes such as "wound healing spreading of cells" and "growth cone" as significant in the context of differential gene expression.

5. Conclusions and Limitations

This study highlights the power of Comparative Analysis of Shapley values (CASh) in revealing complex genetic insights into neurological disorders such as Autism Spectrum Disorder (ASD), Schizophrenia, Bipolar Disorder, and Major Depressive Disorder. CASh has been proven as highly effective in identifying differentially expressed genes, many of which are missed by traditional statistical methods, offering a more nuanced understanding of the molecular mechanisms underlying these conditions. These findings open new opportunities for developing innovative diagnostic and therapeutic strategies that may shed light in the etiology of these complex conditions.

However, several limitations should be considered. The inherent complexity of microarray data—such as noise, batch effects, and variability in sample quality—can introduce biases that affect the accuracy of gene expression analysis, despite the rigorous preprocessing and normalization applied. Additionally, the reliance on public datasets may bring biases related to differences in data collection methods, patient selection, and experimental design, potentially limiting the generalizability of our results. To mitigate these issues, future studies should be conducted to validate the findings by using more diverse cohorts of patients.

Looking ahead, integrating CASh with complementary omics technologies, such as proteomics and metabolomics, promises a more comprehensive view of the pathophysiological processes in brain diseases. This combined approach could significantly improve the development of multi-marker panels, enhancing diagnostic accuracy. Longitudinal studies using CASh could also track disease progression and treatment responses, providing insights into how gene expression evolves over time in relation to disease states.

A further challenge is the computational intensity of CASh, particularly with large datasets. The method requires substantial computational resources, and interpreting Shapley values can be complex. Simplifying the approach—through algorithm optimization or data reduction—would make CASh more accessible for routine clinical and research applications. Additionally, CASh does not account for post-transcriptional modifications or protein-level interactions, which are critical for a complete understanding of disease mechanisms. Future work could address this by integrating CASh with proteomic and metabolomic data to offer deeper insights at the protein level.

Another key limitation is the lack of experimental validation of the identified differentially expressed genes. To confirm the biological relevance of these findings, future studies should incorporate in vitro functional assays, such as gene knockdown or overexpression experiments. Moreover, in vivo studies in animal models would help to further elucidate the roles of these genes in disease mechanisms and assess their potential as therapeutic targets.

Achieving the full potential of CASh will require strong interdisciplinary collaboration. Geneticists, neurologists, oncologists, and bioinformaticians must work together to conduct large-scale studies that validate and refine the gene signatures identified, translating these discoveries into practical clinical applications. By advancing our understanding of the genetic basis of neurological disorders, this research contributes to precision medicine approaches, ultimately improving patient outcomes and reducing the global burden of these conditions.