Identification of Novel Biomarkers in Huntington’s Disease Based on differential Gene Expression Meta-Analysis and Machine Learning Approach

1. Introduction

Huntington’s disease (HD) is an autosomal degenerative disorder in brain that progresses with age [1,2], resulted from the mutation in the huntingtin (HTT) gene in chromosome 4p16.3, where an expansion of CAG trinucleotide repeats (more than 35 CAG) in exon1 of the HTT gene onsets the disease by expressing a toxic HTT protein having expanded polyglutamine (polyQ) amino acids [1-3]. The symptoms of HD can be characterized by a triad of progressive motor, cognitive, and psychiatric/behavioral disturbances [1,4], and are typically categorized into three stages, which are presymptomatic, prodromal, and manifest [5]. The presymptomatic stage is a prolonged stage during which individuals do not demonstrate any indications or symptoms regarding HD despite having expanded CAG trinucleotide mutations [5]. Perhaps, individuals with HD often maintain a typical lifestyle, similar to those without HD, until the onset of the prodromal stage, where individuals with expanded CAG begin to exhibit unspecified or potential motor impairments upon examination and minor but notable cognitive alterations [5]. In the manifest stage, individuals display the triad symptoms of HD, including progressive motor, cognitive, and psychiatric/behavioral alterations [5].

HD is diagnosed based on confirmed motor, cognitive, and behavioral signs and symptoms [6]. These signs and symptoms, specifically, the motor function of HD, are most widely assessed using the Unified Huntington’s Disease Rating Scale (UHDRS), which is based on the diagnostic confidence level (DCL) [6]. The DCL value varies from 1 to 4, where a DCL value of 4 indicates 99% confidence that the motor impairments are due to HD [6]. However, a DCL score of 4 reflects the advanced stage of the disease, indicating that by that time the disease has already progressed to the advanced stage [6]. Hence, detecting any early abnormalities and diagnosing HD during the prolonged presymptomatic stage utilizing DCL score is insufficient [6], as during this stage, individuals do not demonstrate any of the triad of HD symptoms, which would appear as notable clinical signs after decades [5]. Thus, owing to the prolonged duration of the presymptomatic stage [5], clinical staging standards could fail to identify any subtle changes during this stage, highlighting the need for a robust and efficient way to monitor the disease onset as well as the progression and stage of the disease.

Recent research has focused on developing an efficient HD staging approach that could precisely diagnose HD from birth until the medical diagnosis [6]. In 2022, Tabrizi et al. proposed an alternative HD diagnosis system, the Huntington’s Disease Integrated Staging System (HD-ISS), which harnesses biomarkers, besides HD pathophysiology for HD staging [6]. In brief, the HD-ISS systems consist of four stages, 0 to 3, where individuals with stage 0 carry the HTT mutation, but no detectable clinical symptoms [6]. Individuals in stage 1 possess early pathological alterations, whereas stage 2 represents a clinical phenotype, and functional deterioration is associated with stage 3 [6]. However, owing to the stringent landmark assessment criteria of the HD-ISS system, which categorizes HD stages by assessing a measure in at least 100 cohorts, no HD biomarkers have been identified till now [6]. This makes it indispensable to develop novel biomarkers that may satisfy the HD-ISS system. Moreover, developing a prognostic biomarker offers not only a way of diagnosing HD during the presymptomatic stage but also an approach to monitor the disease progression, allowing individuals to take necessary actions to prevent and/or slow down the manifestation of HD [7]. They can also be utilized to monitor the efficacy of a therapy, select additional therapy, and monitor recurrent diseases [8-10].

