Genome-wide association data classification and SNPs selection using two-stage quality-based Random Forests

Background Single-nucleotide polymorphisms (SNPs) selection and identification are the most important tasks in Genome-wide association data analysis. The problem is difficult because genome-wide association data is very high dimensional and a large portion of SNPs in the data is irrelevant to the disease. Advanced machine learning methods have been successfully used in Genome-wide association studies (GWAS) for identification of genetic variants that have relatively big effects in some common, complex diseases. Among them, the most successful one is Random Forests (RF). Despite of performing well in terms of prediction accuracy in some data sets with moderate size, RF still suffers from working in GWAS for selecting informative SNPs and building accurate prediction models. In this paper, we propose to use a new two-stage quality-based sampling method in random forests, named ts-RF, for SNP subspace selection for GWAS. The method first applies p-value assessment to find a cut-off point that separates informative and irrelevant SNPs in two groups. The informative SNPs group is further divided into two sub-groups: highly informative and weak informative SNPs. When sampling the SNP subspace for building trees for the forest, only those SNPs from the two sub-groups are taken into account. The feature subspaces always contain highly informative SNPs when used to split a node at a tree. Results This approach enables one to generate more accurate trees with a lower prediction error, meanwhile possibly avoiding overfitting. It allows one to detect interactions of multiple SNPs with the diseases, and to reduce the dimensionality and the amount of Genome-wide association data needed for learning the RF model. Extensive experiments on two genome-wide SNP data sets (Parkinson case-control data comprised of 408,803 SNPs and Alzheimer case-control data comprised of 380,157 SNPs) and 10 gene data sets have demonstrated that the proposed model significantly reduced prediction errors and outperformed most existing the-state-of-the-art random forests. The top 25 SNPs in Parkinson data set were identified by the proposed model including four interesting genes associated with neurological disorders. Conclusion The presented approach has shown to be effective in selecting informative sub-groups of SNPs potentially associated with diseases that traditional statistical approaches might fail. The new RF works well for the data where the number of case-control objects is much smaller than the number of SNPs, which is a typical problem in gene data and GWAS. Experiment results demonstrated the effectiveness of the proposed RF model that outperformed the state-of-the-art RFs, including Breiman's RF, GRRF and wsRF methods.


