Classification of Alzheimer's disease using robust TabNet neural networks on genetic data

1.   Introduction

Alzheimer's disease (AD) is one common progressive neurodegenerative disease that is accompanied by typical clinical symptoms including memory loss and impairment in daily communications and activities ^[1]. It accounts for approximately 60–80% of all dementia and the etiology of AD is still unknown. The occurrence of AD may be caused by the accumulation of specific proteins in and around neurons ^[2]. And it was thought to be having an estimated 70% of the risk attributable to genetic factors ^[3]. The AD population was expected to be 131.5 million by 2050 worldwide ^[4]. Once it was diagnosed, there were no effective treatments and intervened techniques. Therefore, the early diagnosis of AD has played a great role in delaying the progression of AD ^[5].

In clinical applications, the techniques of AD diagnosis mainly involved neuropsychological scores, cerebrospinal fluid (CSF) testing, neuroimaging examinations, and genetic tests ^[6]. Each approach demonstrated typical capabilities and limitations for AD diagnosis. With wide application in clinical practices, the neuropsychological scores were easily disturbed by subjective and objective factors, such as misrepresentation of family members, cross-cultural differences, etc. ^[7]. While CSF testing was a kind of invasive examination of damaged injury for clinical screening ^[8]. According to the neuroimaging examinations, positron emission tomography (PET) and magnetic resonance imaging (MRI) could present morphological, functional and metabolic imaging features. While the metabolic PET was greatly limited by availability and radiation doses ^[9]. MRI modalities preserved different limitations, such as low sensitivity of structural MRI (sMRI), poor stability of functional MRI (fMRI), and incomplete characterization of brain microstructure for diffusion tensor imaging (DTI) for early AD ^[10–12]. Especially, with the robustness of fast and accurate abilities, genetic testing has become the most potential technique for AD detection ^[13].

Conventional machine learning (ML) was applied to deduce the biological problems in genetic data for the AD intelligence diagnosis ^[14,15]. Ha ^[16] proposed a novel computational framework to predict miRNA-disease associations via matrix factorization with a disease similarity constraint, and got the area under curves (AUCs) of 0.9147 and 0.8905 for global and local leave-one-out cross-validation (LOOCV), respectively; then, a similarity-based matrix factorization framework was presented to identify miRNA-disease associations, yielding better AUCs of 0.9227 and 0.8952 for global and local LOOCV, respectively ^[17]. Once, Oriol et al. ^[18] used the MLs of least absolute shrinkage and selection operator (LASSO), k-nearest neighbor (KNN), support vector machine (SVM) on single nucleotide polymorphism (SNP) data to classify late-onset AD, and the SVM achieved the best performance with the AUC of 0.72. Xu et al. ^[19] attempted the application of SVM to analyze the protein sequences encoded by genes for AD classification with an accuracy (ACC) of 85.70%. Again, Castillo et al. ^[20] implemented the SVM with three genes (PSEN1, PSEN2 and APP) to classify AD and obtained the ACC around 80%. Moreover, Voyle et al. ^[21] used recursive feature elimination-random forest (RFE-RF) to distinguish normal controls (NC), mild cognitive impairment (MCI) and AD with an ACC of 62.70% for gene expression data. Besides, Moradi et al. ^[22] applied linear discriminant analysis (LDA) to analyze blood gene expression profiles and got an AUC of 0.84 for the classification of NC and AD, and 0.80 for the classification of NC and MCI. Up to now, classical ML algorithms including SVM, LDA and RF have been widely used and hold big promise for the early detection of AD ^[23]. However, the tedious feature selection and limited performance of MLs still could not meet the clinical needs ^[24].

