dataesr/scientific-paragraphs-categorization huggingface dataset:
Small M and ef_construct may lead to frequent traps in local minima and thus, low recall while large M and ef_construct may lead to slow speed without proportional improvement in recall (Supplemental Table  S2 ). Therefore, there is a tradeoff between accuracy and speed that should be evaluated first. However, for most users this task would not be necessary because they will work with pre-built databases such as those provided here. Further, the search step against these pre-build databases with query genomes of known taxonomy for evaluating recall and tradeoffs can be performed, within minutes, on any modern laptop with 5-6 GB of memory (Table  1 ). Kmer-based methods for genetic relatedness estimation such as ProbMinHash have lower accuracy between moderately-to-distantly related genomes compare to alignment-based tools (see supplement Note 4 for further discussion).
Our empirical evaluation showed that this relatedness level, for nucleotide searches, is around 78% ANI and 52% AAI for the amino-acid searches (e. g. , ProbMinHash distances do not correlate well with blast-based ANI and AAI at these levels). To circumvent this limitation, we designed a 3-step framework as part of GSearch to classify bacterial genomes that show different levels of novelty compared to the database genomes, with high accuracy. This framework included a search at the universal gene level for deep-branching genomes that are novel at the phylum level (AAI < 52%), for which searching at the entire proteome level is less accurate. Recently, methods that employ kmers that allow mismatches, that is, spaced kmers  39  , have shown promise in accurately estimating genomic relatedness even among distantly related genomes with gains in speed.
To apply spaced kmers to entire genomes, the recently developed "tensor sketch" approaches could be explored in the future to simplify the pipeline for bacterial and viral genomes  40. In the meanwhile, the ProbminHash approach, essentially a Jaccard distance estimation via MinHash-based analysis of kmers, is highly efficiently, and, importantly, can effectively deal with incomplete genomes or genomes of (drastically) different length, an known limitation of MASH-based methods  41. Comparing genomes of different length is not uncommon, e. g. , bacterial genome size can differ by more than two-fold, as can be the case between MAGs of different level of completeness or when searching a short sequence (e. g. , a bacteriophage genome) against a large genome collection (e. g. , the whole viral genome database).
Our own analysis showed that ProbMinHash is robust down to 50% completeness level (Supplemental Table  S6 ), which is also the most commonly used standard for selecting MAGs of sufficient/high quality  42. ProbMinHash is also robust with completed genomes with repeats or gene duplications due to the kmer weighting step. In general, the genome relatedness estimated, or best database matching genomes identified, by GSearch were highly consistent with blast-based AAI results or phylogenetic placement of the genome using GTDB-Tk, particularly for query genomes with close relatives in the database related at the species or genus level (Supplementary File 1, Supplementary Figure  S5 ). For more distantly related query genomes relative to database genomes, classification results of GSearch showed some differences with GTDB-Tk.
These differences were not always possible to assess further for the most correct genome placement but could be due, at least partly, to the incompleteness and/or contamination of query or/and database genomes, which renders the resulting concatenated alignment of universal genes used by GTDB-Tk unreliable  43  as only a few amino-acid positions per gene are used in the final alignment. In contrast, the AAI and ProbMinHash approaches should be more robust to changes of a small number of genes because the entire proteome is considered  17. Graph-based NNS methods achieve good performance compared to tree based and locality-sensitive hashing (LSH) methods. Building a HNSW graph relies on proximity of the database elements; so, if the distances among database elements, in our case genomes, cannot be effectively estimated via hashing algorithms, the navigation in graph will be less efficient (e.
g. , gets trapped in local minima) because the edges to choose from will not be accurate estimations of the relatedness of the corresponding genomes. This is especially problematic for highly sparse/distantly related and diverse datasets, like the viral genome database, in which two phage genomes could often share very little genomic information (kmers). This is confirmed by our own results when using nucleotidelevel search to build the viral graph. Hence, the amino acid level will be much more robust for viral genomes and is the recommended level to use. Finally, the HNSW graph, and graph-based K-NNS in general, can be further improved by adding shortcut edges and maintaining a dynamic list of candidates, compared to a fixed list of candidates by default  44. Graph reordering, a cache optimization that works by placing neighboring nodes in consecutive (or near-consecutive) memory locations, can also be applied to improve the speed of HNSW  45.
Initially, the genetic relatedness among a collection of database genomes is determined based on the ProbMinHash algorithm, which computes the normalized weighted Jaccard distance using the probminhash3a algorithm implemented in the ProbMinHash paper  1. The normalized weighted Jaccard distances are then used as input for building HNSW graphs (note that a distance computation is required only when that genome pair is required for graph building, thus GSearch avoids all vs. all distance computations). Genomes are subsequently recursively added as the nearest neighbors of each node in the built graph file with the same distance computation procedure. The built graph database file is stored on disk. Query genomes are then searched against graph database and subsequently, best neighbors are returned for classification/identification.
In this process, the best neighbor (or neighbors) is also identified based on the smallest normalized weighted Jaccard distance obtained. ProbMinHash Details of differences between ProbMinHash and traditional MinHash can be found in Supplenmentary Methods & Materials. We reimplemented the Probminhash algorithm in Rust to estimate genetic relatedness between any two genomes based on normalized vector (or sequence, kmer profile in our case), perform a greedy search on the proximity graph by comparing the query vector with database vectors under the searching measures (e. g. , cosine similarity or L2 similarity, in our case probminhash distance). Then, the most similar candidates are returned as outputs. The key point for these twostep methods is step 1, to construct a high-quality index graph, which provides a proper balance between the searching efficiency and effectiveness.
To guarantee the searching efficiency, the degree (number of maximum allowed neighbors, denoted as M) of each node is usually restricted to a small number (normally 20~200) while width of search for neighbor during inserting (denoted as ef_construct) is usually a larger number (higher than 1000) to increase the chance to find best M neighbors by increasing the diversity of neighbors due to the large number of neighbors retained. Building graph and searching query against the graph follow very similar greedy search procedures except that there is an extra reverse updating of neighbors list for each vector when inserting database vector (building), one by one, into the existing graph (Figure  1a ). The first phase of the insertion/building process starts from the top layer by greedily traversing the graph in order to find maximum M closest neighbors to the inserted element P in the layer by doing ef_construct times search (Figure  1a ).
After that, the algorithm continues the search from the next layer using the closest neighbor found from the previous layer as entry point, and the process repeats until to the bottom layer. Closest neighbors at each layer are found by a greedy and heuristic search algorithm (Figure  1b  and c ). For building, after searches are finished at the bottom layer for each inserted element, a reverse update step will be performed to update the neighbor list of each node in the existing graph while for searching this is not needed. The overall database building time complexity is O(N*log(N)), where N is the number of nodes in the graph. For searching, since there is no need to reverse update best neighbor list for each node in the graph, time complexity is (only) O(log(N)) (See Supplementary Note 3).
Details of program implementation in Rust There are 2 modules in total: tohnsw and request. Tohnsw is to build graph by gradually inserting genomes into graph while request is to query new genomes against the graph database built in the tohnsw step. Tohnsw starts from reading database genomes and generating kmer profile and sketches for distance calculation. By selecting a random genome as the first genome to insert to the graph, tohnsw module gradually add genomes to existing graph file following HNSW constructing rules mentioned above by computing ProbMinHash distance between genomes. Whenever a genome is going to be inserted into the existing graph, each genome in the graph is associated with a list that stores the M closest neighbors/genomes to the genome and the distance to these neighbors.
Request module will load the graph database and then search query genomes against it to return the best neighbors of each query, following exact the same procedure with building step without updating the database. Both tohnsw and request module are operating in parallel for high performance (see Supplementary Note 6). The GSearch software can be found here: https://github. com/jean-pierreBoth/gsearch GSearch relies on Kmerutils (https://github. com/jean-pierreBoth/kmerutils), which is a Rust package we developed to manipulate genomic fasta files including kmer string compression, kmer counting, filtering using cuckoo filter et. al. Installation guide, manual and pre-built binaries can also be found on the website. We provide static binaries on the release page for major platforms such as Linux and MacOS, with support for different CPU structures, e.
g. Intel x86_64 or ARM64. GSearch program can be run like this : 1) Build a graph database, which can be done running the following command: tohnsw -d. /GTDB_r207 -k 16 -s 12000 -n 128 --ef 1600; 2) Request neighbors of query genomes: request -b. -r. /query_folder -n 50 (--aa). Note that with the --add option in tohnsw module, genomes in the directory will be added to existing graph database, loaded from current directory, thus avoiding building graph database from the very beginning when there are only a small number of new genomes species compared to the current database. However, for larger number of new genome species, rebuild from start is suggested to be able to choose an optimal M and ef_construct to maintain high accuracy.
Prokaryotic classification pipeline The amino-acid level graph showed that closest neighbors were found, with high recall, when the query shared at least 52% AAI to its best neighbor. For more divergent genomes, showing lower than 52% AAI, whole-genome amino-acid level graph loses accuracy and we had to switch to universal, single-copy protein-coding genes. For the nucleotide-level graph, we used kmer=16 for bacteria and archaea to have high specificity for closely related database genomes (e. g. , sharing about 95% ANI). For building the whole-genome amino-acid graph, we used k=7 to have the best specificity without compromising sensitivity, which is also consistent with previous results on classification of amino acid sequences based on kmers 2. For building graph based on universal gene set, we use k=5 because of much smaller total amino acid size.
For further details on the range of kmer to use for bacteria genome and proteome, viral genome and proteome, see Supplemental Notes 2. The proteome of each genome was predicted by FragGeneScanRs v0. 0. 1 for performance purpose as opposed to Prodigal despite small loss in precision (Supplementary Table  S5 )  3. Hmmsearch in the hmmer (v3. 3. 2) software  16  was used to extract the universal gene set for bacteria and archaea genomes (universal gene graph). Note that for viral genomes, this last step was not used because there are no universal single copy genes for viral genomes. Evaluation of the speed and memory requirements for all steps mentioned above were performed on a RHEL (Red Hat Enterprise Linux) v7. 9 with 2. 70 GHz Intel(R) Xeon(R) Gold 6226 CPU. Unless noted otherwise, all 24 threads of the node are available by default.
Distributed implementation and database splitting To accommodate the increasing number of genomes that become available at an unprecedented speed in recent years and will soon reach 1 million or more, we provide an option to randomly split the database into a given number of pieces and build graph database separately for each piece. In the end, all best neighbors returned from each piece will be pooled and sorted by distance to have a new best K neighbor collection returned to the user for each query genome. We hereby prove that in terms of requesting top K best neighbors, the database split strategy is equivalent to non-split database strategy as long as the requested best neighbors for each database piece is larger than or equals to requested best neighbors in the non-split strategy. The underlying reason is that the best neighbors globally are also the best locally  4.
The database split and request will be done sequentially, on one node, without multi-node support. For now, we split GTDB database in to 5 pieces for testing purposes. In theory, a large database can be split into any pieces as long as each piece can be used to build HNSW. In practice, a reasonable way to decide on the number of database pieces to use is so that memory requirement for each piece is equal or smaller than the total memory of host machine. The database split idea has been used in several graph-based larger scale (e. g. , billions) nearest neighbor search tasks in industry  4, 5.
Species database and testing genomes for benchmarking and recall GTDB version 207 was used to build the database for bacteria and archaea genome species 6. The IMGVR database version 3, with species representatives at a ≥95% ANI, was used for for viral database building  7. For fungal genomes, all genomes downloaded from the MycoCosm project (on 24th Jan. , 2022) were used  8. The amino acid sequences of predicted gene on the genomes were obtained using FragGeneScanRs. The Universal Single Copy Gene (USCG) gene set for GTDB genomes were extracted via hmmer software. To test the performance of our pipeline, we specifically chose genomes that are not included in the GTDB database (the database was used for graph building). In particular, the bacterial/archaeal genomes, mostly MAGs, reported by Ye and colleagues  9  and Tara Ocean MAGs (total 8,466 MAGs)  10   gtdb.
msh query. msh. For AAI calculation, the corresponding script of the enveomics package  13  was used: aai. rb -1 query. faa -2 db. faa -p diamond -t 24. Hmmer was used to search for universal single copy gene against pre-built hmm profiles (120 for archaea and 122 for bacteria respectively); the profiles were obtained from the GTDB-Tk software. For viral genomes, FastANI fragment size of 1000 was used instead of 3000 while aai. rb fragment size was 500 instead of 1000 with minimal number of matches of 5. For viral genomes, MASH kmer size of 11 and 7 was used for nucleotide and amino acid levels, respectively. For fungal genomes, we use MUMMER v4. 0. 0 with default parameters for ANI calculation  14. Gene prediction for fungal genomes was performed using GeneMark-ES v2 (--fungus --ES)  15. Kmer size 21 and 11 was used for fungal genomes in MASH for nucleotides and amino acid levels, respectively.
following way: since biological species database are generally sparse because we are far away from sequencing all species in the environment and likely the existence of natural gaps in diversity, a larger top K by HNSW (e. g. , 100) compared to the value used in standard benchmark dataset will offer little, if any advantage, especially when the query genomes are relatively new, e. g. a new family compare to database genomes. Therefore, we use top 5 and 10. Top 5 and top 10 recall are calculated based on top 5 and 10 neighbors found by GSearch and the alterantive tools, and if all top 5 or 10 found by the latter tools were also in top 5 or 10 of our tool, then recall was 100%. Similarly, if only 4.
were used. We randomly selected 1000 genomes/MAGs from Ye's collection and use them as query genomes to test the performance and accuracy of GSearch. To compare with other database search tools for large database e. g. , the viral database, we compare GSearch with PhageCloud 11 , which builds a graph database based on the labels of each viral genome (e. g. , environment source) and its search algorithm is Dashing2 12. To benchmark how GSearch performs compared to ANI/AAI-and MinHash-based tools, we ran FastANI, Diamond blastp-based AAI and Mash to find the best neighbors for the same query genome dataset and evaluated whether or not the best neighbors found by GSearch were the same. FastANI parameters for the bacterial dataset were the following:.
Request/search performance on major CPU platforms for GTDB v207 database for 1000 queries. RHEL v7. 9, Linux v3. 10. 0-1160,all threads used. Parallel package was used to run multiprocess at the same time. FGSrs stands for FragGeneScanRs. Note that in practice only those genomes failed in the Request for nt step (best found is less than 78% ANI) will be used in this step. Only 100 genomes are used for testing hmmsearch because this step is for very new genomes at order level or above and we often do not have that many new genomes in a real-world dataset. Parallel Packages was used to run multiple processes of hmmsearch, one thread per process for hmmsearch.
Acknowledgment This work was supported, in part, by the  US National Science Foundation  (Award No  1759831  and  2129823  to KTK). We want to thank  PACE (Partnership for Advanced Computing Environment) at Georgia Tech  for providing computing resources. We want to thank  Kenji Gerhardt  for helpful discussions on benchmarking against traditional  ANI/AAI  based tools and  Chirag Jain  for discussions on the graph based nearest neighbor search.
were used. We randomly selected 1000 genomes/MAGs from Ye's collection and use them as query genomes to test the performance and accuracy of GSearch. To compare with other database search tools for large database e. g. , the viral database, we compare GSearch with PhageCloud 11 , which builds a graph database based on the labels of each viral genome (e. g. , environment source) and its search algorithm is Dashing2 12. To benchmark how GSearch performs compared to ANI/AAI-and MinHash-based tools, we ran FastANI, Diamond blastp-based AAI and Mash to find the best neighbors for the same query genome dataset and evaluated whether or not the best neighbors found by GSearch were the same. FastANI parameters for the bacterial dataset were the following: fastANI --ql query_path. txt --rl gtdb_path.
txt -k 16 -p 24 --minFrac 3000 -o ANI. txt. GTDB database was split into 50 subsets and each subset was parallelly run on a multi-node supercomputer to reduce memory requirement. MASH parameters were: mash sketch -a (for AA only) -k 21 (7 for AA) -s 12000 -p 24 GTDB/*. fna > gtdb. msh; mash dist -p 24 Recall of AAI-, ANI-and MinHash-based nearest neighbor searching for bacteria/archaea, fungi and viral genomes.
Facing the abundance and diversity of phages, bacteria have developed multiple anti-phage mechanisms. In the past three years, the number of known anti-phage mechanisms has been expanded by at least 5fold rendering our view of prokaryotic immunity obsolete. Most anti-phage systems have been studied as standalone mechanisms, however many examples demonstrate strains encode not one but several antiviral mechanisms. How these different systems integrate into an anti-viral arsenal at the strain level remains to be elucidated. Much could be learned from establishing fundamental description of features such as the number and diversity of anti-phage systems encoded in a given genome. To address this question, we developed DefenseFinder, a tool that automatically detects known anti-phage systems in prokaryotic genomes.
We applied DefenseFinder to >20 000 fully sequenced genomes, generating a systematic and quantitative view of the anti-viral arsenal of prokaryotes. We show prokaryotic genomes encode on average five anti-phage systems from three different families of systems. This number varies drastically from one strain to another and is influenced by the genome size and the number of prophages encoded. Distributions of different systems are also very heterogenous with some systems being enriched in prophages and in specific clades. Finally, we provide a detailed comparison of the anti-viral arsenal of 15 common bacterial species, revealing drastic differences in anti-viral strategies. Overall, our work provides a free and open-source software, available as a command line tool or, on a webserver. It allows the rapid detection of anti-phage systems, enables a comprehensive description of the anti-viral arsenal of prokaryotes and paves the way for large scale genomics study in the field of anti-phage defense.
While many families of systems have been investigated mechanistically, much remains to be uncovered about the anti-viral arsenal at the level of a strain, a species or all prokaryotes. Describing which systems genomes encode will be essential for understanding phage bacteria interactions, in a natural context. Establishing a holistic, genome centric view of the whole anti-viral arsenal of prokaryotes is currently challenging. Most studies tackling the distribution of defense systems focused on one or a few families of systems. They provide important number regarding their abundance in microbial genome. For example, in a dataset of 38 167 genomes, 4 894 encoded CBASS (13%) and 4446 retrons (11%)  6, 17. Systems described in  1  were found at frequencies ranging from 1. 8% (Kiwa) to 8. 5% (Gabija)  1  while RM and CRISPR-Cas are encoded by 74.
2%  18  and 39%  19  of genomes respectively. However, for many systems, the frequency remains to be studied. One of the reasons for such lack of knowledge can be attributed to the absence of a tool dedicated to the genomic detection of known anti-phage systems. Programs exist for the detection of specific systems such as CRISPR-Cas  [20] [21] [22] [23]  , as well as databases of anti-phage systems such as REBASE for RM  24  and databases of defense genes such as PADS  25. However, a single defense gene may not be enough to characterize a functional anti-phage mechanism, therefore it appears more relevant to search for complete systems. The lack of a tool able to find all known anti-phage systems is explained partly by the timeline of the discoveries of the novel systems (mostly since 2018), by their large number (more than 50) and the complex biology of defense systems.
In this study, we developed a tool, DefenseFinder, to detect known anti-phage systems from a genomic sequence. We used this tool, to detect all known anti-phage systems in a database of more than 21 000 complete microbial genomes, describe and analyze their distribution at different phylogenetic scales (from the genome to microbial species and phyla) to provide a systematic and quantitative view of the anti-viral arsenal of prokaryotes.