Random forests on distance matrices for imaging genetics studies

Data sets

Concerning the simulation studies, for the purposes of realism, the simulated genotypes and quantitative phenotypes are built from the public datasets of the International HapMap Project⁷ and the Alzheimers Disease Neuroimaging Initiative (ADNI)⁸ respectively. The genotypic and phenotypic data that is used for the application on Alzheimer’s disease was also retreived from ADNI.

HapMap data

A human haplotype map of over 1.5 million SNPs for 993 subjects from 11 sub-populations across four continents was obtained from the second release of the HapMap phase 3 dataset⁹ of the International HapMap Project.

ADNI data

The ADNI study began at 2004 as a multicenter effort to develop clinical, imaging, genetic and biochemical biomarkers for the early detection and tracking of AD. The initial 6-year phase of the study (ADNI1) included 400 subjects diagnosed with mild cognitive impairment (MCI), 200 subjects with early AD and 200 elderly control subjects. In 2009, ADNI1 was extended with ADNI GO which assessed the existing ADNI1 cohort and added 200 participants identified as having early mild cognitive impairment (EMCI). Finally, in 2011, ADNI2 began, which assesses new participants in addition to those from the ADNI1/ADNI GO cohort. All participants in the ADNI studies were recruited across North America and are followed and reassessed over time to track the pathology of the disease as it progresses.

