Ensemble learning for detecting gene-gene interactions in colorectal cancer

CRC GWAS data and preprocessing

The colorectal (CRC) GWAS case-control data were collected from CRC patients and healthy individuals with matching age, gender, and geographical distributions within the province of Newfoundland and Labrador, Canada. The CRC Transdisciplinary (CORECT) consortium coordinated the genotyping of data. Genotyping was conducted using a custom Affymetrix genome-wide platform (the Axiom CORECT Set) on two physical genotyping chips (pegs) for two datasets with around 1.2 and 1.1 million SNPs (Schumacher et al., 2015). The first dataset has 1,236,084 SNPs and 696 samples with 200 cases and 496 controls and the genotyping rate was 0.997. In GWAS, genotyping rate is computed as the percentage of samples (including both cases and controls) that are successfully genotyped. The second dataset has 1,134,514 SNPs and 656 cases with a genotyping rate of 0.888. Using PLINK (Purcell et al., 2007), a tool for analyzing genetic data, we merged these two datasets based on their common SNPs and obtained a dataset of 265,195 SNPs and 1,152 unique samples. Among the samples, 656 were cases and 496 were controls.

Next, we performed per-sample and per-marker quality control steps and linkage disequilibrium (LD) pruning on the CRC data using PLINK. In the per-sample quality control, sex check was performed, sex chromosomes were excluded, and samples with more than 1% missing genotypes and outlier heterozygosity rate were removed. In the per-marker quality control, SNPs with missing call rates higher than 5%, with minor allele frequencies (MAFs) less than 5%, or with Hardy-Weinberg equilibrium values greater than 0.0001 were removed. Then, we pruned SNPs that were dependent on each other (correlation coefficient r² > 0.6, LD window size 2,000). Missing genotypes were imputed using the most frequent alleles in the population. Last, the genotypes were re-coded using numerical values {0, 1, 2} with 0 standing for homozygous reference, 1 for heterozygous variant, and 2 for homozygous variant. After the steps of quality control and imputation, the final preprocessed and balanced dataset had 186,251 SNPs and 944 samples (472 being cases and 472 being controls).

The original data collection was approved by Memorial University Health Research Ethics Authority (HREA) with the approval number HIC 01.70. Our study was recognized by HREA as the use of secondary data which have already been collected and de-identified, and did not require a clearance.

Feature selection

Most computational methods in informatics find it prohibiting to analyze high dimensional GWAS data (Moore, Asselbergs & Williams, 2010) due to the massive number of attributes, that is, hundreds of thousands of SNPs in the data. It is impossible for machine learning methods to detect interactions by enumerating all combinations of SNPs in a GWAS dataset. Furthermore, the existence of redundant and irrelevant attributes hinders machine learning methods to reveal actual gene-gene interactions in the data. These together set the stage for the necessity of dimensionality reduction or feature (attribute) selection for analyzing GWAS data (Moore, Asselbergs & Williams, 2010).

Feature selection is frequently used as a data filtering step in machine learning when the original data contain noisy or irrelevant features, or attributes, that could compromise the prediction power of learning algorithms (Yu & Liu, 2003). Feature selection methods choose only a subset of the most important features, and thus reduce the dimensionality of the data, speed up the learning process, simplify the learned model, and improve the prediction performance (Dash & Liu, 1997; Guyon & Elisseeff, 2003).

In our previous study (Dorani & Hu, 2018), six feature selection algorithms, including chi-square, logistic regression, odds ratio, and three Relief-based algorithms (ReliefF, Tuned ReliefF (TuRF), and Spatial Uniform ReliefF) were compared based on how they rank the most important SNP attributes that contribute to the phenotypic outcome through gene-gene interactions. We applied the feature selection methods to both simulated and real GWAS datasets and showed that Relief-based methods, specifically TuRF, performed the best in filtering SNPs associated with a disease through interaction effects. This observation was also supported by other independent studies (Beretta & Santaniello, 2011; Urbanowicz et al., 2018). Therefore, in this study, we used the TuRF feature selection method to reduce the CRC GWAS dataset to be of a manageable size such that it can be analyzed using the downstream classification algorithms in a reasonable time frame.

Ensemble learning