Background
The availability of high-throughput genotyping technologies has greatly advanced biomedical research, enabling us to detect genetic variations that are associated with the risk of diseases with much finer resolution than before. With genome-wide genotyping of single nucleotide polymorphisms (SNPs) in the human genome, it is possible to evaluate disease-associated SNPs for helping unravel the genetic basis of complex genetic diseases [1]. SNPs are single nucleotide variations of DNA base pairs, and it has been well established in the genomewide association studies (GWAS) field that SNP profiles characterize a variety of diseases [2]. In light of emerging research on GWAS, hundreds or thousands of objects (with disease or normal controls) are collected, each object is genotyped at up to millions of SNPs. This is a typical problem of the number of SNPs is typically thousands of times larger than the number of objects. The task is to identify genetic susceptibility of SNPs through assaying and analyzing SNPs at the genomewide scale [3].
A number of methods for analyzing of susceptibility of SNPs in GWAS have been proposed in the literature, where each SNP is analyzed individually [4]. However, it is found that only a small portion of the SNPs have main effects on the complex disease traits, but most of the SNPs have shown little penetrance individually. On the other hand, many common diseases in humans have been shown to be caused by complex interactions among multiple SNPs. This is known as multilocus interactions [5].
For dealing with the later challenge, one way of testing the interactions is to exhaustive search the interactions between all SNPs. The approach to test all two-SNPs to see how they are related to diseases is already quite time demanding [6]. Further exhaustive search for higher order interactions becomes computationally impractical because the number of tests increases exponentially with the order of interaction [7]. One approach to overcoming the drawbacks of the large computational cost using traditional statistical test is to first find a small set of high relevant SNPs using univariate tests on each SNP by discarding SNPs with high p-values, and then evaluate the interactions within the high relevant SNP subset [8].
This paper focuses on an approach for selecting informative SNPs, i.e. a small portion of the SNPs that has main effects on the disease, using Random Forests (RF) model [9]. RF has been applied successfully to genetic data in various studies [10][11][12][13][14]. A number of studies has used RFs to rank SNP predictors [15], to predict disease using SNPs [16] and to identify the effects related to diseases [17].
RF is an ensemble method for classification built from a set of decision trees that grow in randomly selected subspaces of data. Each tree is built using a bootstrap sample of objects. At each node, a random subspace of all SNPs is chosen to determine the best split to generate the children nodes. The size of subspace is referred to a parameter mtry that is used in growing the trees. Each node in a grown tree corresponds to a specific best predictor SNP in a subspace with mtry randomly selected SNPs. A grown tree in a forest is represented by a top-down decision tree, in which multiple decision paths from the root to different leaves go through the tree via various nodes. A decision path is a sequence of multiple SNPs including potential interactions between them in terms of hierarchial dependencies. In this way, RF can normally take into account interactions between SNPs (for details, see [18,19]).
A grown RF is able to yield a classification result and a measure of the feature importance for each SNP [18]. Although it is anticipated that RF will help to detect the SNP interactions, the task of selecting the relevant SNPs associated with complex disease in high dimensional genomewide data using RF method still poses significant challenges. In general, the SNP importance measure used to select the relevant SNPs is based on the impact of an SNP in predicting the response. The effectiveness of SNP importance depends on the performance of the grown RF that correctly classifies new objects of given SNPs.
A series of comprehensive studies revealed that the original RF implementation by Breiman is efficient to analyze low dimensional data. Bureau et al. [10] show that RF has worked well in a candidate gene case-control study involving only 42 SNPs. Lunetta et al. [19] show that RF can be applied to simulated data sets with no more than 1000 SNPs. However, it is computationally inefficient to build an accuracy RF model for high dimensional data. As a consequence, RF has rarely been applied on the genomewide level for SNP selection and classification. Specifically, the original RF implementation designed to use a small default SNP subspace size mtry, e.g., log 2 M + 1, is only suitable for low dimensional data, where M is the total number of SNPs. For high dimensional SNP data, there is usually a large number of SNPs that is considered to be irrelevant to the response, and only a small number of SNPs is relevant or informative. The simple random sampling method using a small mtry selects many subspaces without informative SNPs, and the number of objects is usually insufficient to generate numerous nodes to make it up to a good performance. To guarantee the performance of the generated RF model, previous studies recommended to use a relative large mtry in growing the trees of a RF when dealing with high dimensional data such as SNP data in GWA studies. However, the computational cost of the procedure of searching such a good mtry is very high which is dependent on the possible candidates to be searched. In the study by Schwarz et al. [20], a multiple sclerosis case control data set comprised of 325,807 SNPs in 3,362 individuals was used and it took 1 week to generate a full random forest on a server with 82 GHx CPU and 32 GB of memory, where the mtry values to search are 2 √ M , 2 √ M , 0.1M , 0.5M and M . It was found that RFs built by small mtry values for high dimensional SNP data had poor classification performance [21].
In this paper, we propose to use a new approach in learning RFs model using a two-stage quality-based SNP subspace selection method, which is specifically tailored for high dimensional data of GWA studies. The proposed R-F model is computationally efficient to analyze GWA data sets with thousands to millions of SNPs without the need of using a large value of mtry. Furthermore, it is able to deliver a better classification performance than the original RF implementation using large mtry with a large margin. Our idea is to first add shadow SNPs into the original GWA data set. The shadow SNPs do not have prediction power to the outcome. However, they can give an indicator for the selection of informative SNPs. We then apply a permutation procedure to this extended GWA data to produce importance scores for all SNPs. The p-value assessment is used to find a cut-off point that separates informative SNPs from the noisy ones. Any SNP whose importance score is greater than the maximum importance score of the shadow SNPs is considered as important. We then use some statistical measures to split the set of informative SNPs into two groups: highly informative SNPs and weak informative SNPs. When sampling an SNP subspace for building trees, we only select SNPs from these two groups. This maintains the randomness of RFs meanwhile assuring the selection of informative SNPs. The resulting RF model is able to achieve a lower prediction error and avoid overfitting.
We conduct a series of experiments on two genomewide SNP data sets (Parkinson disease case-control data set comprised of 408, 803 SNPs and Alzheimer casecontrol data set comprised of 380, 157 SNPs) to demonstrate the effectiveness of the proposed RF method. To validate the the conjecture that the approach is effective for problems with large M and small N , where N denotes the number of objects, we have conducted additional experiments on 10 other gene data sets with gene expression classification problems. Experimental results show that the proposed RF using two-stage qualitybased SNP sampling can generate better random forests with higher accuracy and lower errors than other existing random forests methods, including Breiman's RF, GRRF and wsRF methods.

