Exploring Promising Biomarkers for Alzheimer’s Disease through the Computational Analysis of Peripheral Blood Single-Cell RNA Sequencing Data

Abstract

Alzheimer’s disease (AD) represents one of the most important healthcare challenges of the current century, characterized as an expanding, “silent pandemic”. Recent studies suggest that the peripheral immune system may participate in AD development; however, the molecular components of these cells in AD remain poorly understood. Although single-cell RNA sequencing (scRNA-seq) offers a sufficient exploration of various biological processes at the cellular level, the number of existing works is limited, and no comprehensive machine learning (ML) analysis has yet been conducted to identify effective biomarkers in AD. Herein, we introduced a computational workflow using both deep learning and ML processes examining scRNA-seq data obtained from the peripheral blood of both Alzheimer’s disease patients with an amyloid-positive status and healthy controls with an amyloid-negative status, totaling 36,849 cells. The output of our pipeline contained transcripts ranked by their level of significance, which could serve as reliable genetic signatures of AD pathophysiology. The comprehensive functional analysis of the most dominant genes in terms of biological relevance to AD demonstrates that the proposed methodology has great potential for discovering blood-based fingerprints of the disease. Furthermore, the present approach paves the way for the application of ML techniques to scRNA-seq data from complex disorders, providing new challenges to identify key biological processes from a molecular perspective.

1. Introduction

Alzheimer’s (AD) disease is a chronic degenerative disease of the brain and is the most common cause of dementia. Its main clinical manifestation is memory impairment, but other brain functions are also affected, and thus, the patient’s social life is seriously affected [1]. Exploring the molecular characteristics of blood cells using single-cell RNA-sequencing (scRNA-seq), data have emerged as a powerful approach to unravelling the molecular basis of AD at the cellular level [2]. It has enabled the identification of transcriptional changes in individual cells, providing unprecedented resolution and sensitivity to detect early changes in disease progression. Furthermore, scRNA-seq can capture the heterogeneity of cell types within the peripheral immune system, which is crucial for understanding the immune response to AD [3,4]. Due to these advantages, scRNA-seq analysis of liquid biopsies, such as blood, is considered a first choice to address the challenges of AD discovery as well as non-invasive blood-based sampling [5]. Although this area may offer promising results, it is in its infancy, and limited studies have been published [6,7,8]. It is worth mentioning that machine learning (ML) analysis has not yet been applied to scRNA-seq data for the purpose of extracting valuable knowledge and identifying potential AD biomarkers. Therefore, there is a significant need for novel approaches to extract insights from such data and advance our understanding of AD.

One major challenge in scRNA-seq data lies in the fact that these are often generated in multiple batches, which can introduce batch effects that confound downstream analyses [9]. In order to address this issue, deep learning (DL) techniques have been developed to integrate scRNA-seq data from multiple batches [10,11,12]. These approaches can effectively remove batch effects, producing a more accurate representation of underlying biological processes. In recent years, there has been a significant advance in the application of deep learning methods for the integration of single-cell RNA sequencing (scRNA-seq) data [13,14]. These methods have shown remarkable progress in addressing the challenges posed by batch effects in scRNA-seq datasets, enabling the integration of data from multiple batches and producing a more accurate representation of underlying biological processes.

During single-cell RNA-sequencing studies, the initial FASTQ file output is transformed into a count’s matrix that provides a summary of the number of molecules for each gene that can be detected in each cell in the dataset [15]. After going through preprocessing procedures, the resulting counts matrix contains tens of thousands of gene expression profiles (features) for each sample, leading to an extremely high-dimensional dataset. Two approaches are commonly employed to handle high-dimensional datasets: dimensionality reduction methods and feature selection [16]. The first entails transforming the initial high-dimensional space into a lower-dimensional space while maintaining pairwise distances between points below a specified threshold. The latter aims to identify the most relevant features that encapsulate relevant information in the dataset. Concentrating on feature selection preserves the gene structure, thus lending interpretability to the studied model.