Many preliminary GWAS employed uni-variable approaches where interactions between multiple variables could be overlooked (Moore & Ritchie, 2004). Moreover, the parametric linear statistical models look at multiple variables but have limitations for detecting non-linear interactions (Moore, Asselbergs & Williams, 2010). Thanks to the intrinsic multi-variable and non-linear properties, tree-based ensemble learning methods have been proved to be a powerful analysis tool for detecting interacting genes in GWAS (Moore, Asselbergs & Williams, 2010; Upstill-Goddard et al., 2012). They can be used to train highly accurate classifiers, as well as to discover new genetic markers by ranking genetic variables based on their importance in classification.

Ensemble learning methods use an aggregation of predictors known as base learners. To produce a final prediction, the predictions of the base learners are weighted and the overall predictions are decided through majority voting for classifications and averaging for regressions. It has been shown that ensemble learning methods are powerful at reducing variance and overfitting by utilizing a collection of diverse base learners such as classification and regression trees (CART) (Breiman, 1996; Dietterich, 2000). In this article, we adopted two tree-based ensemble learning algorithms, random forests and GBM, to explore their power in detecting non-linear gene-gene interactions for GWAS.

Random forests

One popular ensemble learning method is bagging (short for bootstrap aggregating) that uses bootstrapped samples of the training data to train classification or regression models separately (Breiman, 1996). Bagging reduces the variance of an estimated prediction function, and works especially well for high-variance, low-bias procedures, such as trees (Friedman, Hastie & Tibshirani, 2001). The random forests algorithm is a special case of bagging where the variables are randomly selected to determine the optimal split at each node of the tree (Breiman, 2001). Random forests are shown to be a very powerful regression and classification method which utilizes a large collection of possibly uncorrelated decision trees (Breiman, 2001; Szymczak et al., 2009; Ziegler, DeStefano & König, 2007; Schwarz, Inke & Ziegler, 2010). Each tree is grown using the CART methodology (Breiman, 1996).

In random forests, ntree trees are grown independently using bootstrapped training samples. While constructing an individual tree, for each node, mtry (<M) predictor variables are randomly picked from the total M variables in the original data. Then the best variable/partition-point among the mtry variables is picked and used to split the node into two daughter nodes. Each tree grows to its maximum size, that is, when there is only one sample in the leaves. To make a classification for a new testing sample, each tree casts a vote for the predicted class and the majority vote will be the final prediction result. Intuitively, reducing mtry will reduce the correlation between any pair of trees in the ensemble. For classification, the default value of mtry is $\lfloor \sqrt{M}\rfloor$; however, the best value for this parameter will depend on the application problem, and it should be treated as a tuning parameter (Breiman, 2001; Friedman, Hastie & Tibshirani, 2001).

The random forests algorithm is effective in uncovering interactions among genes that do not exhibit strong main effects (Moore, Asselbergs & Williams, 2010). The algorithm has been employed in various studies to predict rheumatoid arthritis risks (Sun et al., 2007), to rank SNP predictors (Schwarz et al., 2007; Sun et al., 2008), and to detect the epistatic effects associated with human diseases (Garca-Magariños et al., 2009; Pan et al., 2014).

Goldstein et al. (2010) indicated that using the random forests algorithm with default settings of hyper-parameters would not yield optimal results for large GWAS datasets. In contrast, tuning the hyper-parameters, including mtry (number of random variables to make best split at each tree node) and ntree (number of trees), and using greater values, work well generally for large GWAS datasets.

In this study, we used a very fast implementation of random forests provided in an R package called “ranger” (Wright & Ziegler, 2017). The “ranger” package provides all functionalities similar to the RandomForest package in R with much greater speed. Therefore, we can use it for the GWAS datasets with a large number of SNPs. We performed parameter tuning on mtry and ntree in our implementation. For mtry we selected the values of {100, 200, 300, 500, 1,000}, and the selected values for ntree were {500, 1,000, 2,000}. Other parameters used the package default settings. Therefore, we had 15 different combinations of mtry and ntree.