One effective method for identifying HD-related gene-expression signatures is high-throughput transcriptome profiling. Over the decades, many RNA microarray studies on HD patients and controls using induced pluripotent stem cell-derived neurones, peripheral blood, and postmortem brain tissue have been published in the Gene Expression Omnibus (GEO). Indeed, some recent studies have already approached integrative bioinformatics analysis on some of those datasets, particularly from the blood [11,12], to identify the potential HD biomarkers using conventional statistical methods. However, identifying biomarkers from all genomic datasets is often unsuitable through conventional statistical methods [13], particularly due to the unsuitability of these methods to elucidate the intricate relationship among gene expressions [13]. Moreover, the large dimensional datasets are complex to visualize using these methods, making it hard to interpret the results [13]. In order to overcome these limitations, scientists have recently utilized an integration approach of machine learning (ML) algorithms with SHapley Additive exPlanations (SHAP) and nomogram analysis, which can systematically prioritize the genes and hence improve the efficiency of identifying a novel prognostic biomarker [14]. Besides, this integrated approach can also predict if the features (genes) positively or negatively correlate with the disease [14].

In this study, we developed a rigorous, multi-phase pipeline to identify reliable biomarkers for HD, integrating meta-analysis with interpretable ML (Figure 1). We first identified a core collection of consistently dysregulated genes by meta-analysis of publicly accessible three postmortem HD-affected brain tissue microarray datasets from the GEO database using a lenient threshold to capture moderate yet consistent expression changes across samples. Subsequently, 19 supervised ML algorithms comprising tree-based classifiers including Decision Tree (DT), ExtraTrees (ET), Random Forest (RF), Gradient boosting (GB), Catboost (CB), ExtremeGradientBoosting (XGB), Adaptive Boosting (AB), Histogram Gradient Boosting (HGB), Light Gradient Boosting Machine (LGBM), as well as other classifiers consisting of Multi-Layer Perceptron (MLP), K-Nearest Neighbors (KNN), Logistic Regression (LR), Stochastic Gradient Descent (SGD), Support Vector Machine (SVM) with different kernels, Naive Bayes (NB), and Bernoulli Naive Bayes (BNB) have been used on the selected gene list to select features that best distinguish HD cases from controls. Finally, by combining the robustness of meta-analyses (stringent criteria) with the power of ML to choose features with high dimensions, we identified two potential biomarkers that obtained excellent predictive performance in independent validation and were strongly associated with HD risk.

2. Materials and Methods

2.1. Gene Expression Datasets

This study used publicly available gene expression profiles from the NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/) using the “GEOquery” package from R [15]. Three datasets (GSE3790_GPL96, GSE3790_GPL97, and GSE26927) were considered and generated from distinct sequencing platforms. The GSE3790 dataset contains 404 human brain tissue samples profiled from two different microarray platforms [16,17]. The first platform, GPL96 ([HG-U133A] Affymetrix Human Genome U133A Array), includes 100 HD samples and 101 healthy controls, whereas the second, GPL97 (Affymetrix Human Genome U133B Array), contains 101 HD and 102 healthy brain tissue samples. Lastly, the GSE26927 dataset, generated using the GPL6255 platform (Illumina HumanRef-8 v2.0 Expression BeadChip), includes 20 samples in total, containing an equal number of HD (n=10) and control (n=10) [18,19]. For external validation, we incorporated the GSE9857 dataset [20], which comprises 18 mouse brain tissue samples profiled using the GPL1261 platform Affymetrix Mouse Genome 430 2.0 Array), including 9 samples from HD models and 9 from control mice. All samples from HD were considered diseased (HD) regardless of their subtypes.

2.2. Meta-Analysis

While the datasets were initially normalized and preprocessed according to standard microarray protocols, we additionally applied log transformation to stabilize the variance and enhance cross-sample comparability. Since the datasets originated from different microarray platforms, batch effect correction was performed based on the Combat algorithm of the “sva” R package [21] to minimize platform-specific biases and ensure uniformity across studies. For meta-analysis, expression matrices were harmonized by mapping probe sets to their corresponding gene symbols; in cases where multiple probes mapped to the same gene, the mean expression value was used to generate a unified gene-level profile. Meta-analysis was conducted using the “MetaVolcanoR” [22] package in R, employing Fisher’s p-value combination method to integrate statistical significance across studies. This approach increases analytical power by aggregating independent P-values, enabling robust identification of DEGs across heterogeneous datasets. Two thresholds were used for gene selection. The genes that fulfilled the lenient criteria of an adjusted p-value < 0.05 and an absolute log₂ fold change greater ≥ 0.85 were chosen for ML modelling. Then we applied a stringent criterion of an adjusted p-value < 0.05 and an absolute log₂ fold change > 1, and the genes that fulfilled the condition were selected and further used to identify key genes from the ML analysis.

