dataesr/scientific-paragraphs-categorization huggingface dataset:
, 2014; Eren et al. , 2015). In order to apply those tools, however, metagenome assembly needs to be performed, which is computationally demanding and introduces additional challenges compared to single genome assembly such as the uneven coverage of contigs (contiguous sequences) from different organisms or the presence of micro-diversity  (Quince et al. , 2017). Together with the ever-growing sequencing capacity, it becomes increasingly demanding to identify which of the available datasets (publicly deposited or locally generated sequence datasets) actually contain sequence data of a given taxon or gene of interest. Instead of assembling short reads into longer contigs, the taxonomic composition of a metagenomic or metatranscriptomic dataset can be assessed using microbiome profilers that classify reads directly.
Traditional phylogenetic tools on the other hand offer several robust evolutionary models for both nucleic and amino acids that theoretically allow for a more reliable taxonomic assignment of sequences, but are slow compared to similarity-based methods, usually prohibiting their application to large metagenome datasets. In addition, short reads generally do not provide enough phylogenetic signal, leading to artifactual inferences  (Matsen et al. , 2010). Several tools have been developed to overcome these barriers by instead placing fragmentary sequences (particularly from amplicon sequencing data) onto a phylogenetic reference tree  (Matsen et al. , 2010; Berger et al. , 2011; Barbera et al. , 2019). Shotgun metagenomic or metatranscriptomic data is often analysed with a focus on gene rather than genome reconstruction, and is then usually called gene-centric.
The goals of the work presented here were to: i Create a computational workflow that could determine the presence of taxa of interest in large short-read datasets based on genecentric assembly and robust phylogenetic inference, especially with the objective of selecting good candidate datasets for metagenomic assembly and genome-resolved metagenomics; ii Create a workflow that uses state-of-the-art methods and is versatile and fast enough to accommodate a broad range of applications, while being modular in order to easily incorporate new approaches; iii Compare the workflow's performance in terms of computational footprint and sensitivity/precision to GraftM, another recently published tool with a similar application.
Implementation PhyloMagnet exploits the idea of gene-centric assembly  (Huson et al. , 2017)  to efficiently screen sequence datasets of short reads for target genes, and to taxonomically classify assembled gene sequences using phylogenetic placement. Below is a description of the analysis steps employed by the pipeline (see also Fig. 1 ), which requires the following inputs: a One or several query short-read sequence data files in FASTQ or FASTA format (potentially 'raw', untrimmed reads, see 2. 3), corresponding to the metagenomic or transcriptomic dataset(s) to query (Fig. 1 :1); b One or several homologous groups of reference proteins, each sequence annotated with its taxonomic affiliation (in the EggNOG format, containing ncbi's taxonomy ID and a unique identifier, e. g. '70448. Q0P3H7').
Alignment to reference protein sequences Identifiers from the EggNOG database, containing orthologous groups of protein sequences from all domains of life with functional annotations  (Huerta-Cepas et al. , 2016b) , can be specified to be used as input reference sequences. Alternatively, sets of homologous protein sequences curated by the user in FASTA format can be used. In order to check for the potential presence of homologs encoding these proteins of interest in the query metagenomes or metatranscriptomes, each of the short-read datasets given as input (see 2. a above) is then aligned to the collection of reference protein sequences using the DIAMOND aligner in blastX mode (Fig. 1.
'Magnetizing' trees and identifying candidate datasets The user can choose to specify taxonomic names (e. g. 'Escherichia') that should be used to filter ('magnetize') the list of all labels, specify a taxonomic rank (e. g. 'family'), or a combination of both. The occurrences of the chosen taxonomic labels are summarized per reference group and metagenomic or transcriptomic datasets in order to assist manual decision of candidate datasets (Fig. 1 :6). The information of how many trees were positive for a taxon of interest can be used as an approximation of coverage (see 4). The user could for example select datasets that display differential coverage for subsequent genome extraction, which often relies on such differences to group genome contigs together  (Albertsen et al. , 2013; Alneberg et al.
Availability PhyloMagnet is an open source software package and released under a GPLv3 license. It is written as a Nextflow (Di Tommaso et al. , 2017) script and available on github (github. com/maxemil/PhyloMagnet). Several functions and utilities are implemented either in python or bash  (Dalke et al. , 2009; McKinney, 2010; Huerta-Cepas et al. , 2016a). All needed dependencies are available as a singularity  (Kurtzer et al. , 2017)  container (singularity-hub. org/collections/978) and the documentation can be found on ReadTheDocs (phylomagnet. readthedocs. io).
Benchmarking To evaluate the performance of the PhyloMagnet workflow and exemplify its potential uses, we performed three benchmark experiments using an in vitro mock community as well as environmental metagenomic and metatranscriptomic sequencing datasets. We chose the datasets such that we could compare the results produced by PhyloMagnet to reference genome mapping data  (Singer et al. , 2016) , genomes extracted from metagenomes with taxonomic annotation  (Delmont et al. , 2018)  and an assembled metatranscriptome  (Frazier et al. , 2017) , respectively. For details on command line parameters see the supplementary methods.
Reference sequences To assess the general taxonomic composition of datasets we used a set of 16 ribosomal proteins (rp16) that are thought to represent reliable phylogenetic markers, as they should be vertically inherited throughout evolution and present in a single copy in most organisms  (Brown et al. , 2015). For this, we downloaded the corresponding sets of unaligned homologous sequences from the EggNOG database v4. 5. 1 (Huerta-Cepas et al. , 2016b). As a second set of reference protein sequences, we used the set of 12 protein coding genes known to be present in chloroplast genomes of Dinophyceae  (Howe et al. , 2008). This phylum of single-celled algae can be found in a wide range of aquatic environment and notably contains coral symbionts within the genus Symbiodinium  (Gómez, 2012).
For each of the genes we downloaded all available curated chloroplast encoded protein sequences for all phyla from UniProt  (Apweiler et al. , 2004)  as well as all available proteins from the Dinophyceae from the same database. All reference groups were aligned using MAFFT E-INS-i (Katoh and Standley, 2013) and reference trees were reconstructed using IQ-TREE (under the LG+G+F model;  Nguyen et al. , 2015).
Datasets The first dataset we selected was the MBARC-26 (Mock Bacteria ARchaea Community), an in vitro mock community of 23 bacterial and 3 archaeal strains (in 24 genera) with finished reference genomes that were pooled and sequenced on an Illumina HiSeq instrument  (Singer et al. , 2016). As the taxonomic classification is dependent on the reference sequences, we added orthologous sequences to the EggNOG rp16 references for those genera missing from the original EggNOG datasets. To avoid using identical sequences as references, we used available genomes from related species within the same genera to expand the rp16 references. The orthologous proteins were identified by performing HMMER (v3. 1b2; Eddy, 2011) searches of the EggNOG rp16 reference alignments against the additional proteomes (see supplementary methods).
These extended rp16 references were then used as references in PhyloMagnet in order to classify the MBARC-26 short-read data. As a second dataset we used several metagenomic datasets from the geographic region 'Southern Ocean' that are part of the metagenomic datasets of the Tara Oceans Initiative, as defined by  Delmont et al. (2018). We used the EggNOG rp16 references to assess taxonomic composition in those datasets and compared the results to the taxonomic classification of the MAGs reconstructed by  Delmont et al. (2018). The authors extracted 375 genome bins from these datasets, but only presented detailed information, including taxonomy, for 13 'non-redundant' MAGs that passed several quality and completeness filters. To be able to compare our classification results to a more extensive set of extracted genome bins, we inferred taxonomic labels for those bins that were not part of the 13 non-redundant MAGs using the tool sourmash that uses k-mer matches to taxonomically classify genomes (Titus  Brown and Irber, 2016).
Finally, we analysed the metatranscriptomes published by  Frazier et al. (2017) , who sequenced mRNA from both healthy corals and such that are affected by so-called 'bleaching', a stress response in which Symbiodinium symbionts are expulsed  (Howe et al. , 2008). We used the chloroplast protein references to search for Dinophyceae (and especially Symbiodinium) sequences in the metatranscriptome data. We then compared the assembled sequences and their placement in the reference trees with the sequences from the metatranscriptomic assembly available at NCBI's GEO database  (Barrett et al. , 2012). Similar to how the assembly was generated, we combined all of the 27 individual datasets by  Frazier et al. (2017)  into a single dataset for this analysis. In order to identify the transcripts in the assembly, we performed a tblastN search, querying the reference sequences against a database of the quality filtered transcripts.
Comparison with GraftM We compared the performance of PhyloMagnet with that of the recently published tool GraftM  (Boyd et al. , 2018, v0. 11. 1). GraftM also places sequences (either unassembled reads or pre-assembled contigs) onto a reference phylogeny using the tool pplacer, for which EPA-ng represents a scalable replacement that is able to handle larger amounts of data. We created GraftM reference packages (gpkgs; containing the reference alignment, tree and the taxonomic annotation) from each of the extended rp16 references using the create command (see supplementary methods). We then used each gpkg to analyse the MBARC-26 dataset and recovered taxonomic classifications of the query sequences. For both tools, we counted the number of genera that were correctly identified in each tree (true positives) as well as the number of genera that were identified even though they were not present in the MBARC-26 mock community (false positives).
Classification of ribosomal proteins in the MBARC-26 dataset We evaluated the performance of PhyloMagnet and GraftM to detect the presence of the 24 MBARC genera (23 of those detectable, as Nocardiopsis was part of the pooled community but not present in the sequence data from  Singer et al. , 2016)  in the metagenomic dataset (Fig. S1  and Table  S1 ). The number of correctly detected as well as falsely reported genera are displayed in Fig. 2. PhyloMagnet correctly identified up to 20 (87%) of the MBARC genera and up to 7 (with an average of 2) false positive genera in all of the 16 trees. In contrast, GraftM identified a maximum of 9 (39%) of the correct MBARC genera while giving up to 14 (with an average of 4) false positives for each tree (Fig. 2 ). Some of the reported false positive and false negative errors of both PhyloMagnet and GraftM could be attributed to closely related and possibly unresolved taxonomic groups such as Escherichia/Salmonella, Thermobacillus/Paenibacillus or possibly Clostridium/Ruminiclostridium.
Another confounding factor might be the well-known disagreement between phylogeny and taxonomy in some cases (e. g. Escherichia/Salmonella; Retchless and Lawrence, 2010). Some taxa with very low abundance in the data (e. g. Corynebacterium and Clostridium) were picked up by GraftM but not PhyloMagnet, which is likely due to the fact that there are not enough reads to reconstruct longer contigs for theses taxa, impeding an identification by PhyloMagnet as we used the default cut-off implemented in the gene-centric assembler. In general we observe a correspondence between the percentage of mapped reads  (Singer et al. , 2016)  and the number of trees a genus was detected in (Fig. S2 ), suggesting that we can use the number of trees as a rough proxy for the abundance of a taxon in a dataset. When comparing results for the full dataset and the subsampled datasets, PhyloMagnet seems to profit immensely from the additional data, likely because the assembler can connect more reads and thus reconstruct more contigs above the length threshold.
In terms of runtime, PhyloMagnet is twice as fast as GraftM when using 10 threads, making more efficient use of available computational resources. It uses, however, significantly more memory due to the requirements of MEGAN that performs the memory intensive sequence assembly, which GraftM does not include (see Fig. S1 ). We tested the performance of PhyloMagnet in a scenario where no sequences from the correct genus or family were available by removing the respective sequences from the rp16 references. PhyloMagnet was then run on each such reference dataset and we assessed wether the correct family or order was recovered. Our results show that if not the right genus (or family) is present but only other sequences from the corresponding family (or order), the right taxon could be recovered in 30% and 20% or the reference trees (see Fig.
S3 ). To further compare the performance of PhyloMagnet with that of a k-mer based metagenome profiling tool, we used Kraken2  (Wood and Salzberg, 2014; Wood et al. , 2019)  to classify the MBARC-26 dataset. Kraken performs the classification at a significantly higher speed and with a very good recovery rate of true positives, but also predicts a significant amount of false positive labels (see Fig. S4 ).
Classification of taxa in the Tara Southern Ocean dataset The PhyloMagnet workflow could identify 65 taxa (families) over the 8 datasets and 16 reference trees, whereof 21 were found in at least 4 trees for at least one dataset (Fig. 3 ). These taxa cover all but one of the taxonomic groups for which genomic bins could be identified by  Delmont et al. (2018, marked  with an asterisk in Fig. 3 ). Most noticeable, the authors of that study recovered 6 non-redundant high-quality MAGs for the family Flavobacteriaceae which could be identified in every single tree for each dataset here as well. Further, the original authors identified one MAG each within the Alteromonadaceae, the Rickettsiales and the Alphaproteobacteria as well as two within Gammaproteobacteria. All of these taxonomic groups were detected by PhyloMagnet except for the Rickettsiales, who could have been misidenfied as Pelagibacterales in the phylogenetic placement.
Alternatively, the genomic bins could have been mislabelled as Rickettsiales and actually belong to the Pelagibacterales, as those two lineages commonly artifactually branch together in phylogenetic trees due to convergent genome streamlining resulting in a similar sequence composition bias  (Roger et al. , 2017; Martijn et al. , 2018; Rodríguez-Ezpeleta and Embley, 2012; Viklund et al. , 2013). It is very likely that the MAGs that were labelled as Gammaproteobacteria by  Delmont et al. (2018)  are actually members of the Piscirickettsiaceae or Porticoccaceae, which were both detected by PhyloMagnet in several individual datasets and the majority of single gene trees. Here, we also recovered the additional taxonomic labels (Porticoccaceae, Rhodobacteraceae, Pelagibacteraceae, Cryomorphaceae) that could be assigned to raw genomic bins (which were not included in the original analyses as they did not pass quality and/or completeness thresholds) from the same study (Table  S2 ).
Identification of chloroplast genes in the Coral bleaching dataset Using PhyloMagnet contigs were reconstructed from the pooled coral bleaching dataset of  Frazier et al. (2017). Using the Phylogenetic placement workflow, contigs classified as Symbiodiniaceae could be identified in 10 out of the 12 chloroplast gene reference trees. The number of contigs that were reconstructed for each gene from the pooled sequencing data of 23 datasets ranged from 2 (psbE) to as many as 169 (psbA), whereas we could identify either one or two transcripts from the corresponding published transcriptome assembly for 9 out of the 12 genes. The two genes for which no contigs could be reconstructed were psbI and petD, both missing in the assembled transcriptome as well, which is likely due to two distinct issues.
Conclusion We have shown that by applying phylogenetic placement methods to protein sequences that were reconstructed from short-read sequencing data, our PhyloMagnet workflow can accurately identify short-read sequence datasets that contain sequences for genes and taxa of interest. We compared PhyloMagnet to a similar tool that does not rely on using a gene-centric assembly approach and demonstrated that PhyloMagnet is faster and has a higher precision and sensitivity (at the price of consuming more memory). PhyloMagnet allows researchers to explore the microbial diversity of a specific clade, or to specifically assess the presence of a metabolic pathway of interest. For example, PhyloMagnet was able to identify several lineages from single gene trees that match the results of a genome-resolved metagenomic study, showcasing how our tool could be used to screen the contents of a metagenomic dataset before applying metagenome assembly and binning methods.
Acknowledgements Several computations were performed on resources provided by  SNIC  through  Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX)  under Project  SNIC 2019-3-28. (http://www. singek. eu/;grant agreement No. H2020-MSCA-ITN-2015-675752 ). LE was funded by the  European Union 's  Horizon 2020 research and innovation programme  under the  Marie Skłodowska-Curie  grant agreement No  704263. TJGE is supported by grants of the  European Research Council  ( ERC Starting  grant  310039-PUZZLE_CELL ), the  Swedish Foundation for Strategic Research  ( SSF-FFL5 ) and the  Swedish Research Council  (VR grant  2015-04959 ).
Moreover, since biology is complex, a simple statistical model such as a linear regression might not be able to find important biomarkers. Those that are found in such experiments are often hard to reproduce, suggesting overfitting. Exploring the solution space and finding true biomarkers are not only statistically challenging, but also computationally expensive. In machine learning terms, biomarker discovery can be formulated as a problem of feature selection: identifying the best subset of features to separate between categories, or to predict a continuous response. In the past decades, many feature selection algorithms that deal with high-dimensional datasets have been proposed. Due to the difficulties posed by high-dimensionality, linear methods tend to be the feature selector of choice in bioinformatics.
A widely used linear feature selector is the Least Absolute Shrinkage and Selection Operator, or Lasso  (Tibshirani, 1996). Lasso fits a linear model between the input features and phenotype by minimizing the sum of the least square loss and an ' 1 penalty term. The balance between the least square loss and the penalty ensures that the model explains the linear combination of features, while keeping the number of features in the model small. However, in many instances biological phenomena do not behave linearly. In such cases, there is no guarantee that Lasso can capture those non-linear relationships or an appropriate effect size to represent them. In the past decade, several non-linear feature selection algorithms for high-dimensional datasets have been proposed. One of the most widely used, called Sparse Additive Model, or SpAM  (Ravikumar et al.
, 2005) , they select the features with the strongest dependence with the phenotype. However, associationbased methods do not account for the redundancy between the features, which is frequent in biological datasets, since they do not model relationships between features. Hence, many redundant features are typically selected, hindering interpretability. This is important in applications like drug target discovery, where only a small number of targets can be validated, and it is crucial to discriminate the most important target out of many other top-ranked targets. To deal with the problem of redundant features,  Peng et al. (2005)  proposed the minimum redundancy maximum relevance (mRMR) algorithm. mRMR can select a set of non-redundant features that have high association to the phenotype, while penalizing the selection of mutually dependent features.
In addition, it uses an ' 1 penalty term to select a small number of features. This results in a convex optimization problem, for which one can therefore find a globally optimal solution. In practice, HSIC Lasso has been found to outperform mRMR in several experimental settings  (Yamada et al. , 2014). However, HSIC Lasso is memory intensive: its memory complexity is Oðdn 2 Þ, where d is the number of features and n is the number of samples. Hence, HSIC Lasso cannot be applied to datasets with thousands of samples, nowadays widespread in biology. A MapReduce version of HSIC Lasso has been proposed to address this drawback, and it is able to select features in ultra-high dimensional settings (10 6 features, 10 4 samples) in a matter of hours  (Yamada et al. , 2018). However, it requires a large number of computing nodes, inaccessible to common laboratories.
Since it relies on the Nystro ¨m approximation of Gram matrices  (Scho ¨lkopf and Smola, 2002) , the final optimization problem is no longer convex, and hence finding a globally optimal solution cannot be easily guaranteed. In this article, we propose block HSIC Lasso: a simple yet effective non-linear feature selection algorithm based on HSIC Lasso. The key idea is to use the recently proposed block HSIC estimator  (Zhang et al. , 2018)  to estimate the HSIC terms. By splitting the data in blocks of size B ( n, the memory complexity of HSIC Lasso goes from Oðdn 2 Þ down to OðdnBÞ. Moreover the optimization problem of the block HSIC Lasso remains convex. Through its application to synthetic data and biological datasets, we show that block HSIC Lasso can be applied to a variety of settings and compares favorably with the vanilla HSIC Lasso algorithm and other feature selection approaches, linear and non-linear, as it selects features more informative of the biological outcome.
Problem formulation Assume a dataset with n samples described by d real-valued features, each corresponding to a biomolecule (e. g. the expression of one transcript, or the number of major alleles observed at a given SNP), and a label, continuous or binary, describing the outcome of interest (e. g. the abundance of a target protein, or disease status). We denote the ith sample by x i ¼ ½x The goal of supervised feature selection is to find m features (m ( d) that are the most relevant for predicting the output y for a sample x.
Block HSIC Lasso In this article, we propose an alternative HSIC Lasso method for large-scale problems, the block HSIC Lasso, which is convex and can be efficiently solved on a reasonably sized server. Block HSIC Lasso employs the block HSIC estimator  (Zhang et al. , 2018)  instead of the V-statistics estimator of Equation (  2 ). More specifically, to compute the block HSIC, we first partition the training dataset into n=B partitions ffðx ' i ; y ' i Þg B i¼1 g n=B '¼1 , where B is the number of samples in each block. Note that the block size B is set to a relatively small number such as 10 or 20 (B ( n). Then, the block HSIC estimator can be written as HSIC b ðf k ; yÞ ¼ B n X n=B '¼1 HSIC v ðf ð'Þ k ; y ð'Þ Þ; where f ð'Þ k 2 R B represents the kth feature vector of the 'th partition. Note that the computation of HSIC v ðf ð'Þ k ; y ð'Þ Þ requires OðB 2 Þ memory space.