Given a training data
where X i are predictor SNPs, Y ∈ Y ∈ {1, 2, . . . , c} is the outcome containing possible classes (diseases), N is the number of training samples (also called case-control objects) and M is the number of SNPs. Random Forests [9] independently and uniformly samples with replacement the training data L to draw a bootstrap data set L * k from which a decision tree T * k is grown. Repeating this process for K replicates produces K bootstrap data sets and K corresponding decision trees T * 1 , T * 2 , . . . , T * k which form a RF. Given an input X = x, letĥ k (x) denote the prediction of class j of input x ∈ R M by the kth tree, the RF prediction is obtained by aggregating the results given by all K decision trees, denoted asŶ, that iŝ where I(·) denotes the indicator function.

Importance score of SNP from a RF
The importance score of SNPs can be obtained in growing trees [9]. At each node t in a decision tree, the split is determined by the decrease in node impurity. The node impurity is the gini index. If a sub-data set in node t contains samples from c classes (c ≥ 2), the gini index is defined as Gini(t) = 1 c j=1p 2 j , wherep 2 j is the relative frequency of class in t. Gini(t) is minimized if the classes in t are skewed. After splitting t into two child nodes t 1 and t 2 with sample sizes N 1 (t) and N 2 (t), the gini index of the split data is defined as The SNP providing smallest Gini split (t) is chosen to split the node. The importance score of each SNP is computed over all K trees in a RF. These raw importance scores can be used to rank the SNPs.

Two-stage quality-based SNP sampling method for subspace selection
The importance scores from a RF only give a simple ranking of SNPs. However, it is very difficult to select informative SNPs because of the noisy nature of the GWA data. For better subspace selection at each node of a tree, we first need to distinguish informative SNPs from noisy ones. Then, the informative SNPs are divided into two groups based on the statistical measures. When sampling the SNP subspace, SNPs from these groups are taken into account. Since the subspace always contains highly informative SNPs which can guarantee a better split at any node of a tree.
In the first stage we build R random forests to obtain raw importance scores and then use Wilcoxon ranksum test [22] to compare the importance score of an SNP with the maximum importance scores of generated noisy SNPs called shadows. The shadow SNPs are added into the original GWA data set and they do not have prediction power to the outcome. From the replicates of shadow SNPs, we extracted the maximum value from each row of the importance score of the shadow SNP and put it into the comparison sample. For each SNP, we computed Wilcoxon test to check whether its mean importance score is greater than the maximum importance score of noisy SNPs. This test confirms that if a SNP is important, it consistently scores higher than the shadow over multiple permutations. Given a significance threshold θ (the default setting is 0.05), any SNP whose p-value is greater than θ is considered to be an irrelevant SNP and is removed from the system, otherwise, the relationship with the outcome is assessed. This method has been presented in [23].
In the second stage, we find the best subset of SNPs which is highly related to the outcome. We now only consider the subset of SNPsX obtained from L after neglecting all irrelevant SNPs and use a measure correlation function c 2 (X, Y) to test the association between the outcome label and each SNP X j . Let X s be the best subset of SNPs, we collect all SNPs X j whose p-value is smaller than or equal to 0.05 as a result from the c 2 statistical test. The remaining SNPs {X\Xs} are added into X w .
We independently sample SNPs from the two subsets and merge them together as the subspace SNPs for splitting the data at any node. The two subsets partition the set of informative SNPs in data without irrelevant SNPs. Given X s and X w , at each node, we randomly select mtry (mtry > 1) SNPs from each group of SNPs. For a given subspace size, we can choose proportions between highly informative SNPs and weak informative SNPs that depends on the size of the two groups. That is mtrys = [mtry × (||X s ||/||X||)] and mtry w = [mtry × (||X w ||/||X||)], where X s and X w are the number of SNPs in the groups of highly informative SNPs X s and weak informative SNPs X w , respectively. ||X|| is the number of informative SNPs in the input GWA data set. These are merged to form the SNP subspace for splitting the nodes in trees. This new sampling method always provides highly informative SNPs for the subspace at any node in growing a decision tree.
The RF algorithm using two-stage quality-based SNP sampling method We now present the random forest algorithm called ts-RF using a new SNP sampling method to generate splits at the nodes of CART trees [24]. The new algorithm is summarized as follows.
(1) Generate the extended data set of 2M dimensions by permuting the corresponding predictor SNP values for shadow SNPs.
(2) Build a random forest model RF from the extended data set and compute R replicates of raw importance scores of all SNPs and shadows with RF . Extract the maximum importance score of each replicate to form the comparison sample of R elements.
(3) For each SNP, take R importance scores and compute Wilcoxon test to get p-value.
(5) The c 2 statistical test is used to separate the highly and weak informative subsets of SNPs X s and X w , respectively.
(7) For each L k , grow a CART tree T k as follows: (a) At each node, select a subspace of mtry (mtry = mtry s +mtry w , mtry >1) SNPs randomly and separately from X s and X w and use the subspace SNPs as candidates for splitting the node.
(b) Each tree is grown nondeterministically, without pruning until the number of SNPs per leaf n min is reached.
(8) Given a X = x new , use Equation (1) to predict new samples on the test data set.