2.3. Gene Ontology (GO)

The ClusterProfiler [23] package in R was used to analyze the Gene Ontology (GO) terms for the identified DEGs to gain further insight into the biological significance of the selected genes, identifying key biological processes, molecular functions, and cellular components associated with differential gene expression.

2.4. Machine Learning (ML) Modeling

For the ML modeling, we selected the batch effect-corrected combined GSE expression profile as the dataset for ML analysis. We randomly split the dataset into training and testing sets, with 80% of the samples used for training and 20% reserved for testing to evaluate the model’s performance. We selected 19 ML algorithms, leveraging the strengths of various model families to identify the optimal predictive model, ensuring robustness and enhancing the reliability of biomarker identification. We applied grid search hyperparameter tuning based on 10-fold stratified cross-validation in the train set for each algorithm to explore the best set of hyperparameters (Table S1). We computed and examined the AUC (Area Under the ROC Curve) to evaluate the effectiveness of the ML models, along with additional metrics including accuracy (ACC), F1 score, Matthews correlation coefficient (MCC), specificity (SPE), and sensitivity (SEN). To interpret the models, we employed Shapley Additive exPlanations (SHAP) using the Python package “SHAP” [24], which provides both global and local interpretations, enabling us to understand the direction and magnitude of influence for each gene on the predictive outcome.

2.5. External Validation

For external validation, a nomogram was built using the rms package in R using the training set resulting from the ML modelling, which was then tested on GSE9857. The predictive model was developed utilizing genes that overlapped between the machine learning analysis and the meta-analysis performed under a stringent threshold. This approach enables us to assess the selected gene biomarkers’ clinical significance and generalizability. To further evaluate the discriminatory power of the final selected genes in distinguishing between HD and normal samples, ROC curve analysis was performed using the qROC package in R.

3. Results

3.1. Analysis Results of Differentially Expressed Genes

The batch effect was removed from three datasets (GSE3790_GPL96, GSE3790_GPL97, and GSE26927) using ComBat, and the uniformity of the datasets was assessed by using principal components analysis (PCA). The results are plotted in Figure 2a, showing that the batch effect was successfully removed from the datasets. Following this, we conducted the meta-analysis using the R package MetaVolcanoR, where we set the adjusted p-value < 0.05 and absolute log2 fold change ≥ 0.85 to identify the DEGs to consolidate differential expression findings across the cohorts and identify genes consistently dysregulated in HD. We obtained 41 DEGs, as shown in the volcano plot (Figure 2b), of which 10 were down-regulated and the remaining 31 were up-regulated in HD. Additionally, we created a heatmap using these 41 DEGs, which is presented in Figure 2c, which depicts the comparison in expression profiles of each gene between datasets and between the two groups, Control and HD.

3.2. Gene Enrichment Analysis

Next, we examined Gene Ontology (GO) analysis to reveal the functional relevance of identified genes (41 DEGs) in the cellular organization and confirm whether prioritizing candidate genes for ML modeling demonstrates critical cellular activities relevant to HD pathology. As shown in Figure 3a, most of the genes are identified as components of synapses, membrane transport, and ion channel regulators. The GO terms analysis of biological process (BP) shows that regulating monoactomic ion transmembrane transporter activity, neurotransmitter release, dendritic arborization, and action potential are the critical pathways that are significantly dysregulated in HD. Some notable genes associated with these pathways are represented in Figure 3b.

3.3. Machine Learning (ML) Analysis

