Gene expression datasets from the prefrontal cortex brain and blood were obtained from GSE33000 (Narayanan et al., 2014) and Alzheimer's Disease Neuroimaging Initiative (ADNI) (adni.loni.usc.edu) (Petersen et al., 2010), respectively. Brain and blood gene expression data were generated using the Rosetta/Merck Human 44k 1.1 microarray (Narayanan et al., 2014) and the Affymetrix Human Genome U 219 array (Petersen et al., 2010), respectively. We used 113 AD and 244 control (CN) samples for the blood dataset. To maintain consistency between the blood and brain datasets, mild cognitive impairment (MCI) samples were excluded from the analysis, as the brain dataset does not contain MCI samples. For the brain dataset, we used 310 AD and 157 CN samples after excluding Huntington's disease samples. The blood and brain datasets consist of 49,386 and 39,328 probes, respectively. The mean imputation was performed on the remaining missing data values in each dataset. We removed the 30% of probes with a low interquartile range (IQR) across samples. To compare gene expression across different platforms, we mapped probe IDs to Entrez IDs using the biomaRt package (Durinck et al., 2009). The probe with the maximum IQR value was selected for multiple probes annotated with the same gene. By using a MinMaxScaler, values in each dataset are preprocessed to fit within the range of 0 and 1. As a result, 13,666 and 12,319 genes remained in the brain and blood datasets, respectively. A total of 8,922 genes were common between the brain and blood datasets, and they were used in this study.

Pathway information was compiled from the GO biological process (BP) (Ashburner et al., 2000) and KEGG (Kanehisa and Goto, 2000) databases obtained from the Molecular Signature Database (MSigDB) (Subramanian et al., 2005). A total of 168 KEGG pathways and 4,026 GO BP were selected. They include ten or more genes among the 8,922 genes used in this study. Using each of the two databases, a pathway information matrix M∈ℝ^m×n is constructed, which is a sparse matrix representing the relationship between pathways and genes (Equation 1), where m is the number of pathways and n is the number of genes.

To aid in identifying transcriptional features regulating AD, we developed a deep neural network-based classification model named PINNet, which predicts AD based on pathway information. Figure 1 shows the structure of the PINNet. Contrary to previous deep learning approaches, PINNet does not rely on operating as a black box but instead focuses on prior biological knowledge from the GO or KEGG databases. PINNet consists of four layers: an input layer, a pathway layer, a fully connected hidden layer, and a softmax output layer, and the model is built to predict AD (see Figure 1 and Methods). Blood or brain transcription datasets are fed into the input layer of the model, and the output layer indicates a sample status (AD or CN). We formulated the structure of the model from genes to biological functional groups, enabling biological interpretation. The use of multiple levels of ontology (ranging from 10 to 1,156 genes per pathway) acknowledges the diverse complexity of biological processes. During training, the expression data of each gene are induced to patterns of pathway-level of activities, enabling in silico investigation of the biological processes underlying transcriptome-disease association. Genes are partially propagated through the same biological process that contains them, giving rise to functional changes at the pathway level, consequently predicting sample status. Using this design, the PINNet embedded in GO includes 4,026 biological process terms, and the corresponding model for KEGG includes 168 pathways. To incorporate the effect of the pathways, we partially masked the weights of the pathway layer using a binary matrix that represents relationships between pathways and genes. As these neurons' connection is designed to represent biological processes, pathway nodes only receive inputs from the genes included in the pathway. On average, one GO BP is connected to 84 input gene nodes and 42 input gene nodes for KEGG. We normalized pathway nodes with the size of pathways to ensure that the number of genes included in the pathway does not affect the importance of the pathway node. In addition, since the pathway node is very sparsely connected, there is a difference in the range of values of the pathway node and the fully connected node. Therefore, layer normalization was applied to the mini-batch when pathway nodes and fully connected nodes were concatenated. Afterward, they pass through one hidden layer and softmax to predict AD.

PINNet is a four-layer deep neural network with biological pathway information as prior knowledge for AD classification. The PINNet comprises four layers: an input layer, a pathway layer, a hidden layer, and an output layer, as shown in Figure 1. Gene expression values of a sample g∈ℝⁿ are allocated to input nodes. Thus, the number of input nodes (n) is the same as the number of genes (n = 8, 922) in the dataset. The pathway layer consists of fully connected nodes f and pathway nodes p. Pathway nodes p represent relations between genes and pathways, and fully connected nodes f represent features for all genes regardless that they are included in the pathway. The number of nodes in p∈ℝ^m was the same as the number of pathways used. For instance, m is 7,470 when using the GO database and 186 when using the KEGG database. In order to represent the relationship between the gene and the pathway, we obtain pathway nodes p as

