dataesr/scientific-paragraphs-categorization huggingface dataset:
Detecting genetic interaction would also greatly improve the explained genetic variance. However, optimal methods for selecting and combining SNPs remain to be developed. Several methods have been proposed to analyze whole genotyping datasets looking for genetic interactions (for a review and their application to CD see  8, 9  ). A search for specific pairwise gene-gene interactions of known genetic factors was initially performed using statistical methods but it failed to identify genetic variants with strong interactions  10, 11. Genome-wide scans looking for two-loci interactions also failed to identify statistically significant epistatic effects  8, 12. Standard methods applied to higher orders of interactions were quickly limited by the multi-testing issue, the size of the datasets and computing power.
For these reasons, more sophisticated machine learning (ML) methods have been proposed in order to capture the whole information of GWAS datasets using a direct pan-genomic approach. To explore the performances of different methods, the receiver operator characteristic (ROC) curve and its maximum Area Under Curve (AUC) are often used to compare the sensitivity and specificity of genetic tests in correctly classifying affected and unaffected individuals  2. The prediction of a specific phenotype such as "CD" or "non-CD" from raw genomic data such as SNPs can be thought in the framework of supervised learning as a binary classification problem. During training, the algorithm learns from a genotyping dataset and adjusts its internal parameters to minimize an error cost function between predicted probabilities of phenotypes and actual patients phenotypes.
After training, the performance of the system is measured on another set of a completely new case/control sample, to evaluate the generalization ability of the proposed algorithm. Linear ML methods including multivariate logistic regression (LR) and sparse penalized methods such as Lasso have been proposed to identify disease associated SNPs  13. Penalized methods perform two main tasks in the same process. First, it identifies a set of SNPs involved in disease prediction. Second, it calculates a weight for each SNP which reflects its contribution to the general model. Using the Welcome Trust GWAS dataset, Abraham et al. showed that penalized models achieved better performance than non-penalized methods even if the maximum AUC was no more than 0. 76  13. Using either GWAS datasets and/or Immunochip datasets, Chen et al.
found an AUC of 0. 80  14. Wei et al. applied a LR with L1 penalty to the Immunochip dataset and obtained an AUC of 0. 86  15. The observed differences are likely explained by the different datasets used and the quality controls (QC) applied  13, 14. Penalized or non penalized LR are not suited to capture non linear interactions between loci because they only capture additive risk contributions. A non-linear model achieving better disease prediction results (in term of AUC score in our case), would allow to evaluate the amount of extra-information related to such interactions and provide some hint about their nature. To explore non linear interactions, ensemble tree-based methods such as random forest and gradient boosted trees (GBT) have been proposed  8, 16, 17. In these methods, which are often providing state-of-the art results for structured datasets, ensembles of decision trees are trained to minimize the prediction error.
Each neuron treats the information from several inputs (like the presence or not of a specific allele) and produces an output which is itself used by another neuron in the following layer. This structure draws multiple levels of representation. Starting from the raw inputs, it forms non-linear combinations of the initial features, each one providing a new representation of the initial data. By adding deeper layers, the networks build higher levels of abstraction and complexity which could be interpreted as biological functions. Deep network are known to amplify relevant inputs and lower the noise in data like images and sounds. However, the application of these approaches to standard tabular datasets does not generally outperform ensemble tree-based methods and it remains an open challenge to design efficient deep learning systems for this type of data.
Despite its ability to explore intricate structures in high-dimensional data  22, 23  , to our knowledge it has rarely been applied to population genetics and GWAS datasets  24. IIBDGC has collected a large dataset from CD patients and healthy controls genotyped for more than 150 thousands genetic variants -mainly single nucleotide polymorphisms (SNPs)-forming the Immunochip panel  3, 25. Using this dataset, we first explored the impact of QC methods, allele coding strategies and marker selections on the test accuracy of Lasso as the reference method. Second, we applied a panel of different penalized logistic regression, GBT and NN methods. Finally, we explored the robustness of the three approaches and their complementarity.
S1 ), genetic variants previously associated to immune disorders being over-represented  25. In the following we call A the major allele and a the minor one at a given site. Due to the biallelic nature of SNPs under study, in the dataset we substituted them by numerical values 0 and 1, and call N SNP the number of SNPs. Quality control and imputation. It is well known that QC and imputation on GWAS data are delicate pre-processing steps for any genotype-phenotype association analysis, and that they can strongly affect results and biological interpretation  [26] [27] [28]. Since one of the aim of this paper is to compare different ML approaches to the classification case-control problem, we first addressed the question of the impact of such an issue, on the AUC scores and on the interpretability of the feature importance selection.
In particular, for this part of the analysis, we have considered two cases: NoQC -All samples and SNPs are kept for the analysis. QC -The IIBDGC dataset is pre-processed by applying the cuts on samples and SNPs missing rates as in  3. In particular, after excluding samples with missing SNP rate greater than 5%, SNPs with sample missing rate greater than 2% and with Hardy-Weinberg equilibrium (HWE) p-value < 10 -10 in controls, we are left with 17966 CD cases, 33985 controls, for 146237 SNPs. Moreover we proposed three different schemes to treat unknown genotypes in the dataset (notice that more sophisticated algorithms devoted to imputation, like for example IMPUTE2  29  , could in principle perform better. However, in this work we focused on comparing different ML algorithms on a given unbiased dataset.
As we show in the Results Section, strategy B3 is sufficient to this aim). For any given allele in the SNP we: Unkw -leave the unknown values and treat them as a separate allele, or Maj -substitute them with the most common allele, or HW -substitute them with a random choice, following a binomial distribution that satisfy the Hardy-Weinberg equilibrium, for the controls (HW c ), or for all samples (HW a ). Coding. The dataset is unphased, namely alleles in each SNP cannot be ordered for a given sample (phase information does not seem to strongly improve the results in similar setups  30  ). Therefore only 3 classes can be associated to each SNP. Different models can be considered, depending on the assumptions on these 3 classes. In the additive model, the SNP genotypes are ordered numbers: without loss of generality one can fix AA = 0, Aa = 1, aa = 2, implying that each additional number of copies of the minor allele increases the risk by the same amount.
This coding method is referred as sum. A dominant model compares AA versus Aa + aa, and a recessive model compares AA + Aa versus aa, giving rise to only two effective classes (0 and 1). Finally, the three classes can be considered independently, if no strong assumptions can be made about dominance or additivity. This can be achieved in two different ways. The first option is to use a One-Hot Encoding on the three classes, giving rise to an effective number of 3 * N SNP features in the dataset. We refer to this coding as the OHE coding. The second possibility is to keep the data as in its raw version (then called raw coding in the following), meaning considering each allele as an independent feature. It is easy to convince themselves that this is equivalent to consider the 3 independent classes (0-0, 0-1, 1-1).
The advantage of this second option with respect to the OHE model is that the number of effective features is now 2 * N SNP. In this paper, we focus on sum, OHE and raw coding. When interested to a not imputed dataset (referred above as Unkw), only categorical features make sense and OHE coding. Since in the Immunochip dataset the alleles for a given SNP are always either both known or both unknown, in this case one can code only the 4 classes 0-0, 0-1, 1-1, U-U (where U is an unknown allele). Data separation, score, cross-validation. Samples have been randomly permuted, in order to obtain a similar case/control ratio on all subsets. Then the dataset has been separated in a Train dataset (34634 samples) and a Test dataset (17317 samples). In both Train and Test sets, the case/control ratio was around 0.
Nevertheless, different contingency matrices and corresponding chi-squared test can be studied, depending on the genetic model taken into account. ML algorithms can in principle deal with genome-wide SNPs. However, dataset with a large number of features are subject to the curse of dimensionality. Therefore, a more efficient strategy consists in first reducing the total number of SNPs to a manageable level via a screening procedure, and look for causal loci among those passing a given threshold  15, 35, 36. In  15  for example, a SNP preselection based on minimal allele frequency (MAF) threshold and p-values on single SNP association tests, is followed by a LR with Lasso regularization. Since in this case the aim is just to reduce the dimensionality of the space in which the data are embedded, less stringent thresholds than in association studies are used.
Models and implementation. In this paper we considered three classes of models for case/control classification: logistic regression (LR), dense neural networks (NN) and gradient boosting on decision trees (GBT). Figure  1  shows a sketch of the different models strategies to associate a probability to a collection of input variables, in our case the probability of developing the disease starting from genotype data. All data and model analyses have been performed in Python, with an extensive use of the library Scikit-learn  37. All details concerning the parameters used in the different models can be found in the Supplementary Information text.
From a mathematical point of view, one of the key ingredient that makes NN so efficient is their ability to integrate differentiable operations well suited to the structure of the data (convolution for images, recurrent units for time series, attention mechanisms for sequences etc. ). In this paper, we used some of the latest tools from deep learning, such as residual connections  38  , batch normalization  39  , and dropout  40  , into fully connected NN. In particular, we studied separately the cases: • Dense NN with one fully connected hidden layer, but with a variable number of neurons. • Dense NN with different numbers of fully connected hidden layers, all composed by 64 neurons. • Dense NN with different odd numbers of fully connected hidden layers, all composed by 64 neurons, with full pre-activated residual blocks  41  All models have been implemented in Python using the library Keras 42 running on top of TensorFlow  43.
To implement these models, we used the corresponding Python packages. Random. In order to correctly identify the properties of the different models, we built a "random" model in which random weights are given to the preselected SNPs. In particular, this model allowed us to address the question of feature importance selection by taking into account biases related to the distribution of SNPs on the Immunochip. The AUC for this model was around 0. 5 for any fold of the dataset.
The main advantage of PFI score is that it can be universally used, independently on the ML model. Moreover, in order to distinguish the dependence of the results on the model from that due to the criterion used to assign the ranking, we considered also a different criterion for LR and LGBM models. For LR we used the absolute value of the weight associated to each feature. The higher this value after training, the most important the corresponding feature was considered. For LGBM we used the 'gain' option already implemented in the library (it measures the average gain of the feature when it is used in trees). A given model, with a fixed criterion, trained on different subsets of the data, give different results for the feature importance scores. To take into account this variability, for the second part of our analysis we created 10 different folds: starting from the original whole dataset, we arbitrarily permuted 10 times the data and re-divided the dataset in Train and Test sets (with the same proportion of samples).
Results We present here the results obtained by applying several ML classification methods to the IIBDGC Immunochip dataset, in terms of AUC scores and features selection after training. We first report the impact of QC and imputation strategies on the performances of the classification algorithms, in the framework of Lasso LR. We then show the results obtained by the same linear models under different regularizations. Next, we describe the results obtained by powerful and popular algorithms, GBT and NN, shallow and deep. Finally, the comparisons between the best features identified by each method are reported. Linear models. Different algorithms have been applied to the genotype/phenotype association problem in the literature. Besides the univariate models (based on the p-values on single SNP association tests), the simplest multivariate analyses use linear models in order to associate a phenotype to a genotype.
The unprocessed available dataset was associated with higher AUC values suggesting that artifacts may affect the results by inflating these values. The bias was likely related to the presence of missing genotypes which may reflect remaining stratification biases (data not shown). Indeed, imputing the values of the missing alleles using the HW method resulted in a large decrease of AUC values. SNP preselection with more stringent QC criteria limited the impact of biases. Also in this case, the HW method was the best choice to limit the impact of missing genotypes. Notice that there was no major effect of genotype imputation using the allele frequencies derived from the whole dataset or healthy controls only. Finally, the coding method also affected the results, the higher AUC values being observed for the sum method which consists in counting the number of rare alleles contributing to the genotype.
Interestingly, the number of retained SNPs in Lasso model depended on the processing method. While a stringent QC process did not affect significantly the number of retained SNPs, the sum coding method retained less SNPs for a better result, suggesting that it extracts more information from genotypes. Qualitatively, the retained SNPs varied a lot according to the pre-processing method. This is expected because SNPs in strong linkage disequilibrium are interchangeable if they carry a shared information. We thus looked at the loci (defined by a regions of 500 kb, see Materials and Methods section) kept by the different models. The vast majority of the predictor loci are shared by the Lasso methods when performed on the QC dataset whatever the imputation and coding strategies (Table  1 ). The loci contributing to the Lasso model can be classified according to their respective weights.
The best SNPs appeared often different from one model to another, whatever the data processing strategies and even for the same set of analyses. This finding argues for a large number of SNPs with comparable and low weights. For homogeneity reasons, in the rest of the analysis we always used the QC/HWc/sum pre-processing strategy. Feature and sample preselection. To explore the impact of SNP preselection, we next performed Lasso analyses on sets of SNP which passed at least one nominal association test with a p-value threshold ranging from 10 -8 to 1 (see Materials and Methods). The mean AUC value was not significantly affected when SNPs with p-values higher than 10 -5 were discarded (Fig. 2A ). However, even retaining only the most strongly associated SNPs (p-values lower than 10 -8 ), the AUC remained higher than 0.
In the original dataset, the case-control ratio was 0. 53. To explore the impact of this unbalanced ratio on the results of Lasso LR, we changed it with alternative values ranging from 0. 53 to 1. 5, by eliminating controls from the dataset. The results were slightly lower, arguing for keeping all the available samples in the analyses (Fig. 2C ). Due to these results, the rest of the analyses were performed on the set of all cases and controls and for preselected SNPs with >. MAF 0 01 and association p-values lower than 10 -4. The discarded features could in principle have a more important effect for non-linear models. Nonetheless, a similar analysis performed for the GBT algorithm LightGBM shown that no significant improvement of the AUC score is obtained when relaxing the constraints on preselection p-values (Supplementary Fig.
S2 ). Therefore, we made the choice to keep the same preselection conditions also for all non-linear models, for consistency and to alleviate the problems related to the curse of dimensionality. Regularization. To explore the impact of alternative regularization methods, we analyzed the dataset with L2 (Ridge) or mixed L1-L2 (ElasticNet) penalized regression methods (Fig. 2D ). AUC values obtained with these methods were very similar to those of Lasso. However, the number of SNPs contributing to the retained models increased from 2575 (Lasso) to 4245 (Elastic net) and 21896 (Ridge). This finding suggests that the alternative methods were efficient to detect additional SNPs with very small effects but this property did not change the global result. Non-linear models. Non linear methods are supposed to be more efficient in detecting non linear epistatic interactions between genotypes when compared to LR methods.
80 also in this case. Gradient boosting trees. GBT are a family of alternative methods proposed to go beyond linear additive models and take into account complex gene-gene interactions. We investigated three different state-of-the-art algorithms (XGBoost, LightGBM and CatBoost): As for NN, the mean AUC were in the range of 0. 80 (Fig. 3D ). In summary, non linear methods did not appear more performant than linear ones arguing for limited epistatic effects in the genotyping dataset. Combining the models. We investigated the possibility of combining different models in order to improve the performances. In particular we tried an ensemble method, which consists in training several classifiers and combine their predictions to check if it can outperform any single classifier. In our case, when using the average as combination rule, by combining LR, ResDN3 and LGBM we obtained =.
Feature importance comparisons between methods. Comparisons with previous GWAS analy- ses. The IIBDGC performed a large GWAS meta-analysis which included the Immunochip dataset  3. As a result, an association at genome-wide significance ( < ⋅ - p 5 10 8 ) were retained for 140 independent CD loci (a locus was defined as a genetic region of 500 kb around the best associated SNP. In this paper we use a slightly different definition, see Materials and Methods) and the corresponding SNPs. We compared the set of feature SNPs identified by Jostins et al. and those of three methods representative of the LR, NN and GBT approaches, respectively Lasso with weight as feature importance score, LGBM with gain and ResDN3 with PFI, as described in Materials and Methods. As shown in Fig. 4A , most of the SNPs with genomic nominal significance contributed to the architecture of the ML models: SNPs with the largest OR in the analysis of 3 also correspond to a peak for LR, LGBM and ResDN3 methods.
S1 ). Far from the origin though, the first more important SNPs for all models and more importantly for non-linear models like ResDN3 and LGBM, deviate from those by GWAS studies. Nonetheless, almost all CD-associated loci reported in the analysis of 3 appear in the best 800 loci for ML models, independently on the model and ranking criterion. Robustness of the ML results. We next questioned the robustness of the ML methods in determining the important features for this case/control classification problem. Therefore we looked at the robustness R of important loci arising by iterating the analyses using the same model trained on different folds of the data (Fig. 5A ). A large proportion of the first tens of loci were common to all the tests. However, the proportion of loci consistently selected as important, was higher for LGBM than for LR and ResDN3 between the first ~25 and ~100 loci.
5C ), when trained on the same fold of the dataset. As controls, we computed also the robustness between every model and the random choice of markers, and shown the mean results over the different results. Among between-methods comparisons, LR with PFI gives the more consistent ranking with LGBM with PFI for  x 175 and with ResDN3 with PFI for  x 175. On the other hand, ResDN3 provided features significantly more different from those of LR with weight and LGBM with gain. This observation suggests that either the measure used to weight the different loci in the NN models was not suitable, or the ResDN3 method was complementary to LR and/or LGBM. Joint results. Because LR with weight, LGBM with gain and ResDN3 with PFI, could be seen as complementary, we looked at their joint results in terms of feature importance selection.
First, we looked at the loci simultaneously present within the first 140 most important of the three ML models, but absent in the GWAS meta-analysis. Two new loci were identified. rs35320439, located on chromosome 2 was in the vicinity of a gene coding for galactose-3-O-sulfotransferase 2 (GAL3ST2) which is over-expressed under TNFα stimulation in goblet cell-like in vitro  50. rs395157 is located on chromosome 5, in the vicinity of a gene coding for Oncostatin M receptor (OSMR). Recently, high pretreatment expression of Oncostatin M was associated with failure of anti-TNF therapy in IBD patients  51. However, retrospectively, rs35320439 and rs395157 were found nominally associated with CD ( < - p 10 10 ) in the studied dataset with OR values of respectively 1. 158 and 1. 141. We then tested the robustness of a method able to integrate the results of two different models.
We defined a list of importance features for combined methods by giving a rank to the loci in the following way. For each selected combination of models/feature selection criteria, we assigned a rank to a given locus accordingly to its order of appearance in the intersection between the ranked list of all considered models. We then evaluated their consistency throughout the training on different folds of the dataset. The mean values and standard errors are shown in Fig. 5D. Similarly to the results shown in Fig. 5A , when considered in combination with LR with PFI, the contribution of LGBM with PFI improved the consistency for ranking loci list of length between ~25 and ~150, while that of ResDN3 for  x 150.