Recently, deep learning (DL) has been used for the intellectual diagnosis of AD ^[25]. Wang et al. ^[26] proposed an ensemble of 3D densely connected convolutional networks (3D-DenseNets) for AD and MCI diagnosis based on 3D MRI images, with the ACCs of 98.75% for AD/NC, 93.55% for AD/MCI, 98.36% for MCI/NC and 97.19% for AD/MCI/NC, respectively. Yu et al. ^[27] proposed a novel tensorizing generative adversarial network (GAN) with high-order pooling to assess MCI/AD based on the T1-weighted MR images, with the ACC of 89.29% for MCI/NC, 85.71% for AD/MCI, and 95.92% for AD/NC; then, a novel multidirectional perception generative adversarial network (MP-GAN) was presented to visualize the morphological features indicating the severity of AD for patients of different stages, and the MP-GAN achieved better classification performance than GAN in terms of AUC, ACC, specificity, and sensitivity ^[28]. Lee et al. ^[29] used deep neural networks (DNN) to distinguish NC and AD based on blood gene expression data, with an AUC of 0.86. Mahendran et al. ^[30] proposed an enhanced deep recurrent neural network (EDRNN) for DNA methylation data to classify AD and NC, and obtained an ACC of 89.40%. Park ^[31] proposed a deep neural network-based prediction model to assess AD and NC based on gene expression and DNA methylation data, with an average ACC of 82.30%. Due to the low ACC of existing DLs, the betterment of DL models is imperative to improve ACC for AD diagnosis ^[32,33].

In our study, an improved DL approach, namely differential genes screening-TabNet (DGS-TabNet), was proposed for AD binary and multi-class classifications. It was characterized by two modules of DGS and TabNet ^[34,35]. With the usage of DGS, the susceptible genetic features were effectively chosen to facilitate the learning rate and enhance the generalization ability for the proposed model. And the main contributions of our study are summarized as follows:

1) Via the comparison with three novel DLs of multi-layer perceptron (MLP), neural oblivious decision ensembles (NODE), TabNet as well as five classical MLs of decision trees (DT), RF, gradient boosting decision tree (GBDT), light gradient boosting machine (Light GBM) and SVM ^[36–43], our DGS-TabNet achieved the best performance with the ACC of 93.80% for binary classification, and with the ACC of 88.27% for multi-class classification.

2) The biological interpretability related to AD pathogenesis for outputted global important genetic features was validated by gene enrichment analyses of Kyoto encyclopedia of genes and genomes (KEGG) and gene ontology (GO) ^[44]. This would help to acknowledge the capability of genetic features in classification models and provide critical clues for the exploration of the biological interpretability of AD-susceptible genes.

2.   Materials and methods

2.1.   Subjects

Two batches of gene expression data were collected from the gene expression omnibus (GEO) datasets (http://www.ncbi.nlm.nih.gov/geo) from the AddNeuroMed consortium. GEO is an open international public repository that archives and freely distributes microarray, next-generation sequencing, and other forms of high-throughput functional genomics data. It also offers a number of web-based tools and strategies to assist users to query, analyze and visualize data. During the genetic testing, batch 1 (GSE63060) was prepared by Illumina Human HT-12 v3 Expression BeadChips, and batch 2 (GSE63061) was processed by Illumina Human HT-12 v4 Expression BeadChips.

This study enrolled 145 patients with AD, 80 MCIs and 104 NCs for batch 1,131 patients with AD, 101 MCIs and 126 NCs for batch 2. The detailed demographics of enrolled participants were listed in Table 1. The enrolled standard of the AD group met with the requirement of a diagnostic and statistical manual of mental disorders (DSM) ^[45]; MCI participants were included according to Petersen diagnostic criteria; NC cases were chosen by complying with the consortium to establish a registry for Alzheimer's disease (CERAD) ^[46].

Table 1.  Demographic information of enrolled participants.

Demographic information	AD		MCI		NC
	Batch 1	Batch 2	Batch 1	Batch 2	Batch 1	Batch 2
Number of subjects	145	131	80	101	104	126
Age (Mean ± SD)	75.40 ± 6.58	77.82 ± 6.64	74.45 ± 5.99	78.06 ± 7.40	72.37 ± 6.34	75.35 ± 6.04
Sex (Male/Female)	46/99	50/81	41/39	43/58	42/62	48/78

---

 | Show Table

DownLoad: CSV

2.2.   Preprocessing

The preprocessing of gene expression data was performed by R (version 4.0.2), with the following steps: 1) The GSE63060 and GSE63061 datasets were downloaded by the getGEO function in the GEOquery package ^[47]. 2) To exclude unannotated probes, available probes with clinical information were kept. Once a gene corresponded to multiple probes, the median of the probes was selected ^[48]. 3) Based on the probes with clinical information, the expression matrices were calculated by the exprs function on initial 16,789 gene features for batch 1 data and 24,899 gene features for batch 2 data, respectively. 4) With the alignment of the two batches, only 16,352 common probes were left in both datasets of batches 1 and 2 ^[49].

