A Pathway Based Classification Method for Analyzing Gene Expression for Alzheimer’s Disease Diagnosis

Abstract

Background:

Recent studies indicate that gene expression levels in blood may be able to differentiate subjects with Alzheimer’s disease (AD) from normal elderly controls and mild cognitively impaired (MCI) subjects. However, there is limited replicability at the single marker level. A pathway-based interpretation of gene expression may prove more robust.

Objectives:

This study aimed to investigate whether a case/control classification model built on pathway level data was more robust than a gene level model and may consequently perform better in test data. The study used two batches of gene expression data from the AddNeuroMed (ANM) and Dementia Case Registry (DCR) cohorts.

Methods:

Our study used Illumina Human HT-12 Expression BeadChips to collect gene expression from blood samples. Random forest modeling with recursive feature elimination was used to predict case/control status. Age and APOE ɛ4 status were used as covariates for all analysis.

Results:

Gene and pathway level models performed similarly to each other and to a model based on demographic information only.

Conclusions:

Any potential increase in concordance from the novel pathway level approach used here has not lead to a greater predictive ability in these datasets. However, we have only tested one method for creating pathway level scores. Further, we have been able to benchmark pathways against genes in datasets that had been extensively harmonized. Further work should focus on the use of alternative methods for creating pathway level scores, in particular those that incorporate pathway topology, and the use of an endophenotype based approach.

INTRODUCTION

The most common form of dementia is Alzheimer’s disease (AD). It is predicted that by 2050, 1 in every 85 people will be living with the disease [1]. No disease modifying treatments are available for AD and existing treatments only provide short-term symptomatic relief in a subset of patients [2]. Additionally, in the early stages (between 2 and 15 years prior to the development of clinical symptoms) the disease is difficult to diagnose. Villemagne et al. and Jack et al. hypothesize that characteristic AD pathology (the presence of amyloid-β (Aβ) plaques and hyperphosphorylated tau tangles in the brain) begins to develop up to 20 years prior to clinical diagnosis [3, 4]. This extended prodromal stage is an important window in which to target treatments that may be able to alter the course of the disease; provided people could be sensitively and accurately diagnosed. Aβ, tau, and phosphorylated-tau levels are indicative of AD pathology in this prodromal period and can be measured in cerebrospinal fluid (CSF) and by positron emission tomography (PET) imaging [5]. The procedures involved in attaining these measurements can be invasive or expensive and require specialized administration, equipment. and expertise. The development of a less invasive, potentially cheaper technique, such as a blood test, would offer significant advantages [6].

Recent studies indicate that gene expression levels in blood may be able to differentiate AD subjects from normal elderly controls and mild cognitive impairment (MCI) subjects with prodromal disease [7–10]. Han et al. provide an overview of studies of gene expression associated with AD-related phenotypes [11]. They state that the blood transcriptome is vital in the disease mechanism of AD and should therefore be investigated further in independent studies of a large sample size. A more general summary of gene expression data in neurodegenerative diseases is given by Cooper-Knock et al. [12]. This review emphasizes the dysregulation in neuroinflammation and intracellular signaling pathways including calcium signaling in AD. The commonality between these reviews is that they both highlight limited replicability at the single marker level. Furthermore, Han et al. report a greater concordance between differentially expressed genes at the pathway level. A pathway-based interpretation of gene expression may therefore prove more robust across different sample populations. Such an approach may also reduce noise and dimensionality. It is important to note that differential gene expression, as described in these reviews, does not necessarily identify genes that will be useful in a classification context.

Although previous gene expression studies in AD have retrospectively identified pathways altered in disease [9], this is the first study to use pathway scores for each individual to build predictive models across the population. This study used Pathway Level Analysis of Gene Expression (PLAGE) to estimate pathway variability across samples in the population by calculating sample-wise pathway scores [14]. PLAGE outperformed other single sample enrichment methods such as ZSCORE, Gene Set Variation Analysis (GSVA), and Single Sample Gene Set Enrichment Analysis (SSGSEA) in a comparison of sensitivity, specificity, and prioritization by Tarca et al. [15]. PLAGE scores have been used in univariate t-testing and unsupervised clustering methods to investigate the pathways involved in oral leukoplakia and those leading to cell proliferation and migration in leukemia [16, 17]. We combine, for the first time, PLAGE scoring with a supervised machine learning approach to build an AD classifier.