Evaluation measures
We used Breiman's method as described in [9] to calculate the average Strength (s), the average Correlation (r) and c/ s2 as performance measures of a random forest. Out-ofbag estimates were used to evaluate the strength and correlation. Given s andρ, the out-of bag estimate of the c/s2 measure can be computed with r/s 2 . The correlation measure indicates the independence of trees in a forest whereas the average strength correspond to the accuracy of individual trees. Low correlation and high average strength result in a reduction of general error bound measured by c/s2 which indicates a high accuracy RF model.
Let D t denote a test data set and N t denote the number of samples in D t . The two measures are also used to evaluate the prediction performance of the RF models on D t . One is the Area under the curve (AUC). The other one is the test accuracy, computed as: where I(·) is the indicator function, y i indicates the true class of x i ∈ D t and Q(x i , j) =

Results on SNPs data sets
We conducted experiments on two genome-wide SNP data sets whose characteristics are summarized in Table 1 "Abbr" column indicates the abbreviation of the genome-wide SNP data sets used in the experiments. The real data Alzheimer disease has been analyzed and reported in Webster et al. [25]. It contained genotypes of a total of 380,157 SNPs in 188 neurologically normal individuals (controls) and 176 Alzheimer disease patients (cases). The genotype data for Parkinson disease patients and controls were published in [26]. This genome-wide SNP consisted 271 controls and 270 patients with Parkinson disease, cerebrovascular disease, epilepsy, and amyotrophic lateral sclerosis. For raw genotype data with phs000089.v3.p2 study accession can be found in NCBI [1] dbGaP repository.
The 5-fold cross-validation was used to evaluate the prediction performance of the models on GWA data sets. From each fold, we built the models with 500 trees and the SNP partition was re-calculated on each training fold data set. We also compared the prediction performance of the ts-RF model with linear kernel SVM, taken from LibSVM [2], the values of regularization parameter by factors C were 2 -2 and 2 -5 , respectively. These optimal parameter C provided the highest validated accuracy on the training data set. The number of the minimum node size n min was 1. The parameters R, mtry and θ for pre-computation of the SNP partition were 30, 0.1M and 0.05, respectively. We used R to call the corresponding C/C++ functions from the ts-RF model and all experiments were conducted on the six 64bit Linux machines, each one equipped with IntelR XeonR CPU E5620 2.40 GHz, 16 cores, 4 MB cache, and 32 GB main memory. The ts-RF and wsRF models were implemented as multi-thread processes, while other models were run as single-thread processes. Table 2 shows the average of test accuracies and AUC of the models on the two GWA data asets using 5-fold crossvalidation. We compare our ts-RF model with the Breiman's RF method and two recent proposed random forests models, that are the guided regularized random forests GRRF model [27] and the weighting subspace random forests wsRF model [28]. In the GRRF model, the weights are calculated using RF to produce importance scores from the out-of-bag data, in which these weights are used to guide the feature selection process. They found that the least regularized subset selected by their random forests with minimal regularization ensures better accuracy than the complete feature set. Xu et al. proposed a novel random forests wsRF model by weighting the input features and then selecting features to ensure that each subspace always contains informative features. Their efficient RF algorithm can be used to classify multi-class data.
The latest RF [29] and GRRF [30] R-packages were used in R environment to conduct these experiments. For the GRRF model, we used a value of 0.1 for the coefficient g because GRRF(0.1) has shown competitive prediction performance in [27]. We can see that ts-RF and wsRF always produced good results with a different mtry value. The ws-RF model achieved higher prediction accuracy when using mtry = √ M . The ts-RF model using mtry = M p outperformed the RF, GRRF, wsRF models and SVM on both GWA data sets, where M p = ||X s || + ||X w || denotes the number of informative SNPs. The RF model requires a larger number of SNPs to achieve better prediction accuracy (mtry = 0.5M). With this size, the computational time for building a random forest is still too high, especially for GWA data sets. It can be seen that the ts-RF model can select good SNPs in the subspace to achieve the best prediction performance. These empirical results indicate that, when classifying GWA data sets with ts-RF built from small yet informative subspaces, the achieved results can be satisfactory.
[1] http://www.ncbi.nlm.nih.gov/ [2] http://www.csie.ntu.edu.tw/˜cjlin/libsvmtools Table 3 shows the prediction accuracy and Table 4 shows the c/s2 error bound of the random forest models with different numbers of trees while mtry = ⌊log 2 (M ) + 1⌋ was fixed on the GWA data sets, respectively. We conducted these experiments to compare the new model with  other random forests models and observed obvious improvement in classification accuracy on all GWA data sets. For the comparison of the c/s2 error bound, the GRRF model was not considered in this experiment because the RF model of Breimen's method [29] was used in the GRRF model as the classifier. The efficient wsRF model [28] and the Breimen's method were used for comparison in the experiment. We used the RF, wsRF and ts-RF models to generate random forests in different sizes from 20 trees to 200 trees and computed the average accuracy of the results from the 5-fold cross-validation. We can clearly see that the ts-RF model outperformed other models in classification accuracy and produced the lowest c/s2 error in most cases on all GWA data sets. The proposed ts-RF model was applied to the Parkinson genome-wide data and assigned a score of importance to each SNP. The resulting list of SNPs was investigated for potential relevance to the Parkinson disease. Table 5 shows the results of the top 25 SNPs that are located within gene regions studied by the previous work. For each SNP, details including the rank value, SNP ID, gene symbol, gene ID, and p-value obtained using Wilcoxon test. The boldface rows in the table are the interesting genes associated with Parkinson disease. The results of this real data analysis validate the findings of GWA studies such as PTPRD, EPHA4 and CAST. Results also give other potential SNPs and genes that may be associated with the Parkinson disease. Specifically, some of these SNPs were found not to be strongly associated with the Parkinson disease by traditional statistical tests because they have relatively high p-value. This provides evidence of the advantages of using the proposed ts-RF model to detect potential SNPs associated with the disease. However, interpreting results and assessing their biological plausibility is challenging. Biologists can perform further investigation to validate their relationship with the Parkinson disease.
In summary, ts-RF is a promising method for applying RF method to high-dimensional data such as GWA data. The application of ts-RF to GWA data may help to identify potential interesting SNPs that are difficult to be found with traditional statistical approaches.