where W_p is weight matrices. The Wp∈ℝm×n is masked by the pathway information matrix M∈ℝ^m×n with element-wise multiplication. Masking weights (W_p°M) produce a sparse network by weighting between a particular pathway and genes. Genes not included in the pathway are set to zero. Note that non-zero weights of W_p are updated by backpropagation. u∈ℝ^m is a normalization vector to prevent the values of pathway nodes are affected by the number of genes in the pathway and is defined as follows:

For the fully connected nodes f, n input nodes are directly connected to the f through a linear layer as follows:

where W_f is a weight matrix. p and f are concatenated after layer normalization ([LayerNorm(p); LayerNorm(f)]), and it passes through the hidden layer and the output layer. Then, output o∈ℝ is obtained. Tanh is used for the nonlinear transforming function, and batch normalization is utilized to decrease the effect of internal covariate shift induced by various weight scales. Dropout (α = 0.3) prevents overfitting of the model. We performed the training process by minimizing a cross entropy loss.

The training dataset was oversampled using SMOTE (Chawla et al., 2002) to avoid data imbalance. SMOTE oversampling creates synthetic instances and generates a new instance in line, a segment of the randomly selected k-nearest neighborhood of the minority class. To train PINNet, we initialized all weights with the Xavier initialization. We optimized the loss using the ADAM optimizer, the ReduceLROnPlateau scheduler, and a stochastic gradient descent algorithm, with a mini-batch size of 64. Note that the same settings are applied to both PINNet and the neural network model used for the comparison. The model was implemented using Python (version 3.7.4, https://www.python.org/) and Pytorch (version 1.8.1, https://www.pytorch.org/) on GeForce TITAN X GPUs.

We interpreted the model to identify highly contributing genes and pathways using the SHAP (SHapley Additive exPlanations) values. The SHAP values for each model were computed using the DeepExplainer SHAP package in Python, which is based on Deep SHAP (Lundberg and Lee, 2017). DeepExplainer calculates the feature importance of a given input to predict the deep neural network. The SHAP value is close to 0 if the feature does not affect the model's expectation. Deep SHAP explains the difference between output y and reference output ȳ by the summation of the difference between inputs x and reference inputs x̄. Reference output ȳ are determined by executing forward passes in the network under reference inputs x̄. Let △x_i be the difference between input x_i and reference input xī and C△xi(l) be the attribution of node i in the lth layer. Then, the difference between output y and ȳ can be explained by the summation of attribution scores of x_i in the input layer l:

Reference inputs and outputs can be chosen for a given problem, and we set the reference to a training set with a balanced label distribution, where the reference output is the average of the outputs of the training set. The attribution score calculated from Deep SHAP can be adapted to an approximation of SHAP values. To obtain the contribution of the genes in the input layer, x_i in Equation (5) is gene expression values in a test set. Furthermore, to compute the contribution of pathway nodes in the pathway layer, the values of pathway nodes p in Equation (3) are used as x_i in Equation (5). For each pathway and brain and blood datasets, ten models are constructed by 10-fold cross-validation. Absolute mean SHAP values are normalized as z-scores, and the z-scores are averaged for the ten models. We defined it as an importance score. Importance scores are used to evaluate the contribution of the genes and pathways. The steps for evaluating and analyzing highly contributed genes are depicted in Figure 1.

We applied PINNet to brain gene expression data (Narayanan et al., 2014) and blood gene expression data from ADNI (Petersen et al., 2010) and constructed four models using two pathway datasets: blood with GO [blood (GO)], blood with KEGG [blood (KEGG)], brain with GO [brain (GO)], and brain with KEGG [brain (KEGG)]. While constructing a sparse pathway network, we normalized the pathway nodes with the size of pathways. To identify an appropriate normalization method, we tried three types of a normalization vector u: (1) u = 1 (without normalization), (2) u = 1/n, and (3) u=1/n, where n is the number of genes in each pathway, and compared the importance score of pathway nodes (Details about an importance score are described in Methods). Without normalization, pathways with a larger number of genes had higher importance scores in all four models. The absolute values of the Pearson correlation coefficient between the numbers of genes in pathways and importance scores were strongly correlated in the blood (GO) (|R| = 0.43) and blood (KEGG) (|R| = 0.30). However, when normalizing with n, |R| < 0.1 was obtained in four models. This implies that the number of genes in the pathway has almost no effect on the importance scores of the pathways through normalization (Supplementary Figure S1). Based on this result, the value of the pathway node was normalized with n. We trained PINNet to predict the disease status of samples (see Methods) and tested its performance by 10-fold-cross-validation. In the cross-validation, each fold preserved the ratio of sample classes, and the ratios of AD to CN were 0.46 and 1.97 in blood and brain, respectively. A part of the training set was used as a validation set so that the ratio of training, validation, and test set was 8:1:1. During the training of PINNet, a validation set was used for early stopping and selection of the number of fully connected nodes and learning rates with grid search (Supplementary Table S1). We selected the hyperparameters that yielded the best AUC on the validation set in each fold. Early stopping was applied if the validation AUC did not increase by more than ten epochs during the training. The number of fully connected nodes in the pathway layer is a hyperparameter with a set of candidate values of 32, 64, and 128, smaller than the number of pathways. The number of nodes in the second hidden layer was fixed as 64. Candidate values for the learning rate are 0.0001, 0.0005, and 0.001. All compared methods were trained and tested on the same data splits.

With the brain gene expressions from the prefrontal cortex area (Narayanan et al., 2014), the area under the curve (AUC) values between true positives and false-positive rates were 0.9744 and 0.9763 with GO and KEGG pathways, respectively (Table 1). The previous study using the same dataset showed the accuracy of around 86.30–91.22% (Cheng et al., 2021) using several machine learning models, reporting that brain gene expressions are valuable to classifying AD patients. This was also confirmed in our study. With the blood gene expressions from ADNI, we obtained AUC values of 0.6355 and 0.6420 with GO and KEGG pathways, respectively (Table 1). Predictions using the blood gene expressions showed relatively lower classification performance than brain expression data. This range of performance was also reported in the previous study (Lee and Lee, 2020).

We compared PINNet with deep neural network (DNN), support vector machines (SVM), and random forest (RF) methods. The DNN model architecture is identical to PINNet with the same number of the first hidden layer, except that it does not apply a pathway binary mask. The DNN model consists of two fully connected layers without a pathway binary mask, and the number of nodes and learning rate of the first hidden layer were selected using grid search in the validation set. Since the range of nodes in the first hidden layer of PINNet was 200 (168 + 32) to 4,154 (4,026 + 128), the candidate nodes in the first hidden layer include 128, 512, 1,024, and 4,096, and the learning rate used the same set of candidate values as PINNet. Similarly, for SVM and RF, the combination of hyperparameters and kernel type with the highest AUC for the validation set was chosen within each fold. Sets of hyperparameter values searched in the validation and subsequently selected best hyperparameter values are provided in Supplementary Table S2. As shown in Table 1, our PINNet outperformed the SVM and RF models, achieving better classification performance on blood datasets. In the brain dataset, the SVM model showed the highest performances, although all models showed high prediction performances. Despite the extra parameters in DNN, which is a fully connected model, the performance of the DNN was lower than PINNet, especially for blood. Additionally, Table 2 shows F1 scores of PINNet and comparing models.

A major objective of the analysis of gene expressions for the disease is to associate changes in the transcript with changes in phenotype. To reveal such associations, we interpreted PINNet models for AD prediction. We identified which genes played an important role in prediction in the models learned from 10-fold-cross validation in Table 1. To determine the contribution of each gene in AD prediction, the feature importance of nodes in the input layer was measured by the SHAP value using the Deep SHAP method, a backpropagation-based deep neural network interpretation algorithm (see the Methods section). For each pathway (GO and KEGG) and each gene expression dataset (brain and blood), we calculated SHAP values from the trained PINNet models for each gene. We calculated the importance score as the average of the normalized SHAP values of the ten models. As a result, we obtained the importance scores of genes for each model of blood (GO), brain (GO), blood (KEGG), and brain (KEGG).

We analyzed the importance scores of genes using two lists of known AD-related genes: AlzGene (Bertram et al., 2007) and DigSee (Kim et al., 2017). The AlzGene contains 681 genes curated by systematic meta-analyses, and 361 genes were shared with the 8,922 genes used in the study. DigSee contains 1,635 AD-related genes extracted from PubMed abstracts using text-mining techniques, and 961 were common with genes used in this study. We defined genes included in AlzGene or DigSee as AD-related genes and the remaining as non-AD-related genes. We used the Wilcoxon test to compare the importance score between AD-related and non-AD-related genes in PINNet and DNN. As shown in Figure 2, AD-related genes screened from DigSee and AlzGene showed significantly higher contributions than non-AD-related genes in PINNet models for all datasets (p-value ≤ 0.001). However, in the case of DNN models, the difference in importance score between AD-related genes and non-AD-related genes is only partially significant for brain data, and there is no difference in blood. In addition, we observed that the importance scores of known AD genes from PINNet were significantly higher than those from the DNN model for both KEGG and GO sources. These results shows the robustness of our approach. Unlike standard DNNs, PINNet captures important biological features based on prior biological knowledge. It represents that PINNet's biological prior knowledge leads to learning disease-related genes as essential features.

Using a protein interaction network obtained from the STRING database, we examined genes with high importance scores. We selected 1% genes with the highest importance scores from the four models. We integrated genes selected in PINNet (GO) and PINNet (KEGG) models for the blood and brain, separately, and obtained 170 and 167 genes in blood and the brain, respectively. The highly contributed genes were mapped onto the STRING database (version 11.0) (Szklarczyk et al., 2019), where interactions larger than a confidence score of 0.7 were used. Cytoscape (Shannon et al., 2003) was used to configure the PPI network. As shown in Figure 3, 152 and 138 genes in the blood and brain were mapped in the largest network, respectively (Supplementary Table S3). We also used ClusterOne (Nepusz et al., 2012) to screen PPI subnetwork modules. ClusterOne identifies overlapping functional modules by detecting dense regions in the protein interactome network. The parameters were set to p-values < 0.05, the minimum size = 5, and an edge weight as a confidence score in the STRING database. If genes of two subnetworks were overlapped by more than 70%, we merged them into a large module. Then, modules were selected if more than one biological term was enriched with GO BP and KEGG using DAVID (Sherman et al., 2009) (Bonferroni corrected p-values < 0.05). Figure 3 shows blue, green, red, orange, and purple modules in the blood and blue and green modules in the brain. In the blood, the blue module was significantly enriched in 218 terms related to proliferation and several signaling pathways, including Ras, FOXO, and PI3K-Akt signaling pathways. The green module was enriched with antigen processing and presentation terms, especially the MHC class II protein. In the red module, the serine family, comprising cysteine, glycine, and serine metabolism-related terms, was enriched. The orange module was related to several metabolic processes, including glycolysis, ethanol catabolic process, and aldehyde catabolic process. The purple module was enriched with several carbon metabolism pathways. In the brain, the blue module was similar to the blue module in blood. It was enriched with 350 terms related to proliferation and several signaling pathways, including TNF, FOXO, and PI3K-Akt signaling pathways. The green module contained 29 enriched terms related to the oxaloacetate metabolic process (Details about enrichment tests are in Supplementary Table S3).

We also identified highly contributed genes that show differences between blood and brain in AD prediction. First, we compared models based on the same pathway database (GO, KEGG) in different tissues with a two-sided Wilcoxon test (see Supplementary Table S4; Supplementary Figures S2, S3). From the results of 10-fold cross-validation, absolute SHAP values with z-score normalization were used for the test, and the p-values were corrected by Benjamini–Hochberg method (BH). We obtained 24 and 291 significantly different genes between blood and brain models (blood (GO) vs. brain (GO) and blood (KEGG) vs. brain (KEGG), respectively) with the adjusted p-values (BH) < 0.05. Next, we checked highly contributed genes (top 1%) that showed a clear distinction between the tissues. As a result, in the brain, CLU, PTPN11, ITPKB, PRKACB, IGF1, and MTHFD2 were more important compared to blood. CLU is related to several pathological states of AD, including a mediator of Aβ toxicity (Foster et al., 2019). PTPN11 (also known as SHP2) interacts with tau in Alzheimer's disease brain (Kim et al., 2019). Overexpression of ITPKB induces the increase of Aβ40 and tau hyperphosphorylation in a mouse model of familial Alzheimer's disease (Stygelbout et al., 2014). PRKACB is involved in the elevation of Aβ and tau hyperphosphorylation levels (Wang et al., 2019). Reduced IGF-1 signaling in the AD mouse model is attributed to a decline in neuronal loss and behavioral impairment (Cohen et al., 2009). MTHFD2 is differentially expressed in AD posterior cingulate astrocytes (Sekar et al., 2015). In the case of blood, although two genes were significantly important compared to brain, there is no known evidence to explain the relevance of these genes and AD in blood.

We further investigated whether genes selected using the importance scores in PINNet can be relevant features for predicting AD. First, we divided all samples into ten sets; eight sets were used for training, a set for selecting genes (named a feature selection set), and a remaining set for testing (see Supplementary Figure S4). Second, after training PINNet using the training set, we selected 892 (10%) genes with the highest importance scores in the feature selection set. Third, a multilayer perceptron (MLP) model was trained with the 892 selected feature genes, and then the remaining test set was used to evaluate the AD prediction performance. We repeated this process ten times. The MLP model consists of four layers, and the numbers of nodes in each layer are {n, n/2, n/4, 2}, where n is the number of input genes. We compared the performance of 10% of selected genes using PINNet with the same number of randomly selected genes and all genes for predicting AD samples. Genes selected using PINNet outperformed the random gene set in blood models (Table 3), suggesting that only 10% of the genes selected from PINNet can represent the whole blood transcriptome dataset.