2.3.   Differential gene screening TabNet (DGS-TabNet)

The pipeline of DGS-TabNet mainly had two modules: one is differential gene screening (DGS) and the other is TabNet, as shown in Figure 1. The first module of DGS was implemented to select the significant differential genes by the reduction of genetic expression data dimensionality; the second module of TabNet distinguished significant differential genes as susceptible genetic features for binary and multi-class classifications, and output the global important genetic features for the further biological interpretability of AD ^[50].

Figure 1.  The pipeline of the DGS-TabNet.

DownLoad: Full-Size Img PowerPoint

2.3.1.   Differential gene screening (DGS)

For expression genetics, DGS was recognized as one of the most popular methods to discover and understand the molecular pathways behind AD ^[51]. In our scheme, it was used to label differential genes and played a crucial role in AD classification.

It was performed by the package of limma ^[52] and the three main steps were listed as follows: 1) The least-squares fit was conducted by the lmFit function, and the comparison matrices were constructed by the makeContrasts function; 2) The log-fold change rate (logFC) and t-statistic were calculated by the contrasts.fit function; 3) The F-statistic was calculated by the eBayes function. The statistical information of the differential genes was summarized by the Toptable function for the parameters of logFC, AveExpr, P.Value, adj.P.value and the B statistic. Here, LogFC provided the logarithmic value; AveExpr provided the average log expression level of the gene; P.Value reflected the gene's false discovery rate (FDR) of the gene; adj. P.Value represented the FDR correction; the B statistic represented the logarithmic ratio of the differential expression gene. Since multiple T tests might lead to an increase of FDR, the adj.P.Value < 0.01 was chosen to maintain the false/true positive ratio.

2.3.2.   TabNet

The TabNet is known as one kind of neural network for the processing of tabular data. It retains the representational automatic learning characteristics for deep neural networks and preserves the capability of providing predictive analyses. In our study, the TabNet consisted of multiple-step decision units and each step included four parts: attentive transformer, feature transformer, split and ReLU function. At the end of TabNet, the global important genetic features were outputted. The workflow of TabNet was summarized as follows:

1) The susceptible genetic features processed by DGS were inputted to the batch normalization (BN) layer for avoiding overfitting and reducing gradient explosion and gradient disappearance ^[53]. The feature calculation was conducted by the feature transformer layers ^[38] (Figure 2). It consisted of two parts: the first part was trained together on all steps, while the second part was trained separately on each step. The gated linear units (GLU) were used in the feature transformer to determine the information passed to the next layer ^[54]. Here, four GLU blocks (two shared and two independent blocks) were used for robust and parameter-efficient learning. Residual connections were used in the layers to maintain the network stability and to avoid the dramatic changes of variances by multiplying $\sqrt{0.5}$ ^[55]. According to the scheme of the attentive transformer layer, it was calculated as:

Figure 2.  The scheme of feature transformer.

DownLoad: Full-Size Img PowerPoint

where a[i-1] was the portion of the split layer divided from the previous step, ${h}_{i}$(•) represented the FC and BN layers, and p[i] denoted the prior scales and was formulated as:

here, the prior scales were used to indicate the application degree of one feature in the previous step.

2) The mask layer was calculated through the attentive transformer ^[33] (Figure 3). The final output Mask (M) was used to weight the input data by using the fully connected (FC) layer, the BN layer and Sparsemax layer. Where the Sparsemax layer assigned feature weights to each feature and the sum of all feature weights equaled to 1.0 ^[56]. To achieve better sparse results, the sparse regularization was calculated and the value of each step was factored into the overall loss, which was calculated as:

Figure 3.  The scheme of attentive transformer.

DownLoad: Full-Size Img PowerPoint

where $\varepsilon$ was generally assigned a very small value. The sparser of M resulted in smaller L, and meant smaller of loss and vice versa.

3) The vectors outputted from the feature transformer layer were divided into two components by the split layer, one was used for computing the final output of the model by ReLU function, and another was implemented to calculate the Mask matrices for the next step.

2.4.   Experimental setup

In order to validate the performance of our proposed algorithm, it was compared with three novel DLs of MLP, NODE and TabNet and five classical ML models of DT, RF, GBDT, Light GBM and SVM on the data set of GEO for binary classification (AD vs NC) and multi-class classification (AD vs MCI vs NC), respectively. During the performance evaluation, five indices of ACC, AUC, recall (REC), precision (PRE) and F1-Score were measured for quantitative evaluation ^[57]. And the interpretability of DGS-TabNet was also evaluated by enrichment analyses of KEGG and GO.

For classical MLs, it was noted that the differential genes were selected by the DGS and boruta algorithm ^[58]. Here, the experiments for classification were carried out on a PC with Windows 10, 64 bit, with 3.00 GHz Intel^® core^TM i9-10980XE processor. The training and testing groups were randomly divided at a ratio of 1:1, and the repeated iterative training was carried out by 10-fold cross-validation. Here, the normalization was conducted by the minmix function, and automatic tuning of the model was conducted by the tune_model function with the iteration of 100.

Similar to MLs, the DGS was also used to screen differential genes for DLs comparison. During the stage of pre-training, the experimental hardware were 2 chips of NVIDIA RTX6000 GPU and the models of DLs were built on the pytorch framework. The training and test groups were randomly assigned at a 1:1 ratio, and the learning rate was set to 1e-3 with the maximum epoch of 200. Besides, the strategies of dropout and early stopping (patience = 5) were used to alleviate the overfitting of our models. The main hyperparameters were listed as: the width of the decision prediction layer (n_d): 32, 44, 52; the width of the attention embedding for each mask (n_a): 10, 14, 36, 38; coefficient for feature reusage in the masks (gamma): 1.0, 1.2, 1.3, 1.5; number of shared GLU at each step (n_shared): 1, 2, 5; sparsity loss coefficient (lambda_sparse): 0.0011-0.0081; optimizer: Adam; number of steps in the architecture (n_steps): 3, 4, 5; number of instances per batch (batch_size): 10. Here, the Bayesian optimization algorithm was used to tune the hyperparameters ^[59,60].

The biological interpretability of DGS-TabNet was conducted by the analyses of KEGG and GO via the package of clusterProfiler ^[61,62]. Here, the enrichment analysis of GO included three aspects of biological process (BP), cellular component (CC) and molecular function (MF). The globally important genetic features outputted by the DGS-TabNet encoder at the best classification were chosen for biological pathways on the two batches of datasets. Here, only the globally important gene features above the weight threshold of 0 were kept for biological interpretability.

3.   Results

3.1.   Results of binary and multi-class classification

The comparison of MLs and DLs for binary classification was presented in Table 2. It was clear that the DGS-TabNet exhibited the best performance with the ACC of 93.80%, AUC of 98.53%, REC of 93.80%, PRE of 93.96% and F1 of 93.79% for batch 2. Besides, the SVM proved to be the best model in the five traditional MLs. The convergence curves of training and testing were shown in Figure 4.

Table 2.  Results of binary classification.

Method	ACC		AUC		REC		PRE		F1
	Batch 1	Batch 2	Batch 1	Batch 2	Batch 1	Batch 2	Batch 1	Batch 2	Batch 1	Batch 2