In this context, we selected the study conducted by Xu et al. (2021) [6] as it provides the most indicative scRNA-seq approach for peripheral five subtypes of immune cells in AD, offering gene expression profiles for tens of thousands of selected cellular samples. The authors employed a variety of analytical processes, including statistical tests, network-based analysis and functional enrichment analysis, and promising results were obtained, revealing abnormal changes in the immune cell composition and transcriptional state of AD. These findings suggest a disturbance in the immune infiltration of the peripheral immune environment at a single-cell level. While statistical approaches and enrichment analysis offer reliable outcomes, the use of ML methods can provide an alternative perspective and uncover hidden insights that may be missed by traditional methods. In addition, when processing scRNA-seq datasets from multiple groups, such as AD patients with different experimental batches in our case, batch effect issues need to be addressed [17]. Herein, we propose an ensemble pipeline that can extract the most significant gene signatures from scRNA-seq data for a case under study. Our method comprised two main steps: first, we employed a cutting-edge deep-learning-based method for scRNA-seq data integration, which accurately removes batch effects; second, we employed an ensemble methodology that utilized three ML-based feature selection algorithms and a voting scheme for ranking the leading genes. Further biological interpretation using enrichment analysis stressed our findings and extracted valuable biological insights relating to the functions of a group of genes.

2. Materials and Methods

2.1. Dataset

To conduct our analysis, we utilized the peripheral blood cell dataset generated by Xu et al. (GEO accession number: GSE181279) [6]. This dataset includes 36,849 peripheral blood mononuclear cells (PBMCs) obtained from three AD patients with an amyloid-positive status and two cognitively normal controls with amyloid-negative status. Our objective was to identify the most dominant genes that differentiated AD patients from healthy individuals.

The major challenge here arose from the nature of single-cell technologies, which create batch effects. This is a systematic technical variation that occurs due to experimental factors such as differences in sample preparation, processing, and sequencing. These effects can lead to false conclusions and mask true biological differences between the cells. Thus, it is important to address them during data integration in order to ensure accurate and reliable results.

2.2. Methodology

Our pipeline consisted of two primary steps: (a) the integration of scRNA-seq data by eliminating batch effects using a deep-learning approach and (b) utilizing an ensemble methodology that combined three ML-based feature selection algorithms and a voting scheme to identify the most dominant genes across single-cell RNA-seq data (Figure 1). Further details for each step are provided below. The data were pre-processed and visualized through the use of Python’s Scanpy library [18].

2.2.1. Deep Learning for Data Integration

The observed raw transcriptomic dataset was divided into five batches, with each batch labeled to indicate the source subject of the cells collected. To create a unified dataset from the five batches, we utilized a state-of-the-art tool for integrating single-cell data. We applied SCALEX [10], a deep-learning technique that integrates single-cell data by projecting cells into a batch-invariant, shared cell-embedding space without the need for model retraining. The core of the SCALEX approach is a variational autoencoder (VAE) that maps distinct dataset batches into a common low-dimensional space free from batch effects. To achieve this, the authors trained a batch-free encoder and a batch-specific decoder within the VAE framework. More specifically, the encoder component of SCALEX is batch-free and extracts biologically relevant latent features, while the decoder is batch-specific and uses batch information and a domain-specific batch normalization layer (DSBN) to reconstruct the original data from the latent features. The encoder maps individual cells to a common space using a generalized projection function that enables the integration of new data without requiring model retraining. Furthermore, SCALEX employs a mini-batch approach where data are sampled from all batches and subjected to batch normalization to adjust for deviation and align with the overall input distribution.

To integrate the data, we first classified the five batches based on the health status of individuals from whom the cells were collected, categorized as a “normal control” “AD” (indicating Alzheimer’s disease). We then merged the batches into a single dataset and employed SCALEX to adjust the batch effects. We used the default parameters for this analysis. After preprocessing and integration/correction, the dataset comprised 35,685 (22,775 AD and 14,074 NC) cells and 2000 genes. Figure 2 shows four data tSNE plots depicting the data before and after integration/batch effect correction with SCALEX, with color-coded batch names and disease types. In addition, we utilized three other cutting-edge techniques for integrating single-cell data, underscoring the capacity of SCALEX (see Supplementary Materials and Figures S1–S3).

2.2.2. An Ensemble Machine Learning Framework for Dominant Genes Identification

Following the first step with the integration of the data batches, we obtained a dataset containing FS. Although the dataset was not considered high dimensional, the large number of genes posed a challenge in conducting a rigorous and reliable analysis. To overcome this challenge, feature selection was employed as a technique to address the issue.

