dataesr/scientific-paragraphs-categorization huggingface dataset:
2020). Earlier, we developed two R packages for efficiently analyzing large-scale (genetic) data, namely bigstatsr and bigsnpr  (Privé et al. 2018). We then specifically derived a highly efficient implementation of penalized linear and logistic regressions in R package bigstatsr, and showed how these functions were useful for genetic prediction with some applications to the UK Biobank  (Privé et al. 2019 ). Here we would like to come back to some statements made in  (Qian et al. 2020 ) and benchmark bigstatsr against snpnet for fitting penalized regressions on large genetic data. We re-investigate the similarities and differences between the penalized regression implementations of packages snpnet and bigstatsr. Through some theoretical expectations and empirical comparisons, we show that package bigstatsr is generally much faster than snpnet, as opposed to what is presented in  Qian et al.
e. setting lots of β j to 0, a priori  (Tibshirani et al. 2012). Therefore the cyclical coordinate descent used to solve the lasso can be performed in a subset of the data only thanks to these strong rules. However, the main drawback of these strong rules is that they 2 require checking Karush-Kuhn-Tucker (KKT) conditions a posteriori, usually in two phases. These KKT conditions are first checked in the ever-active set, i. e. the set of all variables j with β j = 0 for any previous λ. Then, the cyclical coordinate descent has to be rerun while adding the new variables that do not satisfy these KKT conditions (if any). In a second phase, the KKT conditions are also checked for all the remaining variables, i. e. the ones not in the ever-active set. This last step requires to pass over the whole dataset at least once again for every λ tested.
Then, when the available random access memory (RAM) is not large enough to cache the whole dataset, data has to be read from disk, which can be extremely time consuming. To alleviate this particular issue,  Qian et al. (2020)  have developed a clever approach called batch screening iterative lasso (BASIL) to be able to check these KKT conditions on the whole dataset only after having fitted solutions for many λ, instead of performing this operation for each λ. Hence, for very large datasets, the BASIL procedure enables to fit the exact lasso solution faster than when checking the KKT conditions for all variables at each λ, as performed in e. g. R package biglasso  (Zeng and Breheny 2017). A more pragmatic approach in bigstatsr Since our package bigstatsr is not using this BASIL procedure, reading  Qian et al.
(2020)  could give the impression that fitting penalized regression with R package bigstatsr should be slow for very large datasets. We do check the KKT conditions for variables in the ever-active set, i. e. for a (small) subset of all variables only; this first checking is therefore fast. However, KKT conditions almost always hold when p > n  (Tibshirani et al. 2012) , which is particularly the case for the remaining variables in the second phase of checking. Because of this, we decided in  Privé et al. (2019)  to skip this second checking when designing functions big_spLinReg and big_spLogReg for fitting penalized regression on very large datasets in R package bigstatsr. Thanks to this approximation, these two functions effectively access all variables only once at the very beginning to compute the statistics used by the strong rules, and then access a subset of variables only (the ever-active set).
As we will show, this means that fitting penalized regressions using the approximation we proposed in  Privé et al. (2019)  is computationally more efficient than using the BASIL procedure proposed by  Qian et al. (2020) , and yet provides as accurate predictors. Moreover, as bigstatsr uses memory-mapping, data that resides on disk is accessed only once from disk to memory and then stays in memory while there is no need to free memory, contrary to what could be read from  Qian et al. (2020). Only when the ever-active set becomes very large, e. g. for very polygenic traits, memory can become an issue, but this extreme case would become a problem for package snpnet as well. Please refer to the Discussion section of  Privé et al. (2019)  for more details on these matters. In summary, bigstatsr effectively performs only one pass on the whole dataset while snpnet performs many passes, even though the number of passes in snpnet is reduced thanks to the BASIL approach.
(2020) , namely height, body mass index (BMI), high cholesterol and asthma. We define height using field 50, BMI using field 21001, high cholesterol using field 20002 ("Non-cancer illness code, self-reported"). Asthma is defined using field 20002 as well as fields 40001, 40002, 41202 and 41204 (ICD10 codes); please see code for further details (Section "Software and code availability"). For height and BMI, L1-penalized linear regressions are fitted using function big_spLinReg from bigstatsr and using parameter family="gaussian" in snpnet. For high cholesterol and asthma, L1-penalized logistic regressions are fitted using function big_spLogReg from bigstatsr and using parameter family="binomial" in snpnet. We use sex (Field 22001), age (Field 21022), and the first 16 principal components (Field 22009) as unpenalized covariates when fitting the lasso models.
Discussion We have found the BASIL approach derived in  Qian et al. (2020)  to be a clever approach that alleviates the I/O problem of other penalized regression implementations for very large datasets. BASIL makes significant and valuable contributions to the important problem of fitting penalized regression models efficiently. However, we also find that the implementation of BASIL in snpnet is still an order of magnitude slower than our package bigstatsr, which uses a simpler and more pragmatic approach  (Privé et al. 2019 ). Hereinafter we also come back to some statements made in  Qian et al. (2020)  and provide more recommendations on how to best use penalized regression for deriving polygenic scores based on very large individual-level genetic data. This also enables us to highlight further similarities and differences between implementations from snpnet and bigstatsr.
1, 0. 1, 0. 01, 0. 001, and even until 0. 0001) for smaller sample sizes. Note that using a smaller α leads to a larger number of non-zero variables and therefore more time and memory required to fit the penalized regression. In functions big_spLinReg and big_spLogReg of R package bigstatsr, we allow to directly test many α values in parallel within the CMSA procedure. Therefore an optimal α value can be chosen automatically within the CMSA framework, without the need for more passes on the data. 6 Second, for large datasets, one should always use early-stopping. We have not found this to be emphasized enough in  Qian et al. (2020). Indeed, while fitting the regularization path of decreasing λ values on the training set, one can monitor the predictive performance on the validation set, and stop early in the regularization path when the model starts to overfit (Figure  1 ).
For large datasets, performance on the validation sets is usually very smooth and monotonic (before and after the minimum) along the regularization path, then one can safely stop very early, e. g. after the second iteration for which prediction becomes worse on the validation set. This corresponds to setting n. abort=2 in bigstatsr and stopping. lag=2 in snpnet. This is particularly useful because, when we move down the regularization path of λ values, more and more variables enter the model and the cyclical coordinate descent takes more and more time and memory. Therefore, the early-stopping criterion used in both bigstatsr and snpnet prevents from fitting very costly models, saving a lot of time and memory. Third,  Qian et al. (2020)  recommend not to use scaled genotypes when applying lasso to genetic data.
However, using scaled genotypes is common practice in genetics, and is the assumption behind models in popular software such as GCTA and LDSC  (Yang et al. 2011; Bulik-Sullivan et al. 2015). Scaling genotypes assumes that, on average, all variants explain the same amount of variance and that low-frequency variants have larger effects. Speed et al. (2012)  argued that this assumption might not be reasonable and proposed another model: E[h 2 j ] ∝ [p j (1 -p j )] ν , where h 2 j is the variance explained by variant j and p j is its allele frequency. In  Speed et al. (2017) , they estimated ν to be between -0. 25 and -0. 5 for most traits. Note that scaling genotypes by dividing them by their standard deviations SD j as done by default in bigstatsr assumes ν = -1 while not using any scaling as argued by  Qian et al.
It is true that bigstatsr uses an intermediate format, which is a simple on-disk matrix format accessed via memory-mapping. However, package bigsnpr provides fast parallel functions snp_readBed2 for converting from '. bed' files and snp_readBGEN for converting from imputed '. bgen' files, the two formats used by the UK Biobank. For example, it took 6 minutes only to read from the UK biobank '. bed' file used in this paper. We then used function snp_fastImputeSimple to impute by the variant means in 5 minutes only, which is also the imputation strategy used in snpnet. When reading imputed dosages instead, it takes less than one hour to access and convert 400K individuals over 1M variants using function snp_readBGEN with 15 cores, and less than three hours for 5M variants. When available, we recommend to directly read from '.
bgen' files to get dosages from external reference imputation. As for package snpnet, it uses the PLINK 2. 0 '. pgen' format, which is still under active development (in alpha testing, see https: //www. cog-genomics. org/plink/2. 0/formats#pgen). This format is not currently provided by the UK Biobank, and can therefore be considered as an intermediate format as well. Fifth,  Qian et al. (2020)  point out the lack of flexibility of bigstatsr because it handles linear and logistic regressions but not Cox regression, and can only use standardized variables. Although Cox regression is not yet implemented in bigstatsr, it is one of the regression for which it is easy to derive strong rules  (Tibshirani et al. 2012). As for scaling, we have seen how to effectively get unstandardized genotypes using different penalty factors, and we actually recommend to use an in-between scaling, or to test different scaling values within our CMSA framework using the new parameter power_scale in bigstatsr.
Software and code availability All code used for this paper is available at https://github. com/privefl/paper2-PRS/tree/ master/response-snpnet/code. R packages snpnet, bigstatsr and bigsnpr can be installed from GitHub. A tutorial on fitting penalized regressions with R package bigstatsr is available at https:// privefl. github. io/bigstatsr/articles/penalized-regressions. html.
However, numerous molecular studies have revealed that the mitochondrial tree may deviate from the species tree. Indeed, the maternal inheritance of the mtDNA genome can be misleading for species delimitation in mammals, because females and males have usually different dispersal behaviours (female philopatry versus male dispersal)  [5, 6] , and because interspecific hybrid females are generally fertile, whereas hybrid males are often sterile (Haldane's rule), facilitating mitochondrial introgression between closely related species  [7] [8] [9]. To overcome these limitations, most recent taxonomic studies dealing with the delimitation between cryptic mammal species have focused on multi-locus datasets  [10] [11] [12] , as the use of multiple independent DNA markers has been shown to provide a strong and reliable signal for deciphering relationships among closely related taxa  [13, 14].
However, interpreting the results from multi-locus datasets can be difficult, especially when the DNA markers show low genetic variation or conflicting relationships between them. These difficulties have led to the development of a plethora of new methodological approaches for multi-locus species delimitation  [15, 16] , which may be subdivided into three categories: (1) phylogenetic methods, (2) multispecies coalescent (MSC) approaches, and (3) population genetic methods (Table  1 ). Phylogenetic methods were not originally developed for studying species delimitation, but the species monophyly criterion has been widely used since the origin of molecular taxonomy  [17]. For multi-locus datasets, several phylogenetic approaches can be considered: the concatenation of all markers into a supermatrix (although this approach has been widely criticized  [18] ), the separate analyses of the markers, or more sophisticated methods, such as � BEAST  [19]  or Super-TRI  [20].
Based on the coalescent theory, some authors have suggested that species can be delimited without monophyletic gene trees  [21]. The incorporation of the coalescent model  [22]  in certain software (e. g. � BEAST  [19] , BPP  [23]  and STACEY  [24] ) enabled the inference of species limits from multi-locus data by accounting for incongruences among gene trees in the presence of incomplete lineage sorting  [19]. MSC approaches often require prior assignments of samples to populations or taxa and are hence restricted to the validation of proposed delimitations  [25]. Population genetic approaches are generally applied to detect "cryptic substructure" between groups showing very similar phenotypes. The program STRUCTURE  [26]  is probably the most popular approach for Bayesian clustering using multi-locus data.
Nuclear and mitochondrial datasets used for the analyses Seven giraffe datasets were generated for our analyses using the sequences available in the NCBI nucleotide database: of certain genes, from the genome of Bos mutus available on NCBI (SAMN08580377); and (ii) Ovis canadensis, for which the sequences were extracted by BLAST from the genome available on NCBI (CP011888. 1) (length = 17276 nt); 5. the nuDNA-G137 dataset, comprising the alignments of original consensus sequences of 21 introns for the 137 giraffes (length = 16966 nt), which were recovered by detecting heterozygous sites in Geneious R10 (Biomatters, Auckland, New Zealand); 6. the nuDNA-G137O3 dataset, in which the nuDNA-G137 dataset was aligned to the three outgroup species mentioned above (length = 17276 nt); 7. the nuclear haplotype dataset, named nuDNA-GH274, which was inferred from the nuD-NA-G280 dataset in DNASP v5.
0  [45]  using the "Generate Haplotype Data File" option under non-consideration of gaps/ missing data and subsequent exclusion of outgroup taxa (length = 1362 nt; it contains only the sites found to be variable between haplotypes). The datasets were aligned automatically using the MAFFT algorithm implemented in Geneious R10 (Biomatters, Auckland, New Zealand) and subsequently verified by eye. All alignments generated for this study were deposited in Open Science Framework at https://osf. io/9wv86/.
Phylogenetic analyses The mtDNA-GH82O3 and nuDNA-G137O3 datasets were analysed with probabilistic methods. Bayesian inferences were conducted in MrBayes v3. 2. 6  [46]  by calculating the posterior probabilities (PP) after 10 7 Metropolis-coupled MCMC generations with tree sampling every 1000 generations and a burn-in of 25%. Maximum Likelihood (ML) analyses were performed with PhyML v3. 1  [47]  and Bootstrap percentages (BP) were calculated after 1000 replicates. The GTR+I+G substitution model was applied for both methods, as suggested by the Likelihood calculations in jModeltest  [48]  based on the Akaike information criterion  [49]. Bayesian analyses were also performed for each of the 21 introns using the model of DNA substitution selected under jModeltest (Table  2 ).
STRUCTURE analyses Giraffe haplotypes were reconstructed from the nuDNA-G137O3 dataset for each of the 21 introns by applying the PHASE v2. 1 algorithm implemented in the software DNASP v5. 0  [45] , allowing for recombination and reducing the output probability threshold of conserved regions (CT) from 0. 9 by default to 0. 6. Bayesian analyses of genetic admixture were run in STRUCTURE v. 2. 3. 4  [26]  to identify genetically homogeneous groups of individuals (populations of origin, K). The haplotype information deduced by PHASE for each of the 21 introns was used to code individuals sharing the same genotype with a unique integer in the input file  [51]. The analyses were done as recommended by Gilbert et al. [52] , i. e. , number of MCMC generations = 200 000 and burnin = 100 000 generations for K = 1-10 clusters.
[57] , which recognizes the most likely number of distinct clusters by the largest ΔK value, calculated by the rate of change in the log probability of data between successive K values; and (2) the "plateau "method of Pritchard et al. [51] , where the log probability of the data (ln Pr (X|K) was plotted against a range of K values, and the optimal K was selected as the point at which the plot curvature plateaus. A regression curve and gridlines were added to the diagrams generated by STRUCTURE HARVESTER to help in determining the point of plateau. To assess the reliability of the results, CLUMPAK  [58]  was used to display the barplots from K = 1 to 10 for each of the 20 iterations by means of the implemented software DIS-TRUCT  [59].
Analyses of nuclear haplotypes The nuDNA-GH274 dataset (nuclear haplotypes inferred from the concatenated matrix of 21 introns and 137 giraffes) was used to construct a median-joining network  [60]  using PopART v1. 7  [61]. The robustness of haplotype clusters was evaluated by bootstrapping (1000 replicates) under PAUP � v4b10  [50]  using either the Maximum Parsimony (MP) method (heuristic search, faststep option) or the Neighbor-Joining (NJ) method (GTR+I+G model), and under the Maximum Likelihood (ML) criterion using RAxML on CIPRES  [62]  (http://www. phylo. org). PopART was also used to construct a median-joining network for each of the 21 introns. For six introns (i. e. CTAGE5, NUP155, OTOF, PLCE1, RASSF4 and SOS1), individuals with missing allele(s) were not considered in the analysis in order to avoid any distortion of the results.
Allele sharing among individuals was investigated for each of the 21 introns using the haploweb approach  [30]  by building raw haplowebs with the online tool HaplowebMaker (https:// eeg-ebe. github. io/HaplowebMaker/). Haplowebs were constructed without the consideration of singletons (alleles encountered only once in the whole dataset, "remove singletons" option) and default settings for the other parameters. In the networks, curves are illustrating alleles found co-occurring in heterozygous individuals. Thereby, a group of alleles linked together by heterozygotes represents an exclusive allele pool, the corresponding groups of individuals is called a field for recombination (FFR)  [63]. Based on the allele-sharing information per marker provided by HaplowebMaker, a conspecificity matrix  [31]  was built using the online tool CoMa (https://eeg197ebe.
github. io/CoMa/) with no calculation on heterospecific pairs (option 1), i. e. no value was given to absence of sharing. In this matrix, the conspecificity score for node recovered with significant support in the Bayesian analysis (PP � 0. 9), as well as for other nodes discussed in the text, the two values above indicate the Posterior Probability with MrBayes (PP) and the Bootstrap Percentage obtained from the Maximum Likelihood analysis (BP). The three values below were obtained from the SuperTRI analyses of the 21 introns: from left to right: Supertree Bootstrap Percentage (SBP), Mean Posterior Probability (MPP) and the number of markers supporting the node (NRep). The symbol ''-" indicates that the node was not found monophyletic in the analysis, and the letter ''X" indicates that an alternative hypothesis was supported by SBP > 50.
Multispecies coalescent analyses Three coalescent-based approaches were applied to infer species boundaries within the genus Giraffa: (1) the "Species Tree Ancestral Reconstruction" template ( � BEAST  [19] ), (2) the extension of the � BEAST model called "Species Tree and Classification Estimation, Yarely" (STACEY)  [24] , and (3) the Bayesian Phylogenetics and Phylogeography program (BPP v. 3. 2  [23, 64] ) (see specifications for each program in Table  1 ). We estimated the species-tree phylogeny using the coalescent algorithm implemented in BEAST v. 2. 4. 4  [65]  in order to consider an alternative to the traditional concatenated phylogenetic approach (see Kubatko and Degnan  [66]  for caveats concerning concatenation). Inferences were based on the nuDNA-G274O6 dataset using an a priori assignment at the level of individuals, i.
e. by assigning for each of the 137 giraffes two alleles for the 21 introns. We assumed an uncorrelated lognormal molecular clock for all 21 loci. For each marker, we selected the best suited substitution model inferred in jModeltest  [48]  (Table  2 ). Analyses were run with 2x 10 8 generations, with trees sampled every 5000 steps. The. log files were analysed with Tracer v1. 7  [67]  to assess the convergence of model parameters (effective sample size [ESS] > 200). The species tree was summarized as a Maximum Clade Credibility tree in TreeAnnotator v. 1. 10 [68]  after discarding 25% as burn-in. The nuDNA-G274 and nuDNA-G274O6 datasets were further used for species delimitation analyses using the STACEY template implemented in BEAST v2. 4. 4. STACEY represents an improvement to the DISSECT model  [69]  to infer a "species or minimal clusters tree" (SMC) under the birth-death-collapse tree prior and without the requirement of a guide tree.
0 for the tree prior, a lognormal distribution with a mean of 4. 6 and a standard deviation of 2 to the growth rate prior and a uniform distribution for the relative death rate prior with a lower bound of -0. 5 and an upper bound of 0. 5. The dataset was partitioned by the 21 genes, with independent strict clock models and individual assignment of the best suited substitution model to each gene (Table  2 ). Each analysis was run for 2. 5 x 10 8 generations and convergence of parameters was assessed in Tracer v1. 7  [67]. Subsequently, the most supported number of distinct clusters was estimated using SpeciesDelimitationAnalyser v1. 8. 0  [24]  by analysing the species trees with a burn-in of 25% and the default collapse height of 0. 0001. Species delimitation analyses with BPP v3. 2 were based on a reduced dataset comprising only 66 giraffes due to software limitation.
Seventy-one individuals were excluded from the original dataset using the three following criteria: (1) 14 individuals with missing data, (2) 39 individuals sharing the same haplotype and (3) 18 individuals characterized by a long terminal branch in the Bayesian tree. We analysed the support for five taxonomic hypotheses: two species, with two possible geographic patterns (2Sa and 2Sb hypotheses), three species (3S hypothesis), i. e. G. camelopardalis sensu stricto A, G. giraffa and G. tippelskirchi, four species (4S hypothesis), i. e. G. camelopardalis sensu stricto B, G. giraffa, G. reticulata, and G. tippelskirchi, and five species (5S hypothesis), i. e. G. camelopardalis sensu stricto C, G. giraffa, G. peralta, G. reticulata, and G. tippelskirchi (see results for more details). First, we applied the A00 algorithm, the simple MSC model with the species tree fixed to explain the acceptance proportions of MCMC moves  [23]  under the default gamma prior values G  (2, 2000)  for the tau (τ, root divergence time) and theta (θ, genetic difference among taxa).