Abstract Practical application of genomic-based risk stratification to clinical diagnosis is appealing yet performance varies widely depending on the disease and genomic risk score (GRS) method. Celiac disease (CD), a common immune-mediated illness, is strongly genetically determined and requires specific HLA haplotypes. HLA testing can exclude diagnosis but has low specificity, providing little information suitable for clinical risk stratification. Using six European cohorts, we provide a proof-of-concept that statistical learning approaches which simultaneously model all SNPs can generate robust and highly accurate predictive models of CD based on genome-wide SNP profiles. The high predictive capacity replicated both in cross-validation within each cohort (AUC of 0.87–0.89) and in independent replication across cohorts (AUC of 0.86–0.9), despite differences in ethnicity. The models explained 30–35% of disease variance and up to ∼43% of heritability. The GRS's utility was assessed in different clinically relevant settings. Comparable to HLA typing, the GRS can be used to identify individuals without CD with ≥99.6% negative predictive value however, unlike HLA typing, fine-scale stratification of individuals into categories of higher-risk for CD can identify those that would benefit from more invasive and costly definitive testing. The GRS is flexible and its performance can be adapted to the clinical situation by adjusting the threshold cut-off. Despite explaining a minority of disease heritability, our findings indicate a genomic risk score provides clinically relevant information to improve upon current diagnostic pathways for CD and support further studies evaluating the clinical utility of this approach in CD and other complex diseases.

Author Summary Celiac disease (CD) is a common immune-mediated illness, affecting approximately 1% of the population in Western countries but the diagnostic process remains sub-optimal. The development of CD is strongly dependent on specific human leukocyte antigen (HLA) genes, and HLA testing to identify CD susceptibility is now commonly undertaken in clinical practice. The clinical utility of HLA typing is to exclude CD when the CD susceptibility HLA types are absent, but notably, most people who possess HLA types imparting susceptibility for CD never develop CD. Therefore, while genetic testing in CD can overcome several limitations of the current diagnostic tools, the utility of HLA typing to identify those individuals at increased-risk of CD is limited. Using large datasets assaying single nucleotide polymorphisms (SNPs), we have developed genomic risk scores (GRS) based on multiple SNPs that can more accurately predict CD risk across several populations in “real world” clinical settings. The GRS can generate predictions that optimize CD risk stratification and diagnosis, potentially reducing the number of unnecessary follow-up investigations. The medical and economic impact of improving CD diagnosis is likely to be significant, and our findings support further studies into the role of personalized GRS's for other strongly heritable human diseases.

Citation: Abraham G, Tye-Din JA, Bhalala OG, Kowalczyk A, Zobel J, Inouye M (2014) Accurate and Robust Genomic Prediction of Celiac Disease Using Statistical Learning. PLoS Genet 10(2): e1004137. https://doi.org/10.1371/journal.pgen.1004137 Editor: Joshua M. Akey, University of Washington, United States of America Received: August 1, 2013; Accepted: December 8, 2013; Published: February 13, 2014 Copyright: © 2014 Abraham et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Funding: The celiac disease genotype data was generated under WTCCC award WT084743 and Coeliac UK funding. Part of this work utilized the computing resources of the Victorian Life Sciences Computation Initiative (project VR0126). This study makes use of data generated by the Wellcome Trust Case-Control Consortium. A full list of the investigators who contributed to the generation of the data is available from www.wtccc.org.uk. Funding for the project was provided by the Wellcome Trust under award 076113 and 085475. MI was supported by an NHMRC early career fellowship 637400. MI and GA were supported by University of Melbourne funding. JATD was supported by an NHMRC Postgraduate Medical Scholarship. This work was partially supported by the Australian Research Council, and by the NICTA Victorian Research Laboratory. NICTA is funded by the Australian Government as represented by the Department of Broadband, Communications, and the Digital Economy, and the Australian Research Council through the ICT Centre of Excellence program. This work was made possible through Victorian State Government Operational Infrastructure Support and Australian Government NHMRC IRIIS. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