In addition, we performed a 10-fold CV. The original dataset was partitioned into 10 equal subsets, and for each of the 10 iterations, one subset was picked as the testing data and the rest nine subsets served as the training data. Thus, each sample was used in the testing data exactly once, and the CV accuracy was computed as the percentage of correctly classifying all the testing samples. For each parameter combination and each fold of the CV (150 iterations in total), we repeated the random forests algorithm 10 times given the stochastic nature of bootstrapping and selecting variables to grow the trees. Since each execution of the random forests will yield a unique predictive model, we computed the accuracy of classifying each testing sample as the fraction of it being correctly predicted by the 10 models resulted by the 10 executions of the algorithm. Then the average CV accuracy using a specific combination of parameters was computed as the averaged accuracy across all testing samples. The parameter combination with the best average CV accuracy will be picked.

When the best combination of parameters was picked, the predictive models learned using it were then used to rank features according to the importance of the features in classification. For random forests, we used Gini Index as the quantification of feature importance (Breiman, 2017). Gini Index measures the inequality of dividing samples of two classes using a feature at a node of the trees. A higher Gini Index of a feature infers a better ability of using this feature to differentiate the two classes of the samples. Since we had multiple learned predictive models as a result of CV and repetitive runs, the final importance score of a feature was computed as the average Gini Index across all models.

Gradient boosting machine

Another method to generate an ensemble is boosting where the base learners evolve over time and make weighted votes (Friedman, Hastie & Tibshirani, 2001). In a boosting algorithm, many base learners are built and the new learners improve on the previous ones. The learners are trained sequentially, which result in building a “committee” of complex predictors (Friedman, Hastie & Tibshirani, 2001; Friedman, 2001).

Gradient boosting machine is a boosting machine learning algorithm in which a weighted combination of predictors are used to make the final prediction (Friedman, 2001). A set of n.trees base learner trees are created in an iterative fashion, where a new one learns from the previous one. Let F_m(x), where m = 1, 2,…,n.trees, denote the approximation that maps the input variables x to the desired output y. At iteration m, a new approximation F_m(x) is constructed through improving on the previous one F_m−1(x) by adjusting it using the gradient of the loss function ∇L(y, F_m−1(x)).

For implementation, we used an R package called “gbm” (Ridgeway, 2007). Similar to random forests, we chose a range of values for each hyper-parameter in order to search for the best combination. GBM has three main parameters, n.trees (the number of trees), interaction.depth (the complexity of interactions between nodes, i.e., features), and shrinkage (the learning rate or step-size reduction). The testing values of these parameters were as follows: n.trees = {100, 500, 1,000, 2,000}, interaction.depth = {1, 2, 10}, and shrinkage = {0.001, 0.01, 0.1}. Therefore, we had 36 different combinations of the three parameters. Other parameters of GBM, such as n.minobsinnode (minimum number of samples in the tree terminal nodes), bag.fraction (the fraction of the training set samples randomly selected to construct a new tree), and train.fraction (the fraction of samples used to fit a GBM model), were set to default values as 10, 0.5, and 0.5, respectively.

Similar to random forests implementation, we used a 10-fold CV and repeated the algorithm 10 times for each combination of parameters and each fold of the CV. The best combination of parameters found was then used to rank the features. Feature Importance was estimated based on the number of times a feature is selected for splitting and the improvement to the model as a result of each split, weighted by the number of observations each split operates, and averaged over all trees (Friedman, Hastie & Tibshirani, 2001). Again, since we repeated GBM algorithm 10 times for each iteration of the 10-fold CV, the final importance score of a feature was averaged over 10 runs and 10-folds.

Statistical interaction analysis

Both the random forests and GBM algorithms provide importance estimations on each SNP based on their contribution to the accurate classification. Such importance scores by the ensemble algorithms help to identify the top-ranked SNPs that are most possibly associated with the disease of CRC. It has been suggested that tree-based ensemble learning algorithms are able to detect interaction effects among multiple SNP features. Therefore, we employ an information gain method (Hu et al., 2011; Hu et al., 2013a, 2013b) to quantitatively investigate the statistical evidence of pairwise and three-way interactions among the top identified SNPs.

The information gain method considers the genotype of each SNP and the disease outcome as random variables. It first computes the mutual information I(A; D) between SNP A and the disease outcome D as the main effect of the individual SNP A, since it captures how much information of the disease status the genotype of SNP A provides. Then the information gain IG(A; B; D) of combining two SNPs A and B on explaining the disease status D is computed by subtracting the main effects of A and B, that is, I(A; D) and I(B; D) from the total mutual information I(A; B; D) by joining A and B to explain D. Essentially, the information gain IG(A; B; D) captures the gained information of combining A and B on explaining D, that is, the synergistic and non-additive interaction effect between SNPs A and B.