This study used blood expression data from subjects participating in the AddNeuroMed (ANM) and Dementia Case Registry (DCR) studies to develop models of clinical diagnosis. The performance metrics of gene expression and demographic models is compared with those generated using pathway level measures of expression.

MATERIALS AND METHODS

Cohort

ANM is a European multi-center study aiming to develop biomarkers for AD [18]. Subjects with an AD diagnosis as well as those with MCI and healthy controls were recruited from centers based in Kuopio, Lodz, London, Perugia, Thessaloniki, and Toulouse. Details of study design and enrollment are provided by Lunnon et al. [10]. Subjects for the DCR were recruited from the Maudsley and Kings Healthcare Partners, which incorporates the Alzheimer’s Research UK (ARUK) cohort [19] from whom gene expression data has not previously been reported.

The present study used data from 748 subjects: 614 subjects from ANM and 134 subjects from DCR.

Gene expression

Whole blood samples (2.5 ml) were collected after 2 h of fasting into Paxgene Blood RNA tubes (BD) and extracted as in Lunnon et al. [9]. Illumina Human HT-12 Expression BeadChips were used to analyze the whole transcriptome according to the manufacturers protocol. The gene expression analysis was run in two batches at two different sites. Batch 1 contained samples from 356 ANM subjects run on version 3 of the BeadChip, as previously described [9, 10]. Batch 2 contained samples from 411 subjects: 134 from DCR and 277 from ANM run on version 4 of the BeadChip. Samples from 19 subjects were included in both batches. See Fig. 1 for an overview of sample numbers. The raw gene expression data are available as GEO DataSets (Accession number GSE63060 for batch 1 and GSE63061 for batch 2).

Statistical analysis

Data pre-processing

