Volume 11 Supplement 8
Proceedings of the Neural Information Processing Systems (NIPS) Workshop on Machine Learning in Computational Biology (MLCB)
A machine learning pipeline for quantitative phenotype prediction from genotype data
 Giorgio Guzzetta^{1, 2},
 Giuseppe Jurman^{1} and
 Cesare Furlanello^{1}Email author
DOI: 10.1186/1471210511S8S3
© Furlanello et al; licensee BioMed Central Ltd. 2010
Published: 26 October 2010
Abstract
Background
Quantitative phenotypes emerge everywhere in systems biology and biomedicine due to a direct interest for quantitative traits, or to high individual variability that makes hard or impossible to classify samples into distinct categories, often the case with complex common diseases. Machine learning approaches to genotypephenotype mapping may significantly improve GenomeWide Association Studies (GWAS) results by explicitly focusing on predictivity and optimal feature selection in a multivariate setting. It is however essential that stringent and well documented Data Analysis Protocols (DAP) are used to control sources of variability and ensure reproducibility of results. We present a genometophenotype pipeline of machine learning modules for quantitative phenotype prediction. The pipeline can be applied for the direct use of wholegenome information in functional studies. As a realistic example, the problem of fitting complex phenotypic traits in heterogeneous stock mice from single nucleotide polymorphims (SNPs) is here considered.
Methods
The core element in the pipeline is the L1L2 regularization method based on the naïve elastic net. The method gives at the same time a regression model and a dimensionality reduction procedure suitable for correlated features. Model and SNP markers are selected through a DAP originally developed in the MAQCII collaborative initiative of the U.S. FDA for the identification of clinical biomarkers from microarray data. The L1L2 approach is compared with standard Support Vector Regression (SVR) and with Recursive Jump Monte Carlo Markov Chain (MCMC). Algebraic indicators of stability of partial lists are used for model selection; the final panel of markers is obtained by a procedure at the chromosome scale, termed ’saturation’, to recover SNPs in Linkage Disequilibrium with those selected.
Results
With respect to both MCMC and SVR, comparable accuracies are obtained by the L1L2 pipeline. Good agreement is also found between SNPs selected by the L1L2 algorithms and candidate loci previously identified by a standard GWAS. The combination of L1L2based feature selection with a saturation procedure tackles the issue of neglecting highly correlated features that affects many feature selection algorithms.
Conclusions
The L1L2 pipeline has proven effective in terms of marker selection and prediction accuracy. This study indicates that machine learning techniques may support quantitative phenotype prediction, provided that adequate DAPs are employed to control bias in model selection.
Background
Fitting quantitative phenotypes from genomewide data is a rapidly emerging research area, also object of dedicated data contests [1–3]. Given the complexity of the molecular mechanisms underlying many common human diseases, one of the most significant challenges to catch genetic variations associated to functional effects is enabling a modeling approach that is really multivariate and predictive [4]. In particular, it is clear that modeling should be based on patterns of multiple SNPs (with patterns’ structure extending the notion of haplotype) rather than on single SNPs. Attention is thus directed towards machine learning methods that can provide SNP selection simultaneously with the regression model, and manage highorder interactions and correlation effects among features. In this view, a handy offtheshelf solution is the application of the Random Forest method [5], available with fast implementations (e.g. RandomJungle: http://www.randomjungle.org) both for classification (casecontrol studies) or regression (quantitative phenotype fitting). Regarding the haplotype data pattern problem, new kernel functions have been proposed for predictive classification by Support Vector Machines (SVM) in a crossvalidation experimental framework [6].
Given that flexible machine learning methods for genotype data are becoming available, the second top challenge is building around the modeling exercise a framework that controls the sources of variability involved in the process. Lack of reproducibility in GWAS has been investigated and is known to have multiple causes [7]. Some of the technical causes may well transfer to genotype analyses by multivariate machine learning. Specifically, it is critical to consider the risk of selection bias [8, 9] to warrant that predictive values and molecular markers be reproducible across studies on massive genotype datasets. The issue of reproducibility regards the whole sequence of preparatory and preprocessing steps (upstream analysis), model selection, application and validation (downstream analysis).
Baggerly and Coombes [10] proposed a “forensic bioinformatics” approach to revise a highlyinfluential series of medical papers on genomic signatures predicting response to chemotherapeutic agents. Their attempt at reproduction of the original results led to the discovery of a series of fatal flaws on data preparation and application of methods to publiclyavailable microarray and preclinical chemosensitivity data for several cancer cell lines. A series of clinical trials has been suspended as a consequence. For machine learning methods, the stage of model selection is usually the most complex. To overcome variability and bias effects arising from choices hidden in the modeling path, a serious effort has been provided by the FDA’s led initiatives MAQC and MAQCII [11]. In particular, for classifiers of microarray data, the MAQCII consortium has studied how predictivity and stability of biomarkers is associated to the type of adopted Data Analysis Protocol (DAP), intended as a standardized description of all steps in training, model selection and validation on novel data [12]. The type of internal and external validation methods used for selection of the best markers and models results as one of the main effects on predictive accuracy. Interactive effects of choices in the analysis design (e.g. batch size and composition) have been demonstrated also in GWAS in an extension of the MAQCII study [13]. However, limited efforts have been directed to detailed DAPs in the regression framework, and on genotype data in particular. In this work we propose a machine learning regression approach for genometophenotype prediction to improve the use of quantitative phenotypes as target variables in functional genomics. We consider first a standard Support Vector Regression (SVR) algorithm and then the L1L2 regression [14] approach. The machine learning methods are part of a software pipeline that implements a complete DAP for regression on genotype data. The L1L2 pipeline also includes a model selection module based on the concept of stability of ranked lists, previously developed for genomic profiling [15]. A procedure testing for markers highly correlated and proximal on the chromosome, termed saturation, is also provided in the pipeline. We present examples of prediction of quantitative phenotypes on a genomewide dataset of 12K SNPs. The dataset used in this study , which we will refer to as the ”GSCAN dataset”, is publicy available (website: http://gscan.well.ox.ac.uk), courtesy of the Wellcome Trust Center for Human Genetics. Data include familiar, genotype and phenotype information from a population of 4 generations of heterogeneous stock mice [16]. Two quantitative phenotypes were used: the percentage of CD8+ cells (CD8+), and the Mean Cell Haemoglobin (MCH). The number of samples is 1521 for %CD8+ and 1591 for MCH. The results from our methods are compared with those of a Reversible Jump Monte Carlo Markov Chain (MCMC) model adapted to fitting quantitative phenotypes and applied to the same dataset [1].
Results and discussion
Prediction accuracy on GSCAN mice data
Method  CD8+  MCH 

SVR  0.306 (0.2800.333)  0.147 (0.1250.169) 
L1L2  0.316 (0.2830.347)  0.106 (0.0950.116) 
MCMC  0.314  0.111 
Topcorrelated SNPs characteristics
SNP name  chromosome  n  mean distance (bp)  min distance  max distance 

mhcCD8a2  6  16  724067.5  2701  1735204 
rs13478736  6  13  1372896.5  187157  4078235 
mCV24938952  8  12  2824059.5  191536  6126682 
rs6375522  1  9  952187.1  115037  2360403 
rs13482427  15  9  773579.3  244547  1555507 
rs3145663  17  8  1184576.6  293271  2766351 
rs3672987  17  7  3346835.3  1194418  5686326 
rs13476229  1  6  813948.7  328290  1088152 
rs3684143  17  5  292337.2  116575  505713 
rs3678696  17  5  452020.4  77709  742231 
Conclusions
Prediction of quantitative phenotypes from highthroughput genotype data is an emerging research goal with significant applications. It can be envisioned that this predictive modeling problem will evolve into fitting a multidimensional phenotype pattern or a phenotype trajectory. More sophisticated predictive tools still need to be developed to achieve this goal, but it is anyway urgent to deploy experimental setups that can appropriately support model selection and biomarker identification. Here we introduced a framework for the systematic use of machine learning regression methods on wholegenome datasets. Building on results from the FDA’s led MAQCII initiative, the framework includes a pipeline of procedures (defined through a DAP) to avoid selection bias and ensure reproducibility. The application of the pipeline to up to 550 000 features was made feasible by an efficient software implementation, also suitable for high performance computing facilities. A DAP reproducing those of the original study [1] and employing a Gaussian kernel SVR obtained results comparable with the MCMC method. However the model selection solution does not protect from overfitting and cannot directly derive a list of selected SNP. The L1L2 method was as accurate as the reference study despite the use of a more stringent DAP, and it is able to provide an embedded feature selection which has shown to cope well with the problem of recovering correlated variables. An adjuvant therapy to the issue of correlated variables was proposed with the SNP saturation procedure, based on the concept of topcorrelated features. The saturation procedure can be seen as a black box algorithm within the pipeline that automates an analysis by Linkage Disequilibrium after one biomarker is found. SNP saturation also reduces spatial sparsity, because the additional markers are in general close to the topranked markers, as shown on the GSCAN data. The finding opens the possibility of encoding by special kernels feature and spatial correlation together.
The L1L2 pipeline also makes use of a model selection criterium aimed at increasing the stability of the list of candidate markers. As a result, the features selected by L1L2 are compatible with results of a previous GWAS on the same dataset [16]. This study confirms that machine learning approaches may support a more effective and reproducible use of multivariate genotype data for the prediction of quantitative traits [17].
Methods
Machine Learning methods
The solution (naïve elastic net) correctly selects the relevant features, but with biased weights. Zou & Hastie [20] corrected by rescaling the weights. In the approach proposed in [14] and used in this study, the correction is done by a Regularized Least Squares (RLS) regression performed only on the subset of features selected after the optimization of w in equation 1. The optimal weights in the RLS regression are found as
where and refer to the input data and the regression weights restricted to the subset of selected features. Thus, µ and τ modulate the feature selection, whereas the regularization parameter λ of the RLS controls the weight bias. The minimization of Eq. 1 is computed with a modified gradient descent algorithm, which makes use of weight values shrinkage through an iterative thresholding algorithm. For a more detailed description of the method, we refer the reader to the original paper [14]. L1L2 has been applied to define a transcriptomic profile of hypoxia in neuroblastoma cell lines from classification through regression [22]. A type of L1L2 was also used on eQTL datasets to predict the regulatory potential [23]. The L1L2 regression used in this paper is efficiently implemented through functions from the Open Source mlpy package (also available on MLOSS).
Data Analysis Protocols
Preprocessing
Genotype data are generally encoded with {0, 1, 2} (dominant homozygous, heterozygous, recessive homozygous respectively): this representation has the biological meaning of the number of allelic deviations of the SNP from the dominant homozygous. However, different representations may yield different results. We thus preliminarily compared the proposed encoding with a {−1, 0, 1} encoding, a binary {100, 010, 001}, and one based on the relative frequency of each allelic class for given SNP over the total sample population. For SVR, no significant difference in predictive power was detected between the {0,1,2} and the {−1,0,1} encoding; slightly worse performances were obtained with the frequencybased encoding and significantly worse with the binary {100, 010, 001} encoding. Therefore, we kept to the standard {0,1,2} encoding. Recently, Liu et al. [24] have proposed a sparse binary encoding {00, 01, 11} for a bioinformatics application to ancestry inference. This encoding will be tested in future applications of our regression framework. Finally, missing data (SNPs that are not called) were randomly imputed with probability equal to the relative frequency of each allele at that locus in the population. Random imputation is not biologically plausible, as it will introduce Mendelian errors in the samples and does not take into account linkage disequilibrium effects. More accurate imputation methods for uncalled SNPs have been proposed and are reviewed in [25]. However, the proportion of missing data in the considered dataset is small: 0.14% of all data points are uncalled, with a maximum uncalled samples per SNP of 5.2%. Only 11 SNPs of over 12,000 had more than 2% uncalled samples and only 116 more than 1%. Thus, we expect that the errors introduced by the random imputation will hardly impact the predictive ability of the algorithms.
Model selection
where S_{ p } is the symmetric group on p symbols, and τ_{ i } is the permutation of S_{ p } corresponding to L_{ i } (i = 1,2) for a given order of the features F and . For a given set L of partial lists, the Canberra stability indicator is defined as the mean of all the mutual Canberra distances among the elements of L: the choice of the mean is justified by Hoeffding’s theorem on the asymptotic normality of the distribution of Canberra distances. The Canberra distance was chosen as the (dis)similarity measure because of its intrinsic larger penalization of changes of rank in the top position of the ranked lists. For a complete mathematical description and a few application examples see [26, 27]. The model defined by the optimal (µ, τ, λ) in terms of maximal accuracy and marker list stability was then trained and evaluated on each development/validation split. The L1L2 algorithm and its protocol are implemented within the mlpy Python package and run on Kore, the FBK Linux High Performance Computing facility. It was successfully tested on up to 550k features and a few thousands of samples.
List of abbreviations used
 DAP:

Data Analysis Protocol
 GWAS:

GenomeWide Association Study
 LASSO:

Least Absolute Shrinkage and Selection Operator
 MCH:

Mean Cell Haemoglobin
 MCMC:

Monte Carlo Markov Chain
 RLS:

Regularized Least Squares
 SNP:

Single Nucleotide Polymorphism
 SVM:

Support Vector Machine
 SVR:

Support Vector Regression.
Declarations
Acknowledgements
The authors acknowledge funding by the European Union FP7 Project HiperDART and by the Italian Ministry of Health Project ISITAD (RF 2007 conv. 42).
This article has been published as part of BMC Bioinformatics Volume 11 Supplement 8, 2010: Proceedings of the Neural Information Processing Systems (NIPS) Workshop on Machine Learning in Computational Biology (MLCB). The full contents of the supplement are available online at http://www.biomedcentral.com/14712105/11?issue=S8.
Authors’ Affiliations
References
 Lee SH, van der Werf JHJ, Hayes BJ, Goddard ME, Visscher PM: Predicting unobserved phenotypes for complex traits form wholegenome SNP data. PLoS Genetics 2008, 4(10):e1000231. 10.1371/journal.pgen.1000231PubMed CentralView ArticlePubMedGoogle Scholar
 Casci T: Fitting phenotypes. Nature Reviews Genetics 2008, 9: 896–897. 10.1038/nrg2495View ArticleGoogle Scholar
 Cupples LA, Beyene J, Bickeboller H, Daw EW, Fallin MD, Gauderman WJ, Ghosh S, Goode E, Hauser E, Hinrichs A, Kent J, Martin L, Martinez M, Neuman R, Province M, Szymczak S, Wilcox M, Ziegler A, MacCluer J, Almasy L: Genetic Analysis Workshop 16: Strategies for genomewide association study analyses. BMC Proceedings 2009, 3(Suppl 7):S1. 10.1186/175365613s7s1PubMed CentralView ArticlePubMedGoogle Scholar
 Moore JH, Asselbergs FW, Williams SM: Bioinformatics challenges for genomewide association studies. Bioinformatics 2010, 26(4):445–455. 10.1093/bioinformatics/btp713PubMed CentralView ArticlePubMedGoogle Scholar
 Wooten E, Iyer L, Montefusco M, Hedgepeth A, Payne D, Kapur N, Housman D, Mendelsohn M, Huggins G: Application of Gene Network Analysis Techniques Identifies AXIN1/PDIA2 and Endoglin Haplotypes Associated with Bicuspid Aortic Valve. PLoS ONE 2010, 5: e8830. 10.1371/journal.pone.0008830PubMed CentralView ArticlePubMedGoogle Scholar
 Yosef N, Gramm J, Wang Qf, Noble WS, Karp RM, Roded S: Prediction of Phenotype Information from Genotype Data. Communications in Information and Systems 2010.Google Scholar
 Kraft P, Zeggini E, Ioannidis J: Replication in GenomeWide Association Studies. Statistical Science 2009, 24(4):561–573. 10.1214/09STS290PubMed CentralView ArticlePubMedGoogle Scholar
 Ambroise C, McLachlan G: Selection bias in gene extraction on the basis of microarray geneexpression data. Proc Natl Acad Sci U S A 2002, 99(10):6562–6566. 10.1073/pnas.102102699PubMed CentralView ArticlePubMedGoogle Scholar
 Furlanello C, Serafini M, Merler S, Jurman G: Entropybased gene ranking without selection bias for the predictive classification of microarray data. BMC Bioinformatics 2003, 4: 54. 10.1186/14712105454PubMed CentralView ArticlePubMedGoogle Scholar
 Baggerly K, Coombes K: Deriving Chemosensitivity from Cell Lines: Forensic Bioinformatics and Reproducible Research in HighThroughput Biology. Annals of Applied Statistics 2009, 4(3):1309–1334. 10.1214/09AOAS291View ArticleGoogle Scholar
 Shi L, Perkins RG, Fang H, Tong W: Reproducible and reliable microarray results through quality control: good laboratory proficiency and appropriate data analysis practices are essential. Curr Opin Biotechnol 2008, 19: 10–18. 10.1016/j.copbio.2007.11.003View ArticlePubMedGoogle Scholar
 The MicroArray Quality Control (MAQC) Consortium: The MAQCII Project: A comprehensive study of common practices for the development and validation of microarraybased predictive models. Nature Biotechnology 2010, in press.Google Scholar
 Chierici M, Miclaus K, Vega S, Furlanello C: An interactive effect of batch size and composition contributes to discordant results in GWAS with the CHIAMO genotyping algorithm. Pharmacogenomics Journal 2010, in press.Google Scholar
 De Mol C, Mosci S, Traskine M, Verri A: A regularized method for selecting nested groups of relevant genes from microarray data. Journal of Computational Biology 2009, 16(5):677–690. 10.1089/cmb.2008.0171View ArticlePubMedGoogle Scholar
 Jurman G, Merler S, Barla A, Paoli S, Galea A, Furlanello C: Algebraic stability indicators for ranked lists in molecular profiling. Bioinformatics 2008, 24(2):258–264. 10.1093/bioinformatics/btm550View ArticlePubMedGoogle Scholar
 Valdar W, Solberg LC, Gauguier D, Burnett S, Klenerman P, Cookson WO, Taylor MS, Nicholas J, Rawlins P, Mott R, Flint J: Genomewide genetic association of complex traits in heterogeneous stock mice. Nature Genetics 2006, 38: 879–887. 10.1038/ng1840View ArticlePubMedGoogle Scholar
 Wei Z, Wang K, Qu HQ, Zhang H, Bradfield J, Kim C, Frackleton E, Hou C, Glessner JT, Chiavacci R, Stanley C, Monos D, Grant SFA, Polychronakos C, Hakonarson H: From disease association to risk assessment: an optimistic view from genomewide association studies on type1 diabetes. PLoS Genetics 2009, 5(10):e1000678. 10.1371/journal.pgen.1000678PubMed CentralView ArticlePubMedGoogle Scholar
 Chang CC, Lin CJ: . LIBSVM: a library for support vector machines 2001. Software available at [http://www.csie.ntu.edu.tw/~cjlin/libsvm] Software available atGoogle Scholar
 Smola AJ, Schoelkopf B: A tutorial on Support Vector Regression. Statistics and Computing 2004, 14: 199–222. 10.1023/B:STCO.0000035301.49549.88View ArticleGoogle Scholar
 Zhou H, Hastie T: Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society B 2005, 67(2):301–320. 10.1111/j.14679868.2005.00503.xView ArticleGoogle Scholar
 Tibshirani R: Regression shrinkage and selection via the LASSO. Journal of the Royal Statistical Society B 1996, 58: 267–288.Google Scholar
 Fardin P, Barla A, Mosci S, Rosasco L, Verri A, Varesio L: The l1l2 regularization framework unmasks the hypoxia signature hidden in the transcriptome of a set of heterogeneous neuroblastoma cell lines. BMC Genomics 2009, 10: 474. 10.1186/1471216410474PubMed CentralView ArticlePubMedGoogle Scholar
 Lee SI, Dudley AM, Drubin D, Silver PA, Krogan NJ, Pe’er D, Koller D: Learning a Prior on Regulatory Potential from eQTL Data. PLoS Genetics 2009, 5: e1000358. 10.1371/journal.pgen.1000358PubMed CentralView ArticlePubMedGoogle Scholar
 Liu E, Zhang Q, McMillan L, de Villena FM, Wang W: Efficient genome ancestry inference in complex pedigrees with inbreeding. Bioinformatics 2010, 26(12):i199–207. 10.1093/bioinformatics/btq187PubMed CentralView ArticlePubMedGoogle Scholar
 Marchini J, Cutler D, Patterson N, Stephens M, Eskin E, Halperin E, Lin S, Qin Z, Munro H, Abecasis G, Donnelly P, The International HapMap Consortium: A Comparison of Phasing Algorithms for Trios and Unrelated Individuals. American Journal of Human Genetics 2006, 78: 437–450. 10.1086/500808PubMed CentralView ArticlePubMedGoogle Scholar
 Jurman G, Riccadonna S, Visintainer R, Furlanello C: Canberra Distance on Ranked Lists. In Proceedings, Advances in Ranking – NIPS 09 Workshop Edited by: Agrawal S, Burges C, Crammer K. 2009, 22–27.Google Scholar
 Jurman G, Riccadonna S, Visintainer R, Furlanello C: Algebraic Comparison of Partial Lists in Bioinformatics. 2010. ArXiv:1004.1341v1 [stat.ML]Google Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.