Three different feature selection methods were utilized, each with a distinct purpose, and all of them shared a core perspective from a machine learning standpoint. The aim here was to exploit the results of three well-established feature selection methodologies, each designed to identify relevant genes in a distinct manner. More specifically, the first method leveraged differential expression analysis using a logistic regression approach to identify genes that exhibited significant differences in their expression levels across distinct groups [19]. The second approach (called triku) was based on the nearest neighbor search, aiming to identify genes closely associated with the groups of interest [20]. By utilizing the variable importance measure from the random forest algorithm [21], the third approach leveraged the inherent potential of tree-based methods to examine the significance of each gene in a machine learning prediction process. This technique demonstrated effectiveness in identifying relevant features across various datasets, including single-cell RNA-seq data [22]. In our analysis, we utilized the Mean Decrease Impurity (MDI), which measures the reduction in the impurity of the node caused by a variable. The importance of a variable can be calculated as the total reduction in impurity across all trees in which the variable was used.

In the next step, a set of 500 top genes was extracted from each feature selection method. These three lists were merged using a Borda-based aggregation technique [23] to produce a “consensus” list that included all the features identified as significant by the three feature selection methods utilized in this study. The Borda count score is a voting system that assigns a score to each candidate based on their ranking in a list of preferences. The candidate with the highest overall score is considered the winner. The score is calculated by assigning points to each candidate based on their position in the ranking, with the top candidate receiving the highest number of points and subsequent candidates receiving decreasing points. The total score for each candidate is then calculated, and the candidate with the highest score is considered the winner. The Borda count method is commonly used in elections and decision-making processes to determine the most popular or preferred option.

3. Results and Discussion

To achieve a more robust outcome, we utilized three distinct machine learning-based feature selection methods in this study, each with its unique procedure and significance. Regarding the first method, the objective of conducting differential expression analysis was to discern the genes that exhibited differential expression levels across specific conditions. The identification of such genes could provide valuable biological insights into the underlying mechanisms affected by the conditions of interest. Logistic regression is a simple yet powerful algorithm that proves particularly useful for scRNA-seq analyses [19]. This is due to the substantial number of cells typically available in scRNA-Seq experiments, as well as the ability to incorporate transcript information in gene-level testing. By utilizing logistic regression, the contribution of individual transcripts (isoforms) to the gene-level differential analysis could be identified, which enhanced the interpretability of results. Rather than adopting the conventional method of using cell labels to study gene expression, logistic regression was utilized for each gene to predict cell labels using constituent transcript quantification. This approach capitalizes on the ability of logistic regression to highlight the optimal linear combination of the number of isoforms for differential analysis. In contrast to other methods, fitting a linear logistic regression was relatively and computationally efficient, which was particularly important as the scale of scRNA-Seq experiments continued to grow.

Triku is a feature selection algorithm that belongs to the second category. It selects genes that have a localized expression and are expressed in neighboring cells that share similar characteristics (i.e., functions). The idea behind triku is that biologically important genes have localized expressions in groups of similar cells. It distinguishes between genes that are expressed in a few unrelated cells (uninteresting) and genes expressed in a set of similar cells (interesting) by looking at the expression of nearby cells. Triku calculated the distribution of gene expression values considering the expression of neighboring cells, and this helped distinguish between the two cases. Finally, the expected distribution of gene expression values was calculated using a random sample of k cells.

Variable importance (VI) in random forests is a well-known feature selection method that estimates the importance of each feature in a dataset by measuring how much each feature contributes to the accuracy of the model [21]. In the random forest algorithm, multiple decision trees are built on different subsets of the data, and the VI score of each feature can be calculated based on how much the accuracy of the model decreases when that feature is randomly permuted. Features with high VI scores are considered to be more important, as they have a larger impact on the accuracy of the model.

Regarding our outcomes, the tSNE plots clearly indicated a significant change in the scatter pattern of the data after integration with SCALEX, suggesting that the use of SCALEX effectively corrected the batch effects, making the data appropriate for subsequent analyses. Our study capitalized on the results of three well-established feature selection methodologies, each designed to identify relevant genes in a distinct manner. By exploiting the heterogeneity across these methods, we were able to highlight a diverse set of important genes, offering a more comprehensive understanding of the underlying biological mechanisms at play in the analyzed dataset. Specifically, we utilized a differential expression-based approach, the nearest neighbor search strategy, and an embedded method employing the variable importance measure of a random forest algorithm, each contributing unique insights to the overall analysis.