The 41 genes that met a lenient criteria of adjusted p-value < 0.05 and absolute log₂ fold change ≥ 0.85 from the meta-analysis were input into the ML models for training. A total of 426 samples, containing 238 normal and 186 HD, were randomly divided into training sets (N = 339) and testing sets (N = 87) (Figure S2). We then assessed the classification performance of these ML models on the testing sets (Figure 4a). We chose AUC as the primary evaluation metric for selecting optimal models because it can provide a balanced measure of classification quality, particularly in class imbalance, by considering true and false positives and negatives. The ET model performed very well in categorizing HD and control, obtaining the highest AUC score of 0.867 (Figure 4c) in training and 0.853 (Figure 4d) in testing. Moreover, the MCC of 0.607 in the train set and 0.615 in the test set demonstrates its outstanding performance with minimal difference between the two sets. Besides ET, two other models, CB and XGB, also performed well in classifying the disease from the control. The CB model achieved an AUC of 0.862 (Figure 4c) on the training set and an MCC of 0.587, while on the testing set, it displayed strong discrimination ability with an AUC of 0.854 (Figure 4d). The XGB model attained an AUC of 0.858 (Figure 4d) and an MCC of 0.583 on the test data, while achieving an AUC of 0.870 (Figure 4c) and an MCC of 0.576 during training. We selected these three models to identify significant biomarkers through SHAP analysis.

3.4. Selection of Characteristic Genes from Machine Learning (ML) and Meta-Analysis

We applied the SHAP technique to interpret the predictions of the selected top three ML models (Figure 5a, b, and c), as SHAP values can demonstrate the positive or negative impact of a feature variable on the outcome prediction. We first selected the top 15 genes from each model based on mean absolute SHAP values and compared them to find consistently prioritized genes by all three methods. Notably, NETO2, PAX6, BCL11B, GFAP, GABRD, and PHACTR1 repeatedly appeared among the top predictors (Figure 5d), exhibiting their robust and consistent association with the disease classification. To further investigate how different expression levels of genes contributed positively or negatively to the probability of HD, we computed a SHAP summary plot for each model, presenting the relationship between the expression levels of genes and the degree of high and low SHAP values in the testing dataset. On the SHAP summary plot, feature variables with higher SHAP values were closely associated with an increased risk of developing HD. In the top-performing model ET (Figure 5a), we found that PAX6 shows high feature values x (red dots) clustered on the left side, where SHAP values are negative, which means that the high expression of PAX6 is associated with a lower predicted probability of HD occurring. Conversely, low expression (blue points) corresponded to higher SHAP values and thus increased the model’s predicted probability of HD occurrence. This is seen by other top-performing models (CB, XGB) (Figure 5b and c). The GFAP gene, which is at the top of the XGB model (Figure 5c), shows that high expression levels (red dots) are mainly found in the negative SHAP region, indicating that increased GFAP expression is associated with reducing the risk of HD occurring. Conversely, the feature BCL11B shows the opposite impact, where high expression value is on the positive side of the SHAP region, leading to a significant increase in expression, significantly increasing the model’s prediction of HD occurrence. Other genes follow similar patterns. The high expression value (red points) of GABRD was correlated with positive SHAP values, signifying that an increase in GABRD expression positively influences the model’s prediction of HD. In contrast, lower GABRD expression (blue spots) correlated with negative SHAP values, indicating a reduced risk. A similar expression pattern was also observed for PHACTR1. Subsequently, we selected the 5 robust genes that were found to be present by both ML and meta-analysis with stringent criteria (Figure 5e).

3.5. Construction of a Nomogram and Validation by External Datasets