Results on gene data sets
To validate our conjecture that the proposed ts-RF model is effective for GWA data, we have conducted additional experiments on gene data sets. In this experiment, we compared across a wide range the performances of the 10 gene data sets, used in [31,27]. The characteristics of these data sets are given in Table 6. Using this type of Table 3 The prediction test accuracy of the models on the SNP pair data sets against the number of trees K.   data sets makes sense, since the number of genes of these data sets are much larger than the number of patients. For the RF method to obtain high accuracy, it is critical to select good genes that can capture the characteristics of the data and avoid overfitting at the same time.
For the comparison of the models on gene data sets, we used the same settings as in [27]. For coefficient g we used value of 0.1, because GR-RF(0.1) has shown a competitive accuracy [27] when applied to the 10 gene data sets. From each of gene data sets two-thirds of the data were randomly selected for training. The other one-third of the data set was used to validate the models. The 100 models were generated with different seeds from each training data set and each model contained 1000 trees. The mtry and n min parameters were set to √ M and 1, respectively. The prediction performances of the 100 classification random forest models were evaluated using Equation (3). Table 7 shows the averages of the 100 repetitions of the c/s2 error bound when varying the number of genes per leaf n min . It can be seen that the RF, wsRF models produced lower error bound on the some data sets, for examples, COL, BR2, NCI and PRO. The ts-RF model produced the lowest c/s2 error bound on the remaining gene data sets on most cases. This implies that when the optimal parameters such as mtry = √ M and n min = 1 were used in growing trees, the number of genes in the subspace was not small and out-of-bag data was used in prediction, the results comparatively favored the ts-RF model. When the number of genes per leaf increased, so the depth of the trees was decreased, the ts-RF model obtained better results compared to other models on most cases, as shown in Table 7. These results demonstrated the reason that the two-stage quality-based feature sampling method for gene subspace selection can reduce the upper bound of the generalization error in random forests models. Figures 1, 2, 3, 4 show the effect of the two-stage quality-based feature sampling method on the strength measure of random forests. The 10 gene data sets were analyzed and results were compared to those of the random forests by Brieman's method and the wsRF model. In a random forest, the tree was grown from a bagging training data, the number of genes per leaf n min varied from 1 to 15. Out-of bag estimates were used to evaluate the strength measure. From these figures, we can observe that the wsRF model obtained higher strength on the two data sets NCI and BRA when the number of genes per leaf was 1. The strength measure of the ts-RF model was the second rank on these two data sets and it was the first rank on the remaining gene data sets, as shown in Figure 1. Figures 2, 3, 4 demonstrate the effect of the depth of the tree, the ts-RF model provided the best results when varying the number of genes per leaf. This phenomenon implies that at lower levels of the  Table 7 The (c/s2) error bound results of random forest models against the number of genes per leaf nmin on the ten gene data sets.    tree, the gain is reduced because of the effect of splits on different genes at higher levels of the tree. The other random forests models reduce the strength measure dramatically while the ts-RF model always is stable and produces the best results. The effect of the new sampling method is clearly demonstrated in this result. Table 8 shows the average test accuracy results (mean ± std-dev%) of the 100 random forest models computed according to Equation (3) on the gene data sets. The average number of genes selected by the ts-RF model, from 100 repetitions for each data set, are shown on the right of Table 8, divided into a strong group X s and a weak group X w . These genes were used by the two-stage quality-based feature sampling method in growing trees in ts-RF.
The results from the application of GRRF on the ten gene data sets were presented in [27]. From these  Table 8 Test accuracy results (mean ± std-dev%) of random forest models against the number of genes per leaf n min on the ten gene data sets.  Table 8, the GRRF model provided slightly better result on SRB data set in case n min = 15 and PRO in case n min = 1, respectively. The wsRF model presented the best result on BRA and NCI data sets in case n min = 15. In the remaining cases on all gene data sets, the ts-RF model shows the best results. In some cases where ts-RF did not obtain the best results, the differences from the best results were minor. This was because the two-stage quality-based feature sampling was used in generating trees in the ts-RF, the gene subspace provided enough highly informative genes at any levels of the decision tree. The effect of the two-stage quality-based feature sampling is clearly demonstrated in these results.

Conclusion
We have presented a two-stage quality-based random forests for genome-wide association data classification and SNPs selection. The presented approach has shown to be effective in selecting informative sub-groups of SNPs and potentially associated with diseases that traditional statistical approach might fail. The proposed random forests model works well for the data where the number of case-control objects is much smaller than the number of SNPs, which is a typical problem in GWAS.
We have conducted a series of experiments on the two genome-wide SNP and ten gene data sets to demonstrate the effectiveness of the proposed RF model. The top 25 SNPs in Parkinson data set were identified by the proposed RF model including some interesting genes associated with neurological disorders. Experimental results have shown the improvement in increasing test accuracy for GWA classification problems and reduction of the c/s2 error in comparison with other state-of-the-art random forests, including Breiman's RF, GRRF and wsRF methods.