It is important to point out that our analysis had five distinct single-cell RNA-seq studies (three AD and two normal control). To ensure the reliability of our analysis, we had to tackle batch effects that could have resulted in biased or incorrect results when simply consolidating these studies. Without the proper integration of scRNA-seq data, we may have observed a spurious separation between Alzheimer’s and normal cells, leading to the incorrect identification of dominant genes that may have only been relevant due to the batch effects. However, by performing an accurate deep learning-based integration of the data using the cutting-edge SCALEX tool, we were able to enhance the reliability of our findings and extract valuable biological insights.

In addition, our feature selection approach integrated three different machine learning-based strategies to capture the dynamics of gene selection. This integration resulted in a more robust outcome by identifying genes that were dominant from multiple computational perspectives. Moreover, our ranking score, which was based on the Borda count, provided a fair prioritization of the top genes that could distinguish between Alzheimer’s disease and a healthy state. Our in silico pipeline, combined with a thorough biological interpretation, yielded reliable transcriptional biomarkers for AD. However, further in vitro and in vivo validation of these findings is necessary. We conducted an enrichment analysis using Enrichr [24,25] to explore the high-level functions and utilities of the biological system according to the exported genes. The GO function and Reactome [26] pathway enrichment analyses were performed for most of the provided DEGs (Tables S1–S4). Reactome pathway analysis showed that DEGs were mainly enriched in pathways such as the adaptive and innate immune system, the translocation of the immunological synapse, and PD-1 signaling and neutrophil degranulation were also provided (Figure 3A). The enriched GO terms were divided into the biological process (BP), molecular functions (MF), and cellular component (CC) ontologies. The results of the GO analysis revealed that DEGs were mainly enriched in BP, including the regulation of the immune system, interferone-gamma-mediated signaling, antigen processing, the cellular defense response, and the positive regulation of T cell-mediated immunity (Figure 3B). MF analysis indicated that the DEGs were significantly enriched in MHC class II receptor activity, MHC class II protein complex binding, serine-type peptidase activity, ferric ion binding, superoxide-generating NADPH oxidase activator activity, Tat protein binding and serine-type endopeptidase activity (Figure 3C). For CC, the DEGs were enriched in an MHC and MHC class II protein complex with the luminal side of the endoplasmic reticulum membrane, the integral component of the luminal side of the endoplasmic reticulum membrane and endoplasmatic reticulum for the Golgi transport vesicle membrane (Figure 3D).

More precisely, based on our analysis, EFhd2 encodes a conserved calcium-binding protein that is highly involved in neurodegeneration. It is widely expressed in the brain and central nervous system with cell-type specific functions, including microtubule transport, and is more abundant in cortical neurons [27]. Mielenz and Gunn-Moore found that the down-regulation of the protein could enhance synaptic delivery through kinesin-mediated transport and regulate synaptic plasticity through F-actin binding and Cdc42, Rac and Rho activity [28]. Previous studies identified how EFhd2 interacts with pathologically aggregated tau protein and accomplishes the formation of an amyloid structure [29], while recently, it was reported that EFhd2 in mice interact with several intracellular components such as vesicle trafficking modulators, cellular stress response-regulating proteins, and metabolic proteins [30]. Pathway enrichment analysis showed that a EFhd2 was enriched in the RHO GTPase cycle and signaling by Rho GTases and Miro GTPases (Table S5), which are involved in the regulation of cell migration, neuronal development and mitochondrial distribution with neuronal activity [31,32]. MF analysis showed enrichment in cadherin binding, with these proteins’ participation in the development and functional organization of the adult neural tissue, axon elongation and synaptogenesis [33]. The ribosomal protein L11 (RPL11) was also detected, which encodes RPL11. In a recent study, RPL11 was detected in the brain cerebral gray and white matter of AD and control donors [34]. Regarding TGFBR3, increased levels of this gene have been observed through functional genomic analysis in the temporal cortex of AD patients and particularly, its promoter methylation level was downregulated and negatively related to the Aβ level. The promoter hypomethylation of the gene in Aβ accumulation has been executed through the enhancement of β- and γ-secretase activities, while mediators of the synaptic vesicle cycle, calcium signaling pathway and glutamatergic synapse pathways are associated with an expressed protein [35].