Next, we constructed a nomogram using the training dataset of ML to predict the potential risks of developing HD and assess the predictive capacity of five features identified by ML models, including PAX6, GFAP, BCL11B, PHACTR1, and GABRD. As depicted in Figure 6a, the nomogram model assigns total points ranging from 0 to 260, where a point from 148 to 218 corresponds to 0.1 to 0.9 of the risk of developing HD. The calibration curve was constructed using bootstrapping (resampling = 1000). As illustrated in Figure 6b, the calibration closely follows the ideal curve, suggesting the model’s accuracy in predicting risks of developing HD. The calibration curve demonstrated a p-value (S.p) greater than 0.05 and an AUC of 0.938, verifying the model’s accuracy (Figure 6c). We additionally assessed the accuracy of the model using an external validation dataset, GSE9857, and an ROC curve was constructed. As shown in Figure 6d, the ROC curve had an AUC score of 0.938, further confirming our model’s accuracy. Next, we conducted decision curve analysis (DCA) using the training dataset of ML and the external dataset to evaluate the clinical effectiveness of our constructed model again. As shown in Figure 6e, the DCA curve displayed that the nomogram had good net benefits for clinical use.

We used the external dataset to analyze ROC curves of the five genes to further evaluate the discriminative ability to distinguish between HD and control samples. The outcome shown in Figure S3 demonstrates that each of the five genes displays an AUC higher than 0.60 in the external dataset. Among the five genes, the GABRD and PHACTR1 achieved an AUC of 1.00 and 0.925, respectively. This indicates their superior ability to distinguish HD and normal. Due to their superior discriminative performance in the external validation dataset, we selected GABRD and PHACTR1 as the final biomarkers for HD.

4. Discussion

In this study, we have utilized three HD microarray datasets studied postmortem brain tissue to disclose a potential biomarker with prognostic value that can detect disease progression during the prolonged presymptomatic stage and monitor the treatment efficacy. Although efforts from previous studies on biomarker identification mostly rely on analyzing blood transcriptome profiling, the most reliable way to understand the molecular basis of disease processes specific to a particular tissue is by analyzing RNA expression in postmortem tissue, because of its sensitivity to local pathology, superiority to cell type resolution and lack of systemic confounding [25-27].

To find reliable biomarkers for HD, we used a two-step workflow in this study: first, a thorough meta-analysis of HD microarray datasets; and second, supervised ML. In many disease contexts, meta-analysis of microarray data has increased statistical power and reproducibility, providing a more meaningful understanding of the underlying biological processes [28]. Our findings extend this paradigm to the discovery of HD biomarkers. In addition, this approach addresses sample heterogeneity by employing data aggregation, stringent quality control, and statistical modelling that explicitly considers inter-study variations [29-31]. We used two criteria for sorting the candidate DEGs (a relaxed and stringent), where a relaxed criterion (adjusted p-value < 0.05, and absolute log2 fold change ≥ 0.85), which identifies 41 genes subjected to ML modeling. A lowered threshold (relaxed criterion) ensures the model picks up weak but biologically significant signals, keeping a broader range of potentially informative features. These genes may not meet stringent statistical cut-offs but could contribute meaningfully to classification performance. Perhaps, analysis of GO enrichment suggests that the dysregulation of these 41 genes is responsible for alterations in some essential signaling pathways, including synapse function and dendritic arborization, which are known to exacerbate HD pathophysiology [32-36] and contribute to motor and cognitive symptoms. Furthermore, striatal neurones, which are particularly affected in HD, show alterations in action potentials and membrane potentials, including changes in the patterns of action potential firing and a shift towards the depolarisation of resting membrane potentials, which are also associated with motor and cognitive symptoms [37-39].