Introduction Improving the diagnosis of celiac disease (CD), a common immune-mediated illness caused by dietary gluten, remains a clinical challenge [1], [2]. Despite a prevalence of approximately 1% in most Western countries, lack of awareness and failure to implement appropriate serological, histological and genetic testing means that less than 30–40% of those affected by CD are diagnosed [1], [3]–[5]. Undiagnosed CD is associated with reduced quality of life, substantial morbidity, and increased mortality, however, prompt diagnosis and treatment lowers the burden of disease and may reduce the rate of complications such as osteoporosis, autoimmune disease, and malignancy. Optimizing the diagnosis of CD is now recognized as an important goal for clinicians [6]. CD is characterized by a variable combination of gluten-dependent clinical manifestations, CD-specific antibodies and small bowel inflammation (villous atrophy) [7]. Traditional guidelines for the diagnosis of CD rely on demonstrating villous atrophy and improvement of symptoms, laboratory abnormalities, and/or small bowel inflammation upon exclusion of dietary gluten [8]. Current clinical practice is to screen for CD by detecting CD-specific serum antibodies and then confirm the diagnosis by undertaking small bowel biopsy to demonstrate typical villous atrophy. Serologic screening for CD with transglutaminase-IgA antibodies is reported to be highly sensitive and specific for CD (both >90%), imparting a high positive predictive value (PPV) of over 90% when assessing most populations [9], [10], although the PPV can fall to 45–70% in community screening settings [11], [12]. In practice, serological and histological assessments have technical limitations that generate both false negative and false positive diagnoses. A key feature of CD is its strong dependence on the presence of susceptibility genes encoding for HLA DQ2.5, DQ8, and/or half the HLA DQ2.5 heterodimer (typically DQ2.2), seen in approximately 99.6% of all patients with CD [13]. These genes encode immune-recognition molecules which facilitate CD4+ T cell recognition of specific gluten-derived peptides, a critical step in disease pathogenesis [14]–[18]. Recognizing the crucial role of these genes, the latest consensus diagnostic guidelines for CD recommend testing for these HLA heterodimers (HLA typing) as a first-line investigation for asymptomatic individuals identified at-risk of CD, such as 1st-degree relatives of an affected individual or those with suggestive symptoms [7]. However, a major flaw of HLA typing as a diagnostic tool is that a substantial proportion of the community, typically reported to be 30–40%, express HLA DQ2.5, DQ8, and/or DQ2.2, thus making the presence of these HLA types poorly predictive and of low specificity for CD [13]. Indeed, a recent Australian population study revealed that 56% of the community possessed at least one of these CD susceptibility haplotypes [5]. Thus, while HLA typing can exclude CD in the community with high confidence when the susceptibility haplotypes are absent, these haplotypes will be present in 30–56% of the population, the majority of whom would not have CD. Therefore, if assessed as a stand-alone test, HLA typing has exceptionally high sensitivity and negative predictive value (NPV) but very poor specificity and low positive predictive value (PPV) for CD. Since a positive result poorly predicts the presence of CD, HLA typing is not useful as a stand-alone diagnostic tool for CD. While the relative-risk for CD can be stratified based on the HLA subtype (CD risk DQ2.5>DQ8>DQ2.2) [19], these categories have low positive predictive value and do not provide clinically-informative attribution of CD risk [20]; HLA results are therefore interpreted as a binary outcome: CD susceptibility positive or negative. Despite these limitations, HLA typing is now widely utilized in clinical practice and typically determined using polymerase chain-reaction sequence specific oligonucleotide (PCR-SSO) hybridization, which is time and labor intensive, and costly (AU $120/sample, Medicare; in the USA cost varies but is typically US $150/sample or greater). It is important to distinguish between three different approaches to analyzing the HLA region for association with CD. The first approach, currently in clinical practice, is HLA typing, as described above, where the HLA result is considered a binary variable and its utility is to exclude CD. A second approach, such as that taken by Romanos et al., utilizes the same HLA-DQ haplotypes, stratifies individuals into several nominal risk levels then fits a statistical model to empirically estimate the true risk in each group [21], [22]. While HLA-DQ haplotypes may be inferred from typing several HLA SNPs, importantly the HLA SNPs are only used to assign the HLA type and the SNPs themselves are not directly modeled. The third approach, such as that used here, is based on direct concurrent modeling of many thousands of individual SNPs for association with CD in order to produce a more fine-grained predictive “genomic risk score” (GRS). GRSs have been enabled by the advent of genome-wide association studies (GWAS), which perform unbiased testing of many thousands of SNPs for association with CD. Using GWAS, recent studies have identified multiple non-HLA SNP associations with CD [23], [24]. GWAS are primarily concerned with the detection of variants associated with disease in order to gain insight into the disease etiology and genetic architecture. Due to the high number of significance tests, controlling for false positive associations is a major, valid concern. Therefore, SNP-based risk scores have tended to be constructed from the SNPs found to be significantly associated with the disease status [22], [25]. However, due to the stringent multiple-testing corrections utilized in GWAS there may be other SNPs that fail to achieve genome-wide significance but may be predictive of disease status nonetheless and including them in the model could potentially result in higher predictive ability than achievable by models based solely on genome-wide significant SNPs. In contrast to the GWAS approach, the main overriding aim of a GRS from a clinical perspective is to achieve maximal predictive capacity, the inference of genetic architecture is secondary. We have recently designed computational algorithms which efficiently fit L1-penalized multivariable classification models to genome-wide and whole-genome SNP data [26]. Such models were then shown to be preferable to several other methods such as the standard method of summing the per-SNP log odds (polygenic score) [27], mixed effects linear modeling [28], [29], and unpenalized logistic regression, with both better precision for detecting causal SNPs in simulation and better case/control predictive power [30]. These advantages were consistent across several complex diseases, including two British studies of CD. However, the diagnostic implications of penalized models have not been previously examined nor has the robustness of such models in other populations or the advantage over HLA-typing approaches. In contrast to existing studies that examine a small number of genome-wide significant SNPs, we have shown that many more SNPs (potentially hundreds) are required to achieve optimal predictive ability for CD. Further, the standard GWAS approach of considering each SNP separately when estimating its effect size does not consider its correlation with other SNPs. We have shown that unpenalized predictive models based on these top SNPs suffer from lower predictive ability than L1-penalized models since the pre-screening introduces multiple highly correlated SNPs into the model, of which a substantial proportion may be redundant in terms of contribution to the predictive ability. Similar L1-penalized approaches have also recently been successfully applied to inflammatory bowel disease case/control Immunochip data, where models based on several hundred SNPs have led to high predictive ability [31]. Here, we provide a proof-of-concept that the GRS for CD, induced by L1-penalized support vector machine models, are able to achieve a predictive capacity and robustness that provides information not afforded by current diagnostic pathways utilizing HLA typing alone. This GRS has the potential to provide greater clinical diagnostic utility by enabling each individual to be assigned a more informative risk score beyond the simple designation of “CD susceptible” or “CD non-susceptible”, or “high risk” versus “low risk”. To enable useful comparisons between diagnostic approaches, we model the GRS as a stand-alone test to “diagnose” CD, while at the same time acknowledging that real world clinical practice will need to draw upon clinical history, CD-specific serology and small bowel histology to confirm the diagnosis of CD. We assess the predictive power of the GRS both in cross-validation and in external validation, across six different European cohorts, showing that the models strongly replicate. We test our GRS on three other autoimmune diseases: type 1 diabetes, Crohn's disease, and rheumatoid arthritis, finding some predictive ability for T1D status but none for the others, thus largely supporting the specificity of the scores for CD. To overcome limitations of previous studies utilizing GWAS case/control studies, where ascertainment bias incurs substantially higher rates of false positive results, we undertake genomic prediction of CD in “real world” settings where the prevalence of CD is far lower and evaluate the performance of the GRS using PPV and NPV at several levels of CD prevalence. Unlike HLA typing, the GRS allows flexibility in determining who is considered at higher risk for CD by selecting a clinically determined user-specified threshold. We demonstrate how these scores can be practically applied at various prevalence levels to optimize sensitivity and precision. Finally, we show how the model can be calibrated to produce accurate predicted probabilities of disease.