Ras Homolog Family Member C (RhoC), another high-scored gene in our analysis, encodes a member of the Rho family of small GTPases, which are involved in the reorganization of actin cytoskeleton while their expression of some members of this family has been linked with the diminution of human brain malignancy [36]. This was verified through pathways enrichment analysis according to Table S6, which shows mainly Rho GTPases activation pathways as well as semaphorin interactions and sema4D-induced cell migration. In neural stem/progenitor cells, the expression of RhoC was down regulated by the presence of Aβ42, demonstrating its contribution to the migration of these cellular types as a sufficient indicator for NSPCs’ migration capacity [37]. GO analysis demonstrated that RhoC in BP was significantly enriched in the regulation of lipase activity, cortical cytoskeleton organization and mitotic cytokinesis in MF with GTP and guanyl ribonucleotide binding and in CC with neuron projection, cytoskeleton and intracellular membrane and non-membrane-bounded organelles (Tables S7–S9).

LMNA encodes lamin A/C: a two-dimensional matrix of proteins located next to the inner nuclear membrane, which is both made from the LMNA gene through alternative splicing. Deficiencies in A-type lamins account for cardiomyopathy, muscular dystrophy, peripheral neuropathy and progeroid disorders, even if the significance of lamin A in the brain is not clear, probably because A-type lamins are expressed at later stages when cells undergo differentiation. According to previous studies, the expression of lamin A may depend on the epigenetic regulation of LMNA. However, it should be stressed that the adult brain preferentially expresses C-type lamin rather than type-A lamin [38]. Although gene mutations have been related to an increasing number of diseases, often involving skeletal and cardiac muscle, Benedetti et al. indicated a correlation of the gene with patients affected by the neuromuscular disease [39]. Pathways enrichment analysis showed that LMNA was mainly enriched in the depolymerization of the nuclear lamina, deregulated CDK5 neurodegenerative pathways in AD models, and the defective intrinsic pathway for apoptosis and meiotic synapsis, while BP showed enrichment in the regulation of cell aging and mitotic nuclear membrane organization and reassembly (Tables S10 and S11).

Fc Epsilon Receptor Ig (FCER1G), another observed high score gene, was expressed in dissociated mouse dorsal root ganglia neurons, but the function of the neuronal protein was unknown. Liu et al. showed the significant upregulation of the protein in small-sized trigeminal ganglion neurons in the ovalbumin-sensitized mouse model [40]. The results of the Reactome pathway analysis revealed that FCER1G was enriched in platelet adhesion for collagen exposure and dectin-2 family and neutrophil degranulation, while BP showed enrichment in neutrophil degranulation, neutrophil activation and myeloid cell activation when involved in immune response (Tables S12 and S13). Fc Gamma Receptor IIIa (FCGR3A) encodes a receptor for the Fc portion of immunoglobulin G and participates in the removal of antigen–antibody complexes from the circulation, antibody-dependent cellular mediated cytotoxicity, and the antibody-dependent enhancement of virus infections. FCGR3A was recently detected as a significant contributor to peripheral blood nucleated cells derived from AD patients in an integrated pipeline based on single-cell analysis [41]. In BP, it was found to be significantly enriched in peptide cross-linking and the regulation of inflammatory responses.