To refine these candidates, we leveraged the strength of 19 ML classifiers, as not all models are suitable for all data types. For example, training a deep neural network with high-dimensional data of a few hundred patient samples will likely result in a highly overfitted model. Hence, it is necessary to carefully choose the right types of models a priori and not purely rely on brute force computational power [40]. When estimating a model’s generalization performance from observational data, the variability in biomedical datasets is often high, due to both technical and biological sources of variation. To address this challenge, we evaluated each model using cross-validation to avoid hyperparameter-tuning bias and obtain nearly unbiased performance estimates [41]. AUC selected the top three models from these ML models, providing consistent model evaluation. AUC evaluates the model using all thresholds, effectively addressing the effect of changing prevalence (proportion of positive cases). Models’ feature rankings were interpreted using SHAP, a game-theoretically grounded method that provides consistent, locally accurate attributions for any model class [24]. By intersecting the top features across these three models, we distilled a consensus list of six genes, thereby reducing algorithm-specific selection biases and focusing on features with broad predictive relevance. Next, we reinforced reproducibility by re-applying meta-analysis with a stringent fold-change cut-off of 1.0, yielding 13 DEGs. The intersection of these 13 genes with our six consensus ML-selected genes produced a core set of five candidates. Crucially, this set was subjected to an external validation cohort by building a nomogram and ROC curve analysis, in which the final two biomarkers GABRD and PHACTR1, achieved AUCs of 1.00 and 0.98, respectively, confirming their high discriminative power and generalizability. Such performance in an external cohort underscores the robustness of our integrative approach and parallels outcomes seen in other ML biomarker studies with rigorous external testing.

The GABAAR δ subunit, GABRD, was consistently given priority by the ML models, which aligns with the recognised role of GABAergic dysfunction in the development of HD. A previous study identified that GABRD is abundantly expressed in the cerebellum, cortex, and hippocampus, localized in the cytosol and membranes of the neurons [42]. Santhakumar et al observed that animals deficient in tonic GABA currents had decreased survival after an in vitro excitotoxic challenge with quinolinic acid, indicating a neuroprotective function against excitotoxic damage in the adult striatum [43]. More importantly, in experimental settings, GABRD is identified as a blood biomarker for Alzheimer’s disease (AD) detection by untargeted proteomics profiling, where GABRD-positive neuronal extracellular vesicles represent synaptic dysfunction and neuronal death [42].

On the other hand, PHACTR1, which is implicated in dendritic arborization in the gene enrichment analysis, was repeatedly selected, highlighting novel facets of HD biology that warrant further mechanistic study. PHACTR1 encodes protein phosphatase and actin regulator 1 that binds actin and protein phosphatase 1 (PP1) and plays a critical role in the actin cytoskeleton and dynamics. PHACTR1 has been identified in the synaptosome, co-localizing with the synaptic marker synaptophysin, detectable in dendrites at DIV21, and strongly expressed in differentiated neurons within the adult mouse cerebral cortex [44,45]. Dendritic spine formation and synaptic plasticity are processes known to be impaired in HD, and they are supported by actin dynamics in neurons. Although there is a lack of direct mechanistic studies on PHACTR1 in HD, a recent investigation revealed that PHACTR1 is among a limited number of genes that are more prevalent in the striatum and whose dysregulation correlates with the emergence of somatic CAG repeats and neurodegeneration in HD murine models [46].

Although our study, using an integrated bioinformatics and ML model, discloses GABRD and PHACTR1 as novel key biomarkers for HD, it still has certain constraints. Firstly, our work lacked any in vivo and/or in vitro validation experiments; instead, it exploited data from a public database to identify potential biomarkers. Secondly, further research is yet to be completed to understand intricately how these genes contribute to HD progression and corroborate their viability as diagnostic markers. Thirdly, the selected datasets had a limited sample size and were only collected from the GEO database. Finally, the selected datasets were processed using different platforms, and although we removed the batch effects using the “sva” package, they could still contain some degree of noise.

5. Conclusions

Identifying a novel prognostic biomarker is indispensable not only to diagnose HD during the presymptomatic stage but also to monitor the progression of HD, allowing individuals to take necessary actions to prevent and/or slow down the manifestation of HD. Utilizing an integrated bioinformatics and ML approach, this study discloses two genes, GABRD and PHACTR1, as candidate potential biomarkers for HD. The findings of this study can be further studied using in vitro and clinical testing to understand how these genes contribute to HD and to ascertain their efficiency as prognostic biomarkers.