Discussion In this study, we have sought to exploit the strong genetic basis for CD and leverage comprehensive genome-wide SNP profiles using statistical learning to improve risk stratification and the diagnosis of CD. Our models showed excellent performance in cross-validation and were highly replicable in external validation across datasets of different ethnicities, suggesting that the genetic component is shared between these European ethnicities and that our models were able to capture a substantial proportion of it. Importantly, even without explaining a majority of CD heritability, the models were robust and accurate, showing that it is not necessary to explain most of the heritability in order to produce a useful model. The most frequently employed tools to diagnose CD are serology and small bowel histology, but both have limitations. Differences in the sensitivity of antibody recognition of commercially employed CD-specific antigens such as tissue transglutaminase, deamidated gliadin peptides, and endomysial antigen, as well as the human operator performing the assay can all influence findings and affect reproducibility of serological testing [9], [46]–[49]. Serologic testing in children is reported to be less reliable before the age of 4 and up to 50% of children normalize elevated antibodies over time [50], [51]. While small bowel histology remains the ‘gold standard’ confirmatory test, it is dependent upon patients willing and available to undergo endoscopy, adequate sampling by the gastroenterologist, and appropriate pathological processing and interpretation [52]–[54]. The frequencies of false positives and false negatives in CD serology assays vary widely and also partly depend upon what degree of histologic inflammation is considered compatible with CD [52], [54]–[58]. Notably, the accuracy of both serologic and histologic testing for CD is dependent on the ongoing consumption of gluten. It is clear that clinically significant variability exists in serologic and histologic work-up for CD and new tools to improve the accuracy of CD diagnosis would be of benefit to clinicians. Given the strong genetic basis for CD, genomic tools are logical and appealing because they are relatively robust and less subject to the kind of variability seen with serologic and histologic assessment, are independent of age, and do not rely on dietary intake of gluten. A major shortcoming of clinical HLA typing for risk prediction of CD is its poor specificity. HLA testing would result in virtually all CD cases detected but at the cost of approximately 30–56 people incorrectly implicated for each true case of CD. A significant advantage of the GRS approach is that it can be adapted to the clinical scenario in order to maximize PPV and diagnostic accuracy. By promoting accurate clinical stratification, the GRS could reserve invasive and more expensive confirmatory testing for those who would most likely benefit from further investigation to secure a diagnosis, and it would avoid unnecessary procedures in those who are HLA susceptible but unlikely to have CD. This provides both clinical and economic benefits. HLA typing does not provide the flexibility afforded by the GRS and cannot be effectively employed to identify those who would benefit from endoscopy. For instance, if HLA typing were used as a guide for further investigations, at 10% CD prevalence it would generate over five unnecessary endoscopies per correct endoscopy and at 1% CD prevalence it would generate 30–56 unnecessary endoscopies. Small bowel endoscopy is not a trivial undertaking – the procedure is costly (approximately AUD $750–$1000 for the procedure and associated pathology), has potential complications, necessitates a full day off work, and many patients are reluctant to undergo it. The GRS can be used to exclude patients unlikely to have CD with a performance comparable to HLA typing. Testing with these parameters may be useful in the clinical scenario of assessing individuals at average risk of CD. A common example would be when a person has commenced a gluten-free diet prior to assessment for CD by serology or small bowel examination and are unwilling or unable to resume oral gluten intake in order to make testing reliable. This is an increasingly common clinical dilemma as the number of people following a gluten-free diet without adequate initial testing for CD continues to rise. In the United States approximately 30% of the adult population are interested in cutting back or avoiding dietary gluten [59]. The GRS can also be used to stratify the risk for CD in patients who present with suggestive clinical features. These risk factors include having a first-degree relative with CD or problems such as recurrent abdominal pain, bloating, diarrhea or constipation, fatigue, weight loss, unexplained anemia, autoimmune disease (including thyroid disease, T1D, autoimmune hepatitis, rheumatoid arthritis, and Sjogren's syndrome), infertility or early-onset osteoporosis [3], [60]. Supporting the recently revised diagnostic guidelines for CD, which promote HLA testing as the 1st line investigation for higher-risk cases, genetic testing of CD is likely to be more informative in these sub-populations exhibiting higher-than-normal prevalence. While clinical guidelines recommend screening for CD in these high-risk populations [61], testing often poses a diagnostic dilemma as serologic assessment alone cannot confidently exclude a diagnosis, especially given the higher pre-test probability. HLA typing is not particularly informative as the CD HLA susceptibility haplotypes HLA-DQ2.5 and DQ8 are commonly present (manifesting in over 90% of patients with T1D and in 65% in first-degree relatives of individuals with CD) [62], [63]. Stratifying these higher-risk patients based on a GRS will allow improved identification of those where small bowel biopsy is likely to be informative. Thus, a GRS should reduce the number of unnecessary small bowel biopsies in first-degree relatives who carry HLA susceptibility for CD but do not have it. We have found that our CD models had only moderate predictive ability for T1D, which is consistent with previous findings showing some shared genetics between T1D and CD [38]. Despite the substantial overlap of genetic factors for autoimmune disease, the CD models had negligible predictive ability for Crohn's disease and rheumatoid arthritis. These results indicate that our GRS is specific to CD and less likely to incorrectly identify patients with other autoimmune diseases as having CD, but further work is required to determine whether CD can be as confidently predicted in individuals with T1D as it is in non-T1D populations. Another major clinical challenge that may benefit from genomic risk prediction is determining the natural history of potential CD (formerly termed ‘latent CD’) when there is serologic but not histologic evidence of CD, and identifying which patients are more likely to develop overt CD with small bowel inflammation [64]. Current practice is to follow-up all patients with immunologic evidence of gluten intolerance in order to capture those who will eventually develop overt disease. An analogous clinical scenario is that of children with positive CD serology, of whom 50% will fail to develop small bowel changes consistent with CD during follow-up [50], [51]. In both clinical situations, it is reasonable to expect that a GRS can improve risk stratification of such patients for developing overt CD. Of course, environmental factors are important in the development of CD and the exact extent to which environmental versus genetic factors contribute to the development of overt CD remains unknown. Long-term follow-up studies of patients with potential CD will be necessary to establish the role of genomic risk prediction in this important subgroup. Future work will look at optimizing our GRS as a tool to predict CD risk. Validation of our model in real-life practice will be important to confirm the clinical benefit of the GRS in conjunction with serology and/or over HLA typing alone, as well as to what extent other clinical predictors such as sex, age, and family history can contribute to clinically relevant risk prediction. Future prospective studies will enable direct optimization of clinical utility (accuracy, practicality, throughput and cost) afforded by the GRS, for example in conjunction with CD serology. These studies will also provide a rigorous evidence base for suggested clinical guidelines of GRS usage. Importantly, appropriate GRS cut-off levels to maximize diagnostic accuracy (optimal PPV and NPV for each given clinical scenario and CD prevalence) could be obtained by local prospective validation. Such studies can identify the ultimate clinical role for the GRS: whether it can effectively replace HLA typing and also whether it is a stand-alone test or one to accompany CD serology. Hadithi et al showed that in patients at high-risk of CD the addition of HLA typing to CD serology had the same performance as either testing strategy alone [65], but the greater precision of the GRS over HLA typing may better complement CD serology. Understanding where the GRS fits in the diagnostic algorithm to optimize precision and cost-effectiveness will be essential, as is the role it might play in the diagnostic work-up of CD in populations with lower levels of clinical risk. Health economic modeling will address the cost-benefits of using the GRS in the diagnosis of CD, taking into account the cheaper cost of GRS over HLA typing, and include the downstream benefits of potentially reducing endoscopies (substantial cost savings and value to patients from reduced discomfort) as well as potential improvements in quality of life from the detection of CD. Further, it may be that other statistical modeling approaches yield improvements in predictive power, for example non-additive models that consider epistatic interactions between SNPs. Another avenue for improvement is considering each CD subtype separately, recognizing potentially different genetic bases for these conditions. Based on our results, we do not expect substantial improvements from increasing sample size alone, however this will be important for adequately powered studies of lower frequency genetic variants of assumedly greater effect size. In summary, this study demonstrates that simultaneous modeling of all SNPs using statistical learning was able to generate genomic risk scores that accurately predict CD to a clinically relevant degree. This was despite the models explaining only a minority of disease heritability. The GRS better enables clinicians to stratify patients according to their risk of CD compared to HLA typing alone and, we predict, more accurately determines those suitable for confirmatory testing in the form of small bowel biopsy. Reserving this invasive, time consuming and costly procedure for higher-risk cases is likely to improve the accuracy, cost and public acceptance of testing for CD, and by extension, benefit the overall diagnosis of CD in the community. By better prioritizing higher-risk patients for confirmatory testing, genomic risk prediction carries promise as a clinically useful tool to add to the clinician's diagnostic armamentarium. Ultimately, we envisage a clinical scoring algorithm based on the combination of clinical features, serologic, and genetic information that will accurately predict people with biopsy-confirmed CD and perhaps ultimately overcome the reliance on small bowel histology altogether. Further, the costs of genotyping a select number of marker SNPs with a low-plex, high throughput technology are already far lower than the costs of full HLA typing, resulting in a test that is cheaper, more flexible and more precise than HLA typing. More generally, this study demonstrates that statistical learning approaches utilizing SNPs can already produce useful predictive models of a complex human disease using existing genotyping platforms assaying common SNPs and suggests that similar approaches may yield comparable results in other complex human diseases with strong genetic components.

Methods Ethics Statement All participants gave informed consent and the study protocols were approved by the relevant institutional or national ethics committees. Details given in references van Heel et al [23] and DuBois et al [24]. All data was analysed anonymously. Data We analyzed six CD datasets: UK1 [23], UK2, IT, NL, and Finn [24], and IMM [66]. The main characteristics of the datasets are listed in Table 1. In addition we used three WTCCC datasets (T1D, Crohn's, and RA) that have been described elsewhere [30], [37]. UK1 used the Illumina Hap330v1-1 array for cases and Hap550-2v3 for controls, UK2 used the Illumina 670-QuadCustom-v1 for cases and 1.2M-DuoCustom-v1 for controls, the NL and IT datasets used the Illumina 670-QuadCustom-v1 in both cases and controls, and the Finn dataset used the Illumina 670-QuadCustom-v1 for cases and Illumina 610-Quad for controls. The WTCCC data (T1D, Crohn's, and RA) used the Affymetrix 500K array. In all of our models, we used autosomal SNPs only, and did not include the gender as a covariable, as models built separately on the two genders using the same sample size and case∶control balance showed very similar performance in cross-validation on the UK2 dataset (results not shown). For analyses of the MHC region, we defined the MHC as all SNPs on chr6 in the range 29.7 Mb–33.3 Mb. Quality Control For each of the UK1, UK2, IT, NL, and Finn datasets, we removed non-autosomal SNPs, SNPs with MAF<1%, with missingness >1%, and those with deviations from Hardy-Weinberg Equilibrium in controls P<5×10−6. We also removed samples with missingness >1%. We tested identity-by-descent between samples in UK1 and UK2 and removed one of a pair of samples with pi-hat ≥0.05 (either between the datasets or within the datasets). The QC for the IMM Immunochip data has been previously described [66]; we estimated 5763 Immunochip samples to have pi-hat ≥0.125 (PLINK IBS) with any UK2 sample, and those were removed, leaving 10,304 Immunochip samples in total, with 18,252 SNPs shared with the UK2 dataset (post-QC). The QC for the WTCCC data (T1D, Crohn's, and RA) has been previously described [30], [37]. Assessment of Population Structure Effects To assess the impact of potential cryptic population structure, we estimated the top 10 principal components (PCs) for the UK2 with EIGENSOFT 4.2 [67], after removal of regions with high LD (see Text S1 for details). The principal components themselves showed almost no predictive ability (AUC = 0.52), and models trained on all SNPs accounting for these PCs showed indistinguishable performance from the non-adjusted model, both in cross-validation on the UK2 dataset and in external validation on the Finn, NL, and IT datasets (Figure S5), demonstrating that confounding of our UK2 models by population structure was negligible and was not a contributing factor to the high predictive ability. Statistical Analysis We used L1-penalized support vector machines (SVM) implemented in the tool SparSNP [26] (https://github.com/gabraham/SparSNP) as the classifiers. The L1-penalized SVM is a sparse linear model, that is, many or most of the SNPs will receive zero weight in the model, as determined by the L1 penalty. The use of a sparse model fits with our prior expectation that in autoimmune disease most SNPs will not be associated with disease status. The inherent sparsity of the model obviates the need for subsequent filtering of SNPs by weight, in order to decide which ones show strong evidence of association and which are spurious, as would be required in a non-sparse (L2-penalized) model. In addition, in extensive simulation and in analysis of real genotype data, including the two celiac disease datasets UK1 and UK2, we have previously shown the advantage of L1-penalized SVMs over commonly used approaches such as polygenic scores (sum of the log odds), linear mixed models (GCTA), and unpenalized logistic regression [30]. The advantage of sparse models over standard linear mixed models in predicting autoimmune disease has been recently confirmed in type-1 diabetes as well [68]. We have also shown that our L1-penalized SVMs achieved essentially identical performance to L1-penalized logistic regression (glmnet) in cross-validation over the Finnish subset of the celiac disease dataset, while being substantially faster [26]. Unlike single marker approaches that estimate the effect size of each SNP separately, the L1-penalized SVM is a multivariable model, where the estimated effect of each SNP is conditional on all other SNPs, thereby implicitly accounting for the linkage disequilibrium (LD) between SNPs. Besides imposing sparsity, the L1 penalty tends to produce models where one representative SNP is selected out of a group of highly correlated SNPs, while the rest remain with a zero weight, in contrast with L2-penalized or unpenalized models where many or all of these SNPs may receive a non-zero weight. For an in-depth discussion of these issues and the effects of varying LD levels on the performance of multivariable models, see [30]. The L1-penalized SVM model is induced by minimizing the L1-penalized squared-hinge loss over N samples and p SNPs, where x i is the p-vector of genotypes for the ith sample in allele-dosage coding {0, 1, 2}, y are the binary phenotypes {−1, +1}, β is the p-vector of weights, β 0 is the intercept (also called the bias, which is not penalized), and λ is the L1 penalty. We also investigated adding an L2 penalty to the model (elastic-net), however, based on initial cross-validation experiments, we found no advantage in the L2 penalty and subsequently did not use it. All of our models were additive in the allele dosage {0, 1, 2}. The genomic risk score for a new sample x i consisting of p genotypes is then where the continuous value is later thresholded at different values to produce a binary predicted class. The model was evaluated over a grid of penalties, in 10-fold cross-validation, repeated 10 times. The optimal number of SNPs in the model was decided based on the model with the highest average AUC across the replications. The final model was a consensus model, averaged over all 10×10 = 100 models, and containing approximately the number of SNPs determined earlier. Post processing and plotting of the results was performed in R [69], together with the package ggplot2 [70]. Measures of Predictive Performance To quantify the predictive performance of the models in cross-validation and external validation, we employed receiver operating characteristic (ROC) curves (sensitivity versus 1 minus specificity), the area under the ROC curve (AUC) [71], and the proportion of phenotypic variance explained [32]. To quantify predictive performance in different population settings, we used the positive and negative predictive values, which can be estimated as and where “sens” is the sensitivity = TP/(TP+FN), “spec” is the specificity = TN/(FP+TN), and “prev” is the population prevalence. The PPV/NPV are equivalent to the posterior probability of a person having/not having the disease given a positive/negative diagnosis, respectively. When the PPV and precision are estimated in data with identical prevalence (that is, the observed prevalence in the data is identical to the prevalence in the population for which we wish to estimate PPV), they are equivalent. Precision is defined as TP/(TP+FP).

Supporting Information Figure S1. LOESS-smoothed AUC in 10×10-fold cross-validation for the random subsamples of the UK2 dataset, in increasing sample size proportions of the original data (n = 6785). https://doi.org/10.1371/journal.pgen.1004137.s001 (EPS) Figure S2. Results of externally validating the predictive models, trained on UK2 in cross-validation, and tested on the other CD datasets. Legend: Romanos HLA: 3-levels of risk (low, medium, high) [21] based on imputed HLA type (HIBAG); Romanos HLA+57 SNPs (Immunochip only): 3-level HLA risk plus 57 Immunochip non-HLA SNPs [21]; Monsuur HLA SNPs: logistic regression on individual HLA SNPs [36] (5/6 SNPs or proxies thereof were found in the UK2/Finn/NL/IT datasets, 3/6 were found in the subset of UK1 shared with UK2); GRS MHC SNPs: SparSNP run on individual SNPs on chr6 within 29.7 Mb–33.3 Mb; GRS non-MHC SNPs: SparSNP run on individual autosomal SNPs outside MHC; GRS all SNPs: SparSNP run on all autosomal SNPs. https://doi.org/10.1371/journal.pgen.1004137.s002 (EPS) Figure S3. ROC curves for CD model trained on SNP subsets of the UK2 dataset that were assayed for the WTCCC-T1D dataset: All SNPs (76,847 SNPS), MHC SNPs (186 SNPs in the MHC region of chr6, 29.7–33.3 Mb), and Non-MHC SNPs (76,661 SNPs outside the MHC). https://doi.org/10.1371/journal.pgen.1004137.s003 (EPS) Figure S4. Calibration plots, comparing predicted score in 5% quantiles against observed proportions of cases falling within the bin. The score comes from models trained on the UK2 dataset, and tested on the rest of the datasets. The bars show 95% confidence intervals using the Agresti-Coull method for proportions. We randomly split the test datasets into two halves. In the first half, we plotted the original quantiles of the scores and fitted a LOESS smooth to them. We did this for the original case/control data (prevalence of 40%), shown in (a), and for a subsampled version of the data with prevalence of 10% (c). We then used the LOESS smooth to correct the original quantiles, forming a calibrated score, one for each dataset (Finn, IT, NL, UK1), which was then applied to the second half of the data, shown in (b) and (d) for prevalence of 40% and 10% respectively. The second half of the data was not used in the calibration step. https://doi.org/10.1371/journal.pgen.1004137.s004 (EPS) Figure S5. (a) LOESS-smoothed AUC from 10-fold cross-validation for the UK2 model (all autosomal SNPs), accounting for the top 10 PCs (included in training but not in testing). (b) External validation of the best UK2 model that accounted for the PCs (PCs excluded from testing). https://doi.org/10.1371/journal.pgen.1004137.s005 (EPS) Table S1. The predictive model. The SNPs are sorted in decreasing order of the absolute value of their model weight averaged over the 10×10 cross-validation folds. Stability is the percentage of times a SNP was selected to have non-zero weight over the 10×10-cross-validation folds. Intercept: −0.757226. To annotate the SNPs we used Bioconductor 2.12 together with the packages VariantAnnotation 1.6.5 and TxDb.Hsapiens.UCSC.hg18.knownGene 2.9.0. We considered a SNP to be genic if it was annotated to fall inside one of the regions {spliceSite, intron, fiveUTR, threeUTR, coding, promoter} and intergenic otherwise. For intergenic SNPs, we also annotate the nearest gene and the distance to it. All positions are in hg18 coordinates. https://doi.org/10.1371/journal.pgen.1004137.s006 (PDF) Table S2. Summary of screening results at different prevalence levels for the combine dataset Finn+IT+NL dataset, using different cutoffs to declare the samples as disease cases (expressed as % of the population). The smaller the cutoff, the stricter the definition of a disease case. https://doi.org/10.1371/journal.pgen.1004137.s007 (PDF) Text S1. Supplementary Methods. https://doi.org/10.1371/journal.pgen.1004137.s008 (PDF)

Acknowledgments We thank the investigators of the van Heel et al., 2007, Dubois et al., 2010 and Trynka et al., 2011 papers (David van Heel, Cisca Wijmenga, and Lude Franke) for providing the celiac disease data.

Author Contributions Conceived and designed the experiments: MI GA. Performed the experiments: GA. Analyzed the data: GA MI OGB JZ AK. Contributed reagents/materials/analysis tools: GA JZ AK JATD. Wrote the paper: GA MI JATD.