PELI1 (Pellino E3 Ubiquitin Protein Ligase 1) is involved in the negative regulation of the necroptotic process as well as in protein polyubiquitination and promotes responses to lipopolysaccharide. Its expression can be detected in murine brain neural cells, with the highest expression level being in the microglia [42]. It is also involved in microglial Aβ phagocytosis by inhibiting the expression of scavenger receptor CD36, which has been demonstrated as a potent therapeutic target against neurodegenerative diseases. [43]. According to GO analysis, BP PELI1 enrichment is depicted in the Toll signaling pathway, the regulation of the necroptotic process and protein K63- and K48-linked ubiquitination (Table S14). The important role of E3 ubiquitin ligases and deubiquitinating enzymes in the pathogenesis mechanisms of neurodegenerative diseases have been extensively analyzed, such as their association with a-synuclein in Parkinson’s disease, as well as their deficiency in mitophagy and mitochondrial dynamics in relation to Parkin’s response [44]. ZEB2, which encodes the Zinc Finger E-Box Binding Homeobox 2 protein, significantly enriched the regulation of transforming growth factor beta receptor signaling pathway and regulation of organelle organization (Table S15). Upregulated levels of this transcription factor have been detected in the prospective neural plate and induced the expression of Sox2 and Ncam neural-specific genes while participating in the differentiation of granule neurons in the cerebellum and together with miR-200c as maintain the balance between midbrain dopaminergic progenitor proliferation and neurogenesis [45]. NR4A2 (nuclear receptor subfamily 4 group A member 2) gene encodes a member of the steroid-thyroid hormone-retinoid receptor superfamily, which acts as a nuclear receptor and a transcription factor. Mutations in this gene have been associated with disorders related to dopaminergic dysfunction, including Parkinson’s disease, AD and multiple sclerosis. It has been linked with neuroinflammation, neuronal cell death as well as neuroprotection, and its expression has been reduced in PD post-mortem brains in PD patients. It regulates dopaminergic neuronal differentiation and also localizes in microglia and astrocytes [46]. Nr4a2 is induced by neuronal activity and regulates hippocampal synaptic plasticity and hippocampal-dependent memory [47]. The results of GO analysis indicated that NR4A2 in BP was enriched in midbrain dopaminergic neuron differentiation, dopaminergic neuron differentiation, central nervous system neuron differentiation, regulation of apoptotic signaling pathway and neuron migration (Table S16). The functional role of the integrin alpha4 subunit (ITGA4) in neural crest cell migration has been detected [48] as well as an anti-integrin therapeutic strategy against multiple sclerosis has been described [49] and was found significantly enriched in integrin cell surface interactions in immunoregulatory interactions between a lymphoid and a non-lymphoid cell and extracellular matrix organization pathways (Table S17). It participates also in T cell trafficking during various inflammatory responses and in CNS pathologies such as autoimmune encephalomyelitis, while improved cognitive function in mice has been induced by the inhibition of α4β1 integrin [50]. Further ontological analysis and gene-annotation cluster networks using GeneCodis4 were implemented, as shown in Figure 4 (and Figures S4–S7).

Alzheimer’s disease has a complex etiology involving the interaction of genetic and environmental factors. It can be divided into familial and sporadic. The appearance of the familial form is associated with mutations in APP, PSEN1 and PSEN2. The rare familial form of the disease accounts for 1–2% of all dementia cases and is inherited from generation to generation in an autosomal dominant manner; on the other hand, the sporadic form concerns 98–99% of cases and occurs in people over 65 years of age with the main predisposing genes for Alzheimer’s disease being APOE, CLU, CR1 and PICALM [51,52]. Previous studies using scRNA-seq analysis have been used to explore differential cell subpopulations in the PBMCs of AD patients, and specifically expressed genes in B cells were indicated [8] as well as in the post-mortem brains of AD patient tissues, whereas cross-study large-scale transcriptomic approaches were performed to reveal the role of aberrant gene expression in the disease [53,54]. These methods provide well-defined networks of genes that could be potential targets for early disease intervention.

4. Conclusions

Emerging next-generation sequencing technologies, such as scRNA-seq, enable the whole-genome analysis of an unprecedented number of cells, producing large datasets with ultra-high volume and complexity. However, the analysis of such data poses a significant challenge due to its sheer size and complexity. In order to address these challenges, DL and ML techniques can be implemented efficiently. These techniques are well-suited to handle big data due to their inherent nature and can produce robust and reliable outcomes when applied to large and complex datasets such as those generated by NGS technologies. In our study, we proposed an ensemble pipeline that extracted the most significant gene signatures from the scRNA-seq data of PBMCs, based on a cutting-edge deep-learning-based method for scRNA-seq data integration followed by three ML-based feature selection algorithms and a voting scheme with which to rank the leading genes. Compared to other state-of-the-art integration techniques, the applied pipeline, including SCALEX, demonstrated a significantly superior performance making it a reliable tool for large-scale scRNA-seq data integration and providing significant insights within individual cells, which could serve as potent transcriptional markers for unraveling the heterogeneity and complexity of complex diseases.