DT	73.60	87.60	76.29	86.14	92.31	76.19	62.34	97.96	74.42	85.71
RF	76.00	89.15	85.88	96.83	65.38	90.48	73.91	87.69	69.39	89.06
GBDT	75.20	90.70	86.75	97.98	76.92	90.48	67.80	90.48	72.07	90.48
LightGBM	78.40	90.70	85.91	97.47	76.92	87.30	72.73	93.22	74.77	90.16
SVM	88.00	93.80	94.02	98.46	86.54	93.65	84.91	93.65	85.71	93.65
MLP	76.00	70.54	82.90	74.07	76.00	70.54	76.63	71.15	76.14	70.22
NODE	75.20	76.74	87.88	84.60	75.20	76.74	77.84	76.75	73.37	76.73
TabNet	80.00	88.00	89.14	93.89	80.00	88.00	79.95	88.16	79.97	88.36
Proposed	89.60	93.80	95.10	98.53	89.60	93.80	89.72	93.96	89.51	93.79

---

 | Show Table

DownLoad: CSV

Figure 4.  Convergence curves of training and testing for binary classification. (A) batch 1; (B) batch 2.

DownLoad: Full-Size Img PowerPoint

The comparison of transitional MLs and DLs for multi-class classification was displayed in Table 3. Similarly, the proposed approach of DGS-TabNet outperformed the best performance with the ACC of 88.27%, AUC of 94.97%, REC of 88.27%, PRE of 88.43% and F1 of 88.24% on batch 2. It was obvious that the performance of multi-class classification was lower than that of binary classification on both data sets of batches 1 and 2. The convergence curves of training and testing were shown in Figure 5.

Table 3.  Results of multi-class classification.

Method	ACC		AUC		REC		PRE		F1
	Batch 1	Batch 2	Batch 1	Batch 2	Batch 1	Batch 2	Batch 1	Batch 2	Batch 1	Batch 2
DT	53.33	65.92	63.85	78.44	50.23	65.04	52.09	68.55	52.43	66.47
RF	63.64	70.39	80.45	88.98	60.98	70.63	63.92	71.12	62.87	70.62
GBDT	61.21	85.47	79.42	94.53	58.05	84.97	61.54	85.72	60.59	85.53
LightGBM	65.45	81.01	82.32	94.17	63.10	80.40	64.71	81.29	64.94	81.10
SVM	72.73	78.77	87.21	92.21	72.53	78.30	72.91	79.59	72.79	79.00
MLP	55.76	46.37	74.89	60.17	55.76	46.37	58.31	33.89	55.73	33.85
NODE	56.97	54.19	78.27	70.65	56.97	54.19	61.88	52.92	54.70	53.09
TabNet	61.82	49.72	74.41	62.85	61.82	49.72	61.06	52.90	60.22	47.86
Proposed	80.00	88.27	88.77	94.97	80.00	88.27	79.90	88.43	79.68	88.24

---

 | Show Table

DownLoad: CSV

Figure 5.  Convergence curves of training and testing for multi-class classification. (A) batch 1; (B) batch 2.

DownLoad: Full-Size Img PowerPoint

3.2.   Interpretability of the DGS-TabNet

In Table 4, the binary classification of batch 2 showed the best performance. There were the 22 global important genetic features left, and the two largest global important feature weights of AVIL and NDUFS4 were 0.3701 and 0.3077, respectively.

Table 4.  Global important genetic features.

Genetic features	Feature weights	Genetic features	Feature weights
AVIL	0.370100	LOC652864	0.005331
NDUFS4	0.307700	LOC649864	0.005196
REXO2	0.065700	CD79B	0.004798
RPS25	0.053190	WDFY2	0.004436
KCNG1	0.051161	CMTM2	0.003629
SNRPB2	0.029107	RPL3	0.003358
EGFL6	0.027134	RSL1D1	0.003124
SNORD33	0.026315	EIF3E	0.000668
MAP3K6	0.016846	TBC1D2B	0.000256
LOC653702	0.012487	KIAA1160	0.000132
HNRNPH2	0.009302	ASNSD1	0.000002

---

 | Show Table

DownLoad: CSV

The biological pathways enriched by KEGG were found to be related to 19 varied diseases including Ribosome, Coronavirus disease (COVID-19), Huntington disease, Parkinson disease, AD etc. (Figure 6). Especially, the rich factor of AD was about 0.0025 and showed weak correlation with the AD. While the GO analysis enriched a total of 48 BP, 24 CC and 17 MF for biological processes of AD (Figure 7).