For the purpose of applying distance-based Random Forests on a real dataset, we acquired genotypes for 464 subjects from the ADNI database ((99 AD, 154 CN, 211 MCI). Genotyping was performed with the Human610-Quad Bead-Chip microarray, which includes 620,901 SNPs and copy number variations (see Saykin et al. (2010) for details). A separate genotyping of the SNPs for the APOEε4 variant has been performed by ADNI, since these were not included in the above microarray. The subjects were unrelated, of European ancestry, and passed screening for evidence of population stratification using the procedure described in Stein et al. (2010).

The selection and preprocessing steps are taken from Silver et al. (2012), where the same data was used for the identification of possible AD causal gene pathways. In detail, the genotypic data consists from autosomal SNPs with genotyping rate >95% (42,680 SNPs), HardyWeinberg equilibrium p-value >5×10⁷ and minor allele frequency >0.1. For each subject we ended up with a set of 434,271 SNPs. Any missing genotypes were imputed using the same procedure as in Vounou et al. (2012).

The imaging data consists of longitudinal brain MRI scans (1.5 T) from the ADNI database. For each of the subjects, serial brain MR images were acquired at baseline and at 6, 12 and 24 months following initial recruitment. Details about the acquisition protocol and the preprocessing steps that were applied by ADNI on the raw scans can be found in Jack et al. (2008).

Further processing was carried out, as described in detail in Silver et al. (2012). The aim is to provide a measure of brain structural change over time. To achieve that, a linear regression model with an intercept term and time as the independent variable is fitted for each voxel. The coefficient for the slope provides a measure of brain tissue change, relative to time, for each voxel. Correction for age and sex, as well as selection of the most discriminative voxels is performed by analysis of variance (ANOVA) (Silver et al., 2012). The final endophenotype data that we use is in the form of longitudinal maps of ventricular/CSF expansion (reconstructed to 1 mm isotropic voxels), which represent brain tissue loss and are normalised across all the subjects. Each subject’s endophenotype vector consists of 148,023 slope coefficients.

The simulation studies below were performed using information from the subset of CN and AD subjects, while for the application study, MCI subjects were included as well.

Imaging genetics simulation engine

In this section we provide details of the imaging genetics simulation engine. We describe, in order, the genotype simulation, the phenotype simulation, the genetic disease modelling, and lastly the disease classification modelling.

Base genotype simulation

The ADNI dataset with 253 subjects is clearly too small to represent a direct population source for the multiple samples required for the multiple simulation runs. Instead we simulate this base population and its genotype within an individual-based, forward-time, population genetics simulation environment. Specifically we employed the SimuPOP Python package (Peng and Kimmel, 2005) and have made use of the SimuGWAS¹⁰ implementation (Peng and Amos, 2010), originally written for the simulation of samples for GWAS studies. This allows one to evolve the genotype of a given starting population over a specified number of generations with exposure to the standard evolutionary forces of mutation, recombination, migration, and selection.

The pattern of linkage disequilibrium (LD) between markers in a forward-time simulated population is sensitive to the make-up of the starting population. In the absence of data on populations in the distant past, we adopt the HapMap phase 3 subjects as a reasonable proxy. Rather than evolving the entire genotype, it is suffices for our simulation study to restrict our attention to a subset of genetic markers located on a single chromosome. From all 993 subjects in the HapMap dataset, we extract 385 sequential, biallelic, SNP markers from the front end of chromosome 19 with a minimum maf of 0.05. We perform a forward-time evolution of this initial population of 993 individuals over 500 non-overlapping generations using the Wright-Fisher model. We assume a constant population growth rate towards a final population of 10000 individuals. We maintain a recombination rate of 1.0×10^–8, a mutation rate of 1.0×10^–8, and a migration rate of 1.0×10^–3 throughout the simulated evolutionary process.

We denote the genotype of individual i at locus s by the minor allele count x_is∈{0, 1, 2}. From the set of 385 markers, we identify two subsets of 16 SNPs with maf~0.2,¹¹ the purpose of which will be made explicit later in Section 6. Let the two subsets be S_c and S_s, where the subscripts indicate the labels potentially causative and potentially spurious respectively. We emphasise the potentiality in their descriptions to highlight the fact that, depending on the specific genetic model simulated (see later), not all the SNPs from the two sets will be associated to phenotypic modifications.

We consider a single chromosome with 385 markers and denote the genotype of individual i at locus s by the minor allele count x_is∈{0, 1, 2}. We assign two subsets of SNPs, S_c and S_s, which we label potentially causative and potentially spurious respectively. We emphasise the potentiality in their descriptions to highlight the fact that, depending on the specific genetic model simulated (see later), not all the SNPs from the two sets will be associated to phenotypic modifications.

Base phenotype simulation

There is no constraint on the number or type of imaging phenotypes for the simulation model. In this paper we consider two such quantitative phenotypes – a multivariate, vectorial, representation and a matrix representation of imaging phenotypes. These representations are motivated by both the neuroimaging data currently available, and also by the pathology of AD. For instance, the neuronal atrophy in AD patients can be represented as vector of time-averaged rates of volume change of the various ROIs in the brain. Also, as demonstrated in several more recent studies (Sun et al., 2009; Huang et al., 2010; Brier et al., 2012), the number of connections between different ROI is shown to be lower in AD subjects as compared with CN subjects. A similar decrease in connectivity is also observed in longitudinal FDG-PET and resting-state functional MRI (rs-fcMRI) scans of individual AD subjects (Sun et al., 2009; Brier et al., 2012). Although there are exceptions at the beginning of AD progression where connections between certain brain regions increase, it is believed that such localised changes are simply compensatory responses to the overall decline in brain connectivity. As the disease progresses, the exceptional trends reverts to the global pattern of a decrease in connectivity. Therefore, for the second endophenotype, we consider a matrix-valued phenotype representing brain-networks.

We denote the base vectorial phenotype of the kth ROI of individual i by y_ik. Although we generate these vectors from the ADNI data representing voxel-wise ventricular volume expansion rates, the simulation framework is not dependent on the chosen method; indeed any suitable method for generating realistic vectors representing alternative imaging phenotypes is entirely permissible.

We plot the linear regression slope coefficient of the ventricular volume against time and identify, for the 253 ADNI subjects, a 148023-voxel subset of the full-brain MRI scan which is maximally discriminative between AD and CN. Given that we aim to perform multiple RF runs over multiple scenarios, we perform a simple k-means clustering, transforming the 148023-voxel ventricular volume gradients to that of just 400 parcels, simulating a lower-resolution image. This choice is arbitrary except for the fact that it is larger than the sample size, thereby simulating a large-p/small-N scenario.

For each individual in our simulated population, the aim is simply to generate a 400-dimensional vectorial phenotype to represent the base ventricular gradients of the 400 ROIs on which the genetic model modifications can be imputed. We require, for realism, that after implementing the genetic modifications, the multivariate distribution of this phenotype resembles the empirical distribution of the subjects from the ADNI dataset. The practical implementation of the base ventricular volumes is performed in two steps, as follows.

For the first step, we generate random samples from a multivariate Gaussian distribution where the mean and covariance matrix parameters are taken from the sample mean and sample covariance matrix of the 154 CN subjects. Note that because there are fewer subjects than ROIs (154<400), the sample covariance matrix is unstable, requiring, therefore, the use of a James-Stein-type shrinkage estimator (Schäfer and Strimmer, 2005) for the covariance matrix. The distribution of rate gradients across all parcels of cognitive normal (CN) subjects has an average and standard deviation of –0.38 and 2.27 respectively, reflecting an overall decrease and variability in ventricular volume gradients in the base population across all the parcels. We denote the imaging trait of the kth ROI of individual i by y_ik.

For the second step, in anticipation of our genetic effects model where an overall increase in each vector component is observed, we modify the above random samples by reducing the gradient for selected ROIs. As shown in (Silver et al., 2012), the brain matter atrophy is not uniform throughout the brain. We simulate this effect by selecting a subset of n_d=28 randomly selected ROIs which we label disease-affected-ROIs, which corresponds to regions affected by selected SNPs which, in turn, has a measurable on the disease status. We label the ROI-index subset by R_d. Therefore the modification in this second step is

where ζ>0 is a parameter which is tuned according to the disease model (see below).

It is not unreasonable to assume that the volumes, and by extension the rate of volumetric changes, of functionally interconnected brain regions correlate positively (McAlonan et al., 2005). This link between structural and functional aspects can be used to indirectly estimate brain connectivity networks via the atrophy gradients. We simulate the brain connectivity by first considering the full (non-partial) correlations between the same 400 ROI. From the n_CN=154 subjects, we have the empirical covariance matrix

We generate the base covariance matrices for each of our 200 subjects by introducing random perturbations from the unit Normal at each component

with for subject i. We symmetrise the resulting matrix ensuring that it remains positive definite (Higham, 2002).

The graph representation of the network is inferred by identically matching the zero-components of the adjacency matrix of G_i with the zeros of the the sparse inverse covariance estimate (SICE) of (Sun et al., 2009; Huang et al., 2010). Let Θ⁽ⁱ⁾=(Σ⁽ⁱ⁾)^–1 be the inverse covariance matrix. Assuming that the observations follow a multivariate Gaussian, the SICE is obtained by maximising over Θ the L₁-penalised log-likelihood,

where ||Θ||₁ is the sum over the absolute values of all elements in Θ. We use the graphical lasso algorithm (Friedman et al., 2008) to estimate Σ, and hence Θ. The regularisation parameter is set as ρ=1 throughout.

One example of such a simulated brain network in the context on AD is found in Sun et al. (2009) where the authors used the SICE approach on FDG-PET data to uncover connectivity among different brain regions.

The brain connectivity network of each subject is represented by a 400-dimensional square matrix, representing either the covariance between the different ROIs or the adjacency matrix of a 400-node graph. In this paper we consider both cases. The covariance and adjacency matrices are written as and G_i respectively.

Genetic disease modelling

The second step in our simulations is the construction of the disease models. There are two aspects to this – the genetic effects on the endophenotypes, which includes both disease-related and “spurious” non-disease related effects, and the disease classification model.

Let y_β be the simulated baseline vectorial phenotype at the ROI with index β. We assume a simple additive genetic model where the components of the vector modified from the base value are given by

where the effect size w_β is built from multiplicative interaction effects as

P_‧ refers to the interaction causal models above, and x_a={0, 1, 2} is the minor allele frequency at locus a, etc. δ_β is an effect-size parameter which, when increased, with all else equal, monotonically boosts the population disease penetrance levels (see Disease classification models below). For simplicity, we fix δ_β=δ for all affected ROI.

In selected simulation runs, we include the presence of “spurious” SNP pairs and brain ROI (see Figure 2). These are labelled “spurious” because the genetic effects on the brain phenotype do not affect the propensity to develop the disease of interest. In our simulations, we assume that these do not overlap with their “causal” counterparts. Using the identical tuning parameter δ we employ a similar genetic effects model as above.

For the changes to brain networks, we propose an additive genetic model, similar to (24) where the strength of correlations between different ROI is scaled with a simple exponential factor. Let be the covariance matrix representing the brain connectivity. The components of the attenuated representative covariance matrix Σ′* are given by

where γ>0 is the effect-size parameter. In this model, both positive and negative correlations between different brain regions are reduced in magnitude.

The final piece in the imaging genetics simulation engine is the method for assigning a case-control disease status to each subject. This serves the following two purposes. It allows us to establish the power of imaging genetics via a comparison of quantitative vs. case-control RF regressions, and it is used to validate the supervised learning approach to define distances between brain networks.

To minimise the model complexity, we only consider the effects of the increase in vectorial phenotype on the disease classification, neglecting the contribution of the decrease in brain connectivity. Given that a correlation between structural and functional connectivity is built in into our simulations, this is not an unreasonable simplification.

Based on empirical observations of the ADNI data, we have the following design considerations for constructing an AD disease classification model. Firstly, such a model is inherently probabilistic with the propensity to develop AD increasing with the overall ventricular volume gradient of a sub-region of the whole brain. Analysis of the ADNI data in Silver et al. (2012) suggests that this highly discriminating region is made up of ~7% of the total brain area; this corresponds, in our simulation set-up of 400 parcels, to 28 ROI. Secondly, as already mentioned, the distribution of volume gradients in the output of such a model, for both AD and CN classes, must agree with the empirical conditional probability distributions from the ADNI data. Thirdly, the model should allow for a range AD population penetrance values which can be set independently of the effect size. We propose, therefore, the following simple, semi-empirical, Bayesian probability model for AD classification based on the average gradient of the causal region R_d.

Let the represent the average gradient of the causal region R_d for a given subject i, i.e.,

where n_d=28 in this simulation set-up, and the * symbolises the values post-genetic effects. Let the disease status be wholly determined by i.e., we have the following equivalence between conditional probabilities for a subject i begin given an AD classification

where z_i=1, 0 are the AD/CN classification label for subject i respectively. Using Bayes’s Theorem, we write

where p(z_i=1) is the AD population penetrance level. The likelihoods and are empirically determined from the ADNI data. As before, we assume a Gaussian distribution in the ADNI data, with the sample means and for the CN and AD classes respectively. To avoid the situation where the probability of a AD classification given a large, negative, is greater than the probability of a CN classification, and vice versa, we assume that the variances of the normal distributions, representing the respective likelihoods of the CN/AD classes, are identical. We therefore take the average of the two standard deviations, i.e., In Figure 13 we show the empirically determined AD classification probability (28) as a function of for several different AD population penetrance levels.

Figure 13

The empirical AD classification probability as a function of mean brain ventricular volume gradient. The individual curves (green, blue, red, and black) represent a range of AD population penetrance levels (5%, 20%, 35%, and 50% respectively).