The data pre-processing performed in this study is different to that used for the original analysis by Lunnon et al. [9, 10]. The data processing pipeline used in this study aims to address the effects of technical data artifacts in gene expression studies [20]. Raw gene expression data was exported from Illumina’s Genome studio and processed in R (version 3.1.1) [21] using the lumi package [22] and custom in-house pre-processing scripts (GitHub, http://bit.ly/1vjyKNo). Briefly, raw expression data was subject to a model based background correction for bead array [23]. This used negative bead expression levels to correct for background noise. The data was then log base 2 transformed and robust spline normalized in lumi [22]. Outlying samples were iteratively identified using fundamental network concepts and removed, following the methods described by Oldham et al. [24]. To reduce any batch effects we adjusted for technical categorical variables using ComBat [25]. Continuous technical artifacts were accounted for by taking the first principal component across housekeeping and undetected probes and regressing this against technical variables. Variables significantly associated with the first principal component were then regressed against expression for each probe, and the mean adjusted residuals taken forward for all further analyses. Finally, the data was reduced to a subset of probes that could be reliably detected in 80% of samples in at least one diagnostic group. Finally, subjects were excluded where there were discrepancies between the recorded sex and sex determined by the XIST (ILMN 1764573), USP9Y (ILMN 2056795) and EIF1AY (ILMN 1755537 and ILMN 2228976) X- and Y-linked genes.

Demographic data for the ANM and DCR subjects was extracted using CohortExplorer [26].

Pathway level analysis of gene expression (PLAGE)

Gene level expression data were condensed to sample wise, pathway level scores using PLAGE [14]. PLAGE groups genes into pathways defined by the Broad Institute Collection of Curated Pathways [27] and outputs a score, per sample, for each of these sets. We restricted PLAGE to only include pathways with between 10 and 500 genes. The generation of PLAGE scores was implemented through R package ‘GSVA’ and is detailed in Supplementary Methods, section 1 [13].

Data analysis

Clinical diagnosis (AD versus non-demented elderly control) classification models were built using batch 1 gene expression data. Variable selection was performed using recursive feature elimination (RFE) and the creation of a tolerance set using the ‘pickSizeTolerance’ function in R. This function finds a smaller set of variables while maintaining model accuracy [28]. Three Random Forest (RF) models were built, the first of which was a model based on demographic data alone (demographic model) [29]. The demographic variables included were those that were significant in the batch 1 population: sample collection site, age, years in full time education, and APOE status (defined as the presence of any number of ɛ4 alleles) (Table 1). Two further models were built based on thesedemographic variables and gene level data (gene model) or PLAGE scores (pathway model). The purpose of the demographic model is to provide a comparator for the gene and pathway models. If models that include blood expression information (as well as demographics) are no more informative than demographic variables alone there is no benefit in including this information. All model building was performed in the statistical software R (Version 3.1.1) using the ‘caret’ package [28].

Each model was used to predict the diagnostic status of subjects in batch 2. Model statistics including accuracy, sensitivity, and specificity were generated and compared between the demographic model, gene model, and pathway model. Receiver Operator Curve (ROC) analysis was also performed in batch 2 data using R packages ROCR and pROC [30, 31].

Full details of model building are provided in Supplementary Methods section 2.

Additionally, variable importance (determined as the change in Gini index) was examined in the pathway model by permutation testing. The idea of permutation testing is to break the association between outcome (in this case diagnosis) and predictor variables. When the model is re-built based on this permuted data any significant association is spurious. Therefore, by comparing the true variable importance of a pathway to the variable importances that arise by chance in the permuted data we can assess how significant our result is. To achieve this here we used 1000 permutations of the demographic variables (including diagnosis) and for each permutation built a RF model. The importance measures of each pathway were then compared to that of the original model to generate an empirical p-value. A p-value of less than 0.05 was considered significant.

The validity of the pathways selected in the pathway model was also investigated in a similar way. A random set of pathways (of the same size as the final pathway model) were selected, and used to build a RF model. The fact that this set is random breaks the association between predictor variables and outcome. This process was repeated 1000 times and the accuracies across all models compared to the accuracy of the true model to create an empirical p-value.

RESULTS

Cohort demographics

Table 1 gives an overview of the demographics of subjects included in the two batches of gene expression data.

Data pre-processing

As a result of pre-processing 12 samples in batch 1 and 49 samples in batch 2 failed quality control (QC) and were removed. The majority of these samples failed QC as they were identified as outliers. Additionally, some samples were removed because the sex of the individual recorded in the clinical database did not match the biological sample (2 samples in batch 1 and 7 in batch 2).

Samples from 19 subjects were present in both batch 1 and batch 2. Samples from 14 of these individuals passed QC in both batches; only data from batch 1 was used and the other was discarded. Correlation between the two batches was at least 0.9 for all individuals (Supplementary Figure 3). Batch 2 gene expression data contains subjects from the DCR whereas batch 1 does not. This study used the same protocols, staff, and facilities as the London sample collection site within ANM. Principal components analysis (PCA) was performed across the batch 2 gene expression data from DCR and ANM subjects from London. The first three principal components (accounting for >40% of variation) were linearly regressed against the study the individual was enrolled in (DCR or London ANM) and found to be non-significant. Therefore, it was deemed appropriate to group DCR subjects with London ANM, allowing the model trained in batch 1 data to be simply applied to batch 2 data.

After data processing only subjects with either an AD diagnosis at all visits or control status at all visits were analyzed further: 207 subjects in batch 1 and 236 in batch 2.

Only gene probes that mapped between the version 3 and version 4 chips used to generate batches 1 and 2, respectively, were used for analysis (5212 probes). The Broad Institute Collection of Curated Pathways matched these probes to 834 pathways [27].

Data analysis

Demographic model

The following demographic variables were associated with case/control status in our cohorts (Table 1): age, sex, APOE status, years in full time education, and sample collection site. These variables were therefore used in multivariate modeling using RFE. The optimal cross-validated accuracy was found when including all variables; calculation of a tolerance set excluded the variable representing the Lodz sample collection site. Variable importance scores showed age as the most important covariate followed by years in full timeeducation and then APOE status and sample collection site. In batch 2 test data the model achieved an accuracy of 0.69, sensitivity of 0.53 and specificity of 0.84.

The area under the ROC curve was 0.77 (see Table 2 and Fig. 2).

Additionally, a model that did not contain the sample collection site was built. The aim was to create a model based on demographics that would be available to clinicians. This model had a slightly decreased accuracy in comparison to the demographic model but outperformed the pathway model and gene model in accuracy, sensitivity and area under the ROC curve at 0.67, 0.68 and 0.76, respectively. Interestingly, the specificity of the model was lower than all others at 0.67 (Table 2).

Gene model

The top 5% of variables from the bootstrapped variable importance calculations (261 variables) were carried forward to the RFE model building. The optimal cross-validated accuracy from RFE in the gene model was found for all of the 261 variables; calculation of a tolerance set reduced this set to only 13, excluding all demographic variables. For a list of genes, see Table 3.

In batch 2 test data, the gene model accuracy was lower than that of the demographic model and equal to the pathway model. The sensitivity, specificity, and area under the ROC curve of the gene model lay between the demographic and pathway models at 0.59, 0.75, and 0.72, respectively. (see Table 2 and Fig. 2). Note that the pathway model showed higher sensitivity while specificity and AUC ROC were higher in the demographic model.

Pathway model

The top 5% of variables from the bootstrapped variable importance calculations (42 variables) were carried forward to the RFE model building. The optimal cross-validated accuracy from RFE in the pathway model was found for 40 of the variables; calculation of a tolerance set reduced this set to only 16 variables (Table 4): 14 pathways, age, and years in full time education.

Permutation tests of variable importance were performed to assess the size of effect relative to that observed under the null hypothesis of no association. Of the 14 pathways, two achieved nominal significance with a p-value <0.05 and are indicated with a  ^* in Table 4 (Supplementary Figure 5).

Additionally, we compared the model accuracy of 1000 models comprising 16 random variables from the pathways, age, sex, APOE status, and years in full time education. This yielded a p-value of 0.082 indicating that, statistically, the final model does not perform significantly better than a model of random pathways (Supplementary Figure 4).

In batch 2 test data, the model accuracy was lower than that of the demographic model at 0.66, however, the sensitivity was higher at 0.61. Both specificity and area under the ROC curve were lower than the demographic model at 0.70 and 0.73, respectively (see Table 2 and Fig. 2).

There is minimal overlap in genes between the different pathways included in the final pathway model. This is illustrated by the sparse percentage overlap map shown in Fig. 3 and supports the idea that each pathway is contributing an independent signal to the model. Of the 13 genes included in the gene model, only four of them (ILMN 1776104, ILMN 1784286, ILMN 1695645, and ILMN 2048326) appear in any of the pathways in the pathway model.

Misclassification

We discovered that 22% of controls used in the training data had reported memory complaints deemed not serious enough to reflect a change in diagnosis. By studying misclassification rates split by AD subjects, control subjects, and control subjects with memory complaints, we see that the most well classified group in the gene model was those subjects with memory complaints whereas in the pathway model it was control subjects (see Supplementary Figure 6). We also demonstrated that time since disease onset is not related to misclassification of AD subjects and control subjects with memory complaints in the test data (Supplementary Figure 7).

DISCUSSION

In this study we investigated whether AD cases could be differentiated from control subjects using gene expression data analyzed at the pathway level. We were particularly interested in confirming whether pathway level information created a more robust predictor of case/control status than expression data at the gene level as recent reviews of AD studies have suggested [11]. Our results, using subjects from the ANM and DCR cohorts, show similar model performance in a pathway model compared to a gene and demographic only model. In this study, we do not find improved prediction of AD diagnosis using pathway level information using the PLAGE method to calculate pathway scores. However, the robustness of pathway based approaches for AD biomarker discovery should be tested in other gene expression data from different populations andplatforms.

The fourteen pathways included in the final pathway model focused around DNA repair, immune response, and regulation of cellular activities. Of particular interest to AD, two pathways containing genes from the calpain gene family were included in the final fourteen pathways. It is thought that amyloid peptides interfere with calpain activity leading to deregulation of the CDK5 gene and in turn hyperphosphorylation of the tau protein. This promotes the death of neurons [27]. It is encouraging that we have seen relevant pathways in our final pathway model. Overall, the pathways are similar to those identified by Lunnon et al. who studied overall pathway differences using an identical raw dataset that was processed differently [9]. As we would expect, 12 out of 13 of the genes in the final gene model were present in the genes used for modeling by Lunnon et al. The data had been processed slightly differently emphasizing that these signals are robust to alterations in processing and modelingmethods.

RF models are commonly used in biomarker studies [9, 33]. However, it has been shown that they exhibit variable selection bias being more likely to select continuous variables or those with many categories [34]. Additionally, the presence of correlated predictors (as is common in gene expression studies) can add further bias [35]. Strobl et al. aimed to address these issues with an ensemble-learning algorithm based on conditional inference trees; Conditional RF (CRF) models [36, 37]. We attempted to use this methodology in the present study. We hypothesized that the creation of an unbiased predictor may highlight different pathways and genes to those previously discovered, potentially allowing greater predictive ability. However, the process of creating a CRF model was computationally expensive even when using high performance computing resources. Model building considering the 834 pathways and 5,212 genes was consequently infeasible. Work to improve the efficiency of this method would be computationally beneficial and would allow the use of alternative variable importance measures. Measures such as mean decrease in accuracy and conditional mean decrease in accuracy would be an improvement over biased variable importance measures such as the Gini index, which was used in this study.

This study used the Broad institute collection of curated pathways to generate the pathway model and excluded less well-curated gene sets. This method was chosen due to its performance in a comparison study [15] and due to ease of application through the GSVA R package. It may be beneficial, although potentially computationally costly, to create pathway level scores that also reflect pathway topology and thus add further detail to the model. Such methods have been created by Pyatnitskiy et al. [38] building on the work of others [15, 39]. The method detailed by Pyatnitskiy et al. does not depend on predefined gene sets as used in this analysis. However, it is also unable to control the number of genes in a pathway; a potential benefit of using PLAGE. A further limitation of existing pathway approaches is that they often ignore information on the direction of change for each gene within a pathway. This would be an interesting area for further method development.

The creation of a demographic model that excluded sample collection site led to a drop in accuracy. Although RNA extraction and analysis were performed at one site the blood collection may vary by location. We aimed to correct for batch effects occurring in extraction and analysis in the pre-processing. This highlights that although sample collection sites within multi-center studies are following the same protocols major technical differences can still arise and remain after QC steps including batch correction. As much as possible, these differences should be quantified during extraction. Standardization for future biomarker development will aid this. It is possible that the sample collection site effect we see is driven by genetic differences between sites for some genes (expression quantitative trait loci). For a biomarker to have clinical utility it should be robust to such differences. However, in early exploratory work we are more likely to find results of interest if technical data artifacts are not creating a barrier.

The models created in this study all achieved an accuracy of approximately 70% with the pathway model having test sensitivity and specificity results of greater than 60% . The pathway model and gene model did not outperform a model of demographics alone. Any potential increase in concordance from the novel pathway level approach used here has not lead to a greater predictive ability in these datasets. However, we have only tested one method for creating pathway level scores. Further, we have been able to benchmark pathways against genes in datasets that had been extensively harmonized. It is reassuring to see that pathways perform similarly to genes and further work is now needed to see if pathway concordance is more easily detected using other methodological approaches and in data generated by independent groups andplatforms.

Furthermore, we found that the heterogeneity of control subjects may be leading to reduced predictive accuracy and suggest that the use of an endophenotype may be beneficial in future work.

Conclusions

We have used subjects from the ANM and DCR studies to investigate case/control classification using gene and pathway level expression data. We hypothesized that a model built on pathway level data may be more robust than a gene level model and consequently perform better in test data. However, a pathway level model built using scores and a gene level model performed similarly to each other and to a model based on demographic information only. Further work should focus on the use of alternative methods for creating pathway level scores, in particular those that incorporate pathway topology, and the use of an endophenotype based approach.