Similarly, when three SNPs, A, B, and C, are considered, the information gain IG(A; B; C; D) is computed by subtracting all the main effects and pairwise interaction effects from the total mutual information I(A; B; C; D) between joining A, B, and C together and the disease status D. Such an information gain metric measures the pure three-way synergistic interaction among three SNPs on explaining the disease outcome.

Permutation testing can be used to assess the significance level of a computed information gain value. In each permutation, the disease status labels of the data are randomly shuffled, and the information gain values are computed again for the permuted data. In this study, we performed a 1,000-fold permutation testing, that is, 1,000 permutations were collected to provide a null distribution of the assumption that there was no association between the genotypes and the disease outcome. Then the p-value was assessed for each computed information gain by counting the portion of permuted datasets that had a greater value of the information gain than the real data.

Functional enrichment analysis

We also annotated the identified top SNPs with functional information using online resources including the ENSEMBL (http://www.ensembl.org) and the National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/) databases. The databases provide biological information on the allele, chromosome, and gene information for each SNP.

Then we used the Database for Annotation, Visualization, and Integrated Discovery (DAVID) bioinformatics tool (Dennis et al., 2003; Huang, Sherman & Lempicki, 2008) to perform a functional enrichment analysis using the genes from the identified top SNPs. Functional enrichment analysis is a method to identify gene functional categories or associations with diseases that are over-represented in a set of genes. We compared the set of identified genes to the background of the entire human genome, in order to test which functional categories were significantly over-represented. The enriched functional categories may help better understand the disease of CRC. DAVID is widely used for enrichment analysis, specifically for exploring the functions of genes. It is a web-based tool which takes a list of genes as an input and produces annotation tables and charts to show diseases and relevant gene ontology (GO) terms enriched by a given list of genes.

Filtered SNPs using TuRF

We ran the TuRF feature selection method on all 186,251 SNPs. Each SNP was assigned a score based on its contribution to the disease status. Since TuRF is a multi-variable feature selection algorithm, such a contribution may include both the individual main effect or the effect of interacting with other SNPs. Fig. 1 shows the distribution of the TuRF scores for all the SNPs.

To reduce the computational overhead of the downstream ensemble learning analysis, we then picked a score threshold and only included the potential SNPs with TuRF importance scores higher than the threshold. The threshold mean + 3SD was chosen which filtered the top 2,798 SNPs. The feature subset was an appropriate size such that the training of the ensemble learning algorithms could be finished in a reasonable time frame, meanwhile, enough features were included in the subsequent analysis to identify possible gene-gene interactions.

Optimized parameters for the ensemble learning algorithms

The parameter tuning for the random forests algorithm included adjusting the values of two parameters, mtry = {100, 200, 300, 500, 1,000} and ntree = {500, 1,000, 2,000}. Thus, there were 15 combinations of parameter values. Recall that we ran the algorithm 100 times for each combination, the average CV accuracy and area under the curve (AUC) were computed over 100 runs.

Figure 2 shows the comparative results for the 15 configurations of the random forests algorithm. The average CV accuracy and AUC are shown as a function of the parameter ntree, and different curves represent using different values of the other parameter mtry. The highest CV accuracy of 75% and AUC of 0.84 were achieved when mtry = 100 and ntree = 2,000. That is, the algorithm performed the best when ntree was set to the maximum value and mtry was set to the minimum value. This suggests that greater values of ntree and lower values of mtry are preferable for analyzing GWAS data. Both increasing the number of trees and reducing the number of randomly picked variables to construct each node of a tree enhance the diversity of trees in the ensemble. The result suggests that the genetic explanation of complex diseases may be highly heterogeneous, that is, a large number of diverse classification models that use different combinations of SNP attributes are needed to approximate the mapping from the genotypes to the disease outcome.

The GBM algorithm had three parameters, that is, n.trees = {100, 500, 1,000, 2,000}, interaction.depth = {1, 2, 10}, and shrinkage = {0.001, 0.01, 0.1}. Again, for each of the 36 combinations of parameter combinations, 10 runs were collected for each of the 10-fold CV, and the average CV accuracy and AUC were computed over 100 runs.

Figure 3 shows the comparative results for the parameter tuning. Each sub-figure uses a shrinkage value, and shows the average CV accuracy (left) and AUC (right) as a function of the parameter n.trees. The highest CV accuracy was 74% and the CV AUC was 0.82, achieved when n.trees was set to 2,000, interaction.depth was 10, and shrinkage was 0.1. This suggests that, for GWAS data analysis, GBM performs better when greater values of interaction.depth, n.trees, and shrinkage are used.

Top-ranked SNPs

After parameter tuning, we chose the optimal parameter combinations of the random forests algorithm, that is, mtry = 100 and ntree = 2,000, and of the GBM algorithm, that is, n.trees = 2,000, interaction.depth = 10, and shrinkage = 0.1, to assess the importance of the SNP features. In random forests, the Gini Index importance of a feature was estimated based on the impurity of splitting testing samples on the tree node using the feature. Since we repeated the algorithm 100 times (10 times for each fold of the 10-fold CV), the final importance score (between 0 and 1) of a feature was computed as the averaged impurity score over 100 runs. The averaged importance score of each feature using GBM was computed similarly.

Figure 4 shows the scores of the 2,798 SNPs assessed by both algorithms. We see that the scorings by random forests and GBM are highly correlated, which indicates the consistency of the two algorithms. Since a higher score indicates a higher importance of a SNP in classification, the top-right corner of the figure presents the SNPs ranked high by both the random forests and the GBM algorithms. The bottom-left corner of the figure includes a dense cluster of SNPs with low importance scores by both algorithms. We draw a line to separate SNPs into two clusters. The 44 SNPs on the right side of the separating line are the most important SNPs identified by both algorithms, and will be investigated subsequently on their biological functions and statistical interactions. Among them, SNP rs3760948 from gene ARRDC5 is ranked as the most important by both algorithms.

These 44 SNPs and their chromosome and gene information are listed in Table 1. Gene information of 29 SNPs was found using the ENSEMBL and the NCBI databases. In the table, 15 SNPs without any gene information are from non-coding DNA regions. The p-values of SNPs and their ranking among all 186,251 SNPs are also shown in the table. Figure 5 further shows the distribution of the 44 identified SNPs in each chromosome. The comparison of coding vs. non-coding regions is also shown. There were no identified SNPs in chromosomes {10, 13, 15, 16, 17, 22}. Chromosome eight contains six identified SNPs in which four of them belong to coding regions. Chromosome two has four out of five identified SNPs from the coding regions of DNA.

Statistical interactions among the identified SNPs

We calculated the information gain between all $\left(\begin{array}{l}44\\ 2\end{array}\right)$ pairs of the 44 identified SNPs. A 1,000-fold permutation test was performed to assess the significance level of each pair. Table 2 lists the 17 most significant SNP pairs (p < 0.02) with their information gain measures and the p-values. The maximum information gain was 1.3% for the interaction between SNPs rs2010907 and rs3760948 (p = 0.002), meaning that the synergistic interaction effect between the two SNPs explained 1.3% of the disease outcome. Six pairs of coding SNPs, that is, SNPs from genes, had significant interaction effects.

We then calculated the three-way information gain among all $\left(\begin{array}{l}44\\ 3\end{array}\right)$ possible trios of the 44 identified SNPs. Table 3 lists the 16 most significant (p ≤ 0.001) SNP trios with their information gain measures and p-values. The strongest three-interaction was found among SNPs from genes OTOP1, EPHX2, and LINC with a strength of 2.55% and a significance level of p < 0.001, suggesting that the pure three-way synergistic effect among the three SNPs can explain 2.55% of the disease outcome. There were six trios of significant interacting SNPs all from coding regions of DNA.

Enriched gene functional terms

We submitted the identified 29 genes to the DAVID software, and chose categories of Disease, GO, and Pathway for a functional enrichment analysis. We set the gene count threshold as five, that is, only functional categories including more than 4 of the 29 identified genes were considered, and the significance cutoff as 0.05, that is, Fisher’s exact test p < 0.05. Table 4 lists the most significantly enriched functional categories. There were 16 enriched terms, four of which were diseases, and the rest were GO terms. The most significantly enriched term was the disease tobacco use disorder with 14 out of 29 genes in the category and a significance level of 3.9 × 10⁻⁵.