Figure 6.  The KEGG enrichment analysis of global important genetic features.

DownLoad: Full-Size Img PowerPoint

Figure 7.  The GO enrichment analysis of global important genetic features. The x-axis and y-axis stand for RichFactor and GO items, respectively.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

Till now, a variety of classical MLs and novel DLs have been used to distinguish AD on gene expression data. The MLs generally relied on feature engineering to downscale high-dimensional gene expression features to achieve a more accurate classification of AD. While it could not be generalized by the limited performance in clinical practices. To this purpose, our proposed DL algorithm was implemented to improve the performance of AD diagnosis on tabular gene data, and even was used to exploit the intrinsic association of gene expression features with AD pathology. It was verified that our proposed DGS-TabNet presented the best performance for AD classification on public gene expression data.

Since the gene expression data usually contained redundant and irrelevant features, the purpose of DGS was generally used to select the most relevant subset of features from whole features. It was clear that the feature selection decided the performance of classifiers and the DGS applied to genetic data had aroused the concerns of many studies ^[63]. During the performance evaluation, the traditional SVM performed better performance than four other MLs. This was because the SVM had the advantage of being robust against noise in genetic data and deducing the solution unique for linearly separable problems ^[64]. For DLs, the DGS-TabNet presented the best performance for AD classification. It was due to the fact that the DGS module could effectively select the deferential genes for biological interpretability between normal phenotypes and disease phenotypes. Without reliable feature extraction, the MLP and NODE preserved low accuracies for AD classification.

Learned from the KEGG and GO enrichment analyses, the globally important genetic features were used to explore the biological interpretability of the proposed DL model. It was proved that one of the global important genes of AVIL actively participated in the protein encoding of the gelsolin/villin family of proteins and it might play a great role in the development of neuronal cells ^[65]. Whereas the neuropathological features of AD included neurogenic fiber tangles with hyperphosphorylated tau proteins ^[66]. These processes included the proliferation, differentiation and maturation of neural stem cells and the regulation of their synaptic and neurotransmission-related processes through interactions. It was suggested that the AVIL gene might be potentially relevant to the formation of AD. According to another NDUFS4 gene, it was thought that the mutations of NDUFS4 would lead to mitochondrial complex I defects ^[67]. During the KEGG enrichment analysis, the gene NDUFS4 was reported to correlate with 12 biological neurodegenerative pathways, including Huntington's disease, amyotrophic lateral sclerosis, and AD, etc. The NDUFS4 enriched in 1 BP, 3 CC3 and 3 MF was proved to be correlated with mitochondrial NADH dehydrogenase and the mitochondrial integrity declined with age and affected a variety of brain functions, such as memory, learning and sensory processes ^[68]. Other 20 global important genes were also associated with molecular pathways involved in the development of AD such as neuroinflammation, oxidative stress, defects in mitochondrial dynamics and function, cholesterol and fatty acid metabolism, and impairment of glucose energy pathways in the brain ^[69]. These enrichment analyses suggested that the globally important genes were potentially highly relevant to the pathogenesis of AD, and it might serve as strong support for the biological interpretability of the proposed model.

Generally, there were several limitations in our study. First, there was no gold standard genetic data available for reliable method evaluation. Due to the data inconsistency, more large sample data sets should be used to verify the generalizability and robustness of our proposed DL algorithm. Secondly, for the betterment of robustness, more types of DL models should be investigated to assess the biological interpretability of AD. In the future, the AD prediction on longitudinal genetic data would be an important prospect for AD diagnosis.

5.   Conclusions

It was concluded that the proposed DGS-TabNet could well deduce the susceptible genes, and the biological interpretability of susceptible genes also revealed the potential possibility of being AD biomarkers.

Acknowledgments

This study was funded by the National Natural Science Foundation of China (Nos. 61971275, 81830052, and 82072228), the grants of the National Key Research and Development Program of China (2020YFC2008700) and Shanghai Municipal Commission of Science and Technology for Capacity Building for Local Universities (23010502700).

Conflict of interest

The authors declare there is no conflict of interest.