dataesr/scientific-paragraphs-categorization huggingface dataset:
doe. gov/data-and-tools/software-tools/bbtools/) or specialized algorithm such as SortMeRNA  (22)  to align reads against a representative set of rRNA database. Yet, the major concern remains that unassembled shortlength sequences do not contain all the informative regions of the 16S rRNA genes, thus reducing the accuracy of taxonomic assignments to species level  (23). Indeed, depending on the 16S rRNA region carried by each read, taxonomical level of affiliation may vary from phylum to species, impairing the true representation of the microbial community diversity. Converting short reads into full-length or nearly full-length rRNA gene sequences typically yields a more detailed taxonomic resolution, at the species or even strain level  (24) , providing a unique opportunity to thoroughly analyse microbial species that have never been identified before.
Pipelines, including filtering, reconstruction of full-length SSU rRNA genes and their classification, may represent useful tools to facilitate data analysis for scientists. PhyloFlash  (25)  is a compilation of tools to rapidly filter and reconstruct the SSU rRNAs and explore phylogenetic composition from metagenomics or metatranscriptomics datasets. To filter SSU reads from Illumina metagenomics dataset, phyloFlash proposes the generalpurpose mapper BBMap and the rRNA-specialized tool SortMeRNA. To assemble full-length SSU rRNA sequences, the general-purpose genomic assembler SPAdes  (26)  and the reference-based assembler EMIRGE  (27)  are proposed. Taxonomic identification of assembled sequences is done using VSEARCH  (28). Efficiency of 16S rRNA sequence reconstruction is the pivotal step of this approach.
The efficiency of RiboTaxa was determined using a synthetic complex community of human gut microbiota and compared to existing 16S-based metagenomics classifiers (i. e. Centrifuge  (32) , Kraken2, METAXA2, phyloFlash, SPINGO, BLCA, MEGAN). Later, RiboTaxa was applied to various real metagenomics datasets (i. e. ocean, soil, human gut) but also to gene capture by hybridization datasets from previous works to reveal novelties not detected by current bioinformatics pipelines.
Tool benchmarking on microbial MOCK community Microbial MOCK community. The microbial MOCK community was composed of 21 bacterial and 7 archaeal species (Supplementary Figure  S1 ), which was shotgun sequenced in paired-end (2 × 300 bp) MiSeq runs (Illumina) by  Gasc and Peyret (33). Microbial strains were selected from the Leibniz Institute DSMZ collection provided the availability of their genome in GenBank database. The abundance of each species was then defined based on the 16S rRNA copy number per genome and the number of genomes in the mixture. Raw Illumina reads from the shotgun sequencing (under accession number SRR5381736) were downloaded from NCBI open access Sequence Read Archive (SRA) using NCBI SRA Toolkit v2. 9. 1 (https://www. ncbi. nlm. nih. gov/ sra). Fastq-dump v2. 8. 2 (https://www.
Database preparation and indexing. Unless specified otherwise, SILVA small subunit (SSU) ribosomal RNAs (16S/18S) 138. 1 released on August 27, 2020 (http:// www. arb-silva. de) was used for rRNA read extraction, full-length 16S rRNA gene reconstruction and taxonomic classification. SILVA 138. 1 database (Fasta format) was clustered at 97% (NR97) with VSEARCHcluster fast v2. 7. 0, converted from RNA to DNA alphabet, and any ambiguous bases were replaced by random base characters using fix nonstandard chars. py. Specific database indexes were then built for each tool used in this study. The NR97 filtered and fixed database was indexed for BBMap with bbmap. sh (https://jgi. doe. gov/dataand-tools/software-tools/bbtools/), for EMIRGE  (27)  and MetaRib  (30)  with Bowtie (29) and for SortMeRNA  (22)  and MATAM  (31)  with matam db preprocessing.
py, for Kraken2/Bracken (34) with kraken2-build (kmer = 35) and for Centrifuge  (32)  with centrifuge-build. SSU rRNA reads extraction. BBMap v38. 87 (https://jgi. doe. gov/data-and-tools/software-tools/bbtools/): The input reads were aligned (mapped) against the indexed SILVA 138. 1 database, with minimum identity 70% by default, retaining all ambiguous alignments if there are multiple best scoring. Output was written in Fastq format, retaining all read pairs where at least one read could be aligned. SortMeRNA v2. 1b  (22) : Since SortMeRNA accepted only one file as input for the reads, forward and reverse paired-end reads were interleaved into a single Fastq file using merge-paired-reads. sh (SortMeRNA v2. 1b). The reads were then aligned by SortMeRNA against the indexed SILVA 138. 1, using min lis of 2 (number of candidate alignments) and E-value cut-off of 1 by default.
Aligned and unaligned reads were reported with options -aligned and -other, respectively and were written in Fastq format. To keep the order of paired-end reads -paired-in option was used and an overall statistics file (-log) was generated. After filtering, the aligned Fastq file containing the SSU rRNA reads was separated into forward and reverse paired-end reads using unmerge-paired-reads. sh (SortMeRNA v2. 1b). Assembly of nearly full-length SSU rRNA sequences. EMIRGE v0. 61. 1  (27) : Two available scripts were tested, namely Emirge. py and Emirge amplicon. py designed for metagenomics data and PCR amplicon data respectively. Both scripts were tested with the same following parameters. Average insert size of paired-end input was estimated from the initial mapping with Bowtie using mean size.
py (https://gist. github. com/timoast/ af73c0e9fac00187ee49). EMIRGE was run with 120 iterations with -max read length = 300, -insert mean = 500, -insert stddev = 100 and -join threshold = 1. Once all the iterations were completed, emirge rename fasta. py (EMIRGE) was run on the last iteration folder to convert Bam files into a single Fasta file containing all the reconstructed SSU rRNA genes and their relative abundances. During sequence assembly, the join threshold parameter influences sequence reconstruction as if two candidate sequences share ≥ join threshold value over their bases with mapped reads, then both sequences are merged into one for the next iterations. Emirge amplicon. py was run using the same above parameters with -join threshold = 0. 97 (as proposed in the original publication of EMIRGE) and the results were compared.
Identification of 16S variants present in the microbial mock community. To build a 16S rRNA database of the microbial mock community, 16S rRNA sequences were downloaded from the NCBI (https://www. ncbi. nlm. nih. gov/gene/) for each microbial strain and concatenated into a single Fasta file (excluding low-quality variants containing Ns in sequence). The sequences were first deduplicated to get rid of strictly identical sequences and then clustered at 100% using VSEARCH v2. 7. 0. This 16S mock database was used in downstream affiliation of reconstructed sequences using blastn of NCBI BLAST + v2. 11. 0  (21)  with an identity cut-off fixed at 99% to differentiate between microbial variants. Only the top hit alignment between the subject and the query sequence was targeted using HSP = 1 parameter.
Taxonomic affiliation. To classify nearly full-length reconstructed SSU sequences, Kraken2 v2. 0. 8-beta  (14) , Centrifuge v1. 0. 3-beta (32), mothur v1. 33 (classify. seqs) (35) , SPINGO v1. 3  (16) , BLCA v2. 1  (12) , QIIME 2's Sklearn classifier (classify-sklearn plugin, QIIME2 v2020. 8) (36)  and RDP Classifier v11. 5 (http://rdp. cme. msu. edu/) (37)  were used. Kraken2 and Centrifuge had their own 16S databases (SILVA 138. 1 SSU NR99 sequence set) from which indexes were created. Similarly, for SPINGO, the SILVA 138. 1 SSU database was indexed using the provided makefile (https://github. com/tiendu/SPINGO updated/ tree/master/SPINGO SILVA/database). Unlike other tools, BLCA does not require indexed database files. Instead, the SILVA SSU 138. 1 database was provided to BLCA as a BLAST formatted library obtained from the Fasta file using the makeblastdb utility from NCBI BLAST + v2.
11. 0, together with a taxonomy file using write taxonomy. py (https://github. com/FOI-Bioinformatics/flextaxd/tree/ master/flextaxd/modules). Pre-trained full-length SSU SILVA 138. 1 database was downloaded for mothur (https://mothur. org/wiki/silva reference files/) and QIIME 2's Sklearn classifier (https://docs. qiime2. org/2020. 11/ data-resources/) and for RDP Classifier v11. 5, pre-trained full-length 16S rRNA training set 18 released on 14 August 2020 (https://sourceforge. net/projects/rdp-classifier/files/), was used. Taxonomic classification using all the above tools was done with default parameters and a confidence cut-off set at 0. 7.
Tool benchmarking on a synthetic complex human gut microbiota Synthetic complex human gut microbiota. To mimic real dataset analysis, tool benchmarking was then done on a synthetic complex human gut community comprising of 100 microorganisms as described by  Lu and Salzberg (38). To simulate this microbial community, 100 genomes (complete genomes, scaffolds, contigs) in Fasta format were downloaded from NCBI. Barrnap v0. 9 (https://github. com/ tseemann/barrnap) was used to extract all SSU rRNA sequences present in each genome using default parameters. They were clustered at 100% using VSEARCH to remove duplicates and used as references. ART simulator v2. 5. 8  (39)  was then run on each individual genome file to generate 20X coverage of synthetic paired-end Illumina reads with length 250 bp using the following parameters: art-illumina -p -ss MSv3 -f 20 -l 250 -m 300 -s 100.
To calculate theoretical abundance of each microorganism, the number of 16S rRNA reads was filtered from individual paired-end file using SortMeRNA and was divided by the sum of the 16S rRNA reads present in the community of 100 microorganisms. Finally, the synthetic metagenomics dataset was produced by concatenating all forward and reverse Fastq files into a single forward (HG R1) and a single reverse (HG R2) Fastq file respectively. The paired-end human gut files were trimmed using bbduk. sh (k-mer = 21, https://jgi. doe. gov/data-andtools/software-tools/bbtools/, BBTools v37. 62) as described above. For SSU rRNA sequence reconstruction, MATAM, emirge. py and emirge amplicon. py were provided with SSU reads filtered using SortMeRNA while MetaRib was provided with unfiltered high-quality reads.
RiboTaxa pipeline RiboTaxa pipeline is written in bash and is available on GitHub (https://github. com/oschakoory/RiboTaxa). It can easily be installed using miniconda3 (https://docs. conda. io/ en/latest/miniconda. html). Running conda virt env. sh will install all the necessary tools required by RiboTaxa to analyse raw Illumina metagenomics reads. To avoid conflicts between required dependencies and those in existing environment, RiboTaxa uses a virtual conda environment. RiboTaxa pipeline includes SortMeRNA v2. 1b and EMIRGE v0. 61. 1 tools, both of which need indexed databases of their own. For that, parameters for each tool are predefined in the configuration file, indexDB arguments. conf, of RiboTaxa and in-dexDB RiboTaxa. sh is run to create all necessary database indexes required by RiboTaxa to analyse raw metagenomics sequences.
Any SSU database in Fasta format can be used, however, we recommend using SILVA SSU database as it is the most updated. RiboTaxa pipeline is illustrated in Figure  1. The inputs for RiboTaxa are shotgun metagenomics singled-end or paired-end files, which have been generated by an Illumina sequencer. The sequence data can either be in uncompressed Fastq format or in compressed Fastq. gz format. Several samples can be handled in the same folder. RiboTaxa will extract the name of the sample from the name of the Fastq file. By running Pipeline RiboTaxa. sh, RiboTaxa starts by removing Illumina adapters and trims Illumina reads using bbduk. sh (BBTools v37. 62). Quality control is performed before and after trimming with FastQC v0. 11. 9 and summary statistics are reported in a standalone html file using MultiQC v1.
11 (https://github. com/ewels/MultiQC/ releases/tag/v1. 11). For SSU rRNA gene reconstruction, RiboTaxa uses MetaRib (version from 13 November 2019) and EMIRGE v0. 61. 1 (Emirge amplicon. py). MetaRib uses unfiltered high-quality reads for rRNA genes assembly while EMIRGE uses SSU reads that have been extracted using the rRNA-specialized tool SortMeRNA v2. 1. The reconstructed sequences are then used as reference onto which unfiltered high-quality reads are mapped using BBmap v38. 87 to calculate the relative abundance. Prior to taxonomic affiliation, the reconstructed rRNA gene sequences obtained from MetaRib and EMIRGE are clustered at 97% using VSEARCH v2. 17. 0 and abundances from BBmap are summed up for sequences sharing the same cluster. Finally, the clustered sequences are classified taxonomically using QIIME 2's classify-sklearn plugin.
METAXA2 v2. 1. 2 (19) : The function metaxa2 was run using default settings with option -t = bacteria, archaea to filter SSU reads of bacteria and archaea only and to discard any other SSU reads (e. g. mitochondrial, chloroplast). Classifying filtered sequences at different taxonomic levels was done using metaxa2 ttt with default parameters and sequence abundance was estimated using metaxa2 dc. Filtered SSU reads were written in Fasta format while taxonomic identification and read counts were written in txt format. Relative abundance was calculated by dividing each taxon count by the total number of reads classified at the species or genus level. BLCA v2. 1  (12) : Taxonomic classification was performed by running 2. blca main. py with the formatted database and taxonomy file on the interleaved paired-end HG files, in Fasta format (reformat.
sh from https://jgi. doe. gov/dataand-tools/software-tools/bbtools/, BBTools v37. 62). SPINGO v1. 3  (16) : SPINGO was run on the interleaved paired-end HG files, in Fasta format, using default param-eters (k = 8, bootstrap = 10). The output taxonomy file was then summarised using the provided python script, spingo summary (https://github. com/GuyAllard/SPINGO) MEGAN v6. 21. 0  (11) : For taxonomic profiling, MEGAN takes as input DAA (DIAMOND alignment archive) formatted files. Thus, to generate DAA files, DIAMOND (40) v2. 0. 14 was used. First, the latest version of NCBI nr protein database was downloaded (on 19 July 2022) from https://ftp. ncbi. nlm. nih. gov/blast/db/FASTA/nr. gz and indexed using diamond makedb with default parameters. The paired-end HG files were interleaved into a single fastq using reformat.
sh (BBTools v37. 62) and reads were aligned against the pre-built nr database using diamond blastx, specifying the parameter -f 100 to output alignments in DAA format. Finally, the meganizer tool, daa-meganizer, was applied to the DAA file for taxonomic binning based on the SILVA taxonomy mapping file obtained from https://software-ab. informatik. uni-tuebingen. de/ download/megan6/welcome. html. Relative abundance was also calculated using the number of aligned reads output by DIAMOND. phyloFlash v2. 0 (25): Since phyloFlash offers two tools for rRNA reads extraction (BBMap or SortMeRNA) and full-length sequence reconstruction (SPAdes or EMIRGE), we tested several associations between these tools. Default tools (BBMap + SPAdes) were used to produce phyloFlash (BS) results. PhyloFlash was also run using SortMeRNA and EMIRGE (phyloFlash (SE)) or using BBMap and EMIRGE (phyloFlash (BE)).
For all conditions, cluster identity was set at 100%, the maximum length of reads was 300 bp and the taxonomic classification was set at genus and species level. Output from SPAdes and EMIRGE were written in Fasta. Abundance estimates and summary statistics were reported in a standalone Html file. All reconstructed sequences were filtered using bbduk. sh (BBTools v37. 62) to discard sequences below 800 bp. For taxonomic affiliation, phyloFlash classified the last common ancestor (LCA) consensus using the SILVA taxonomy and species-level classification was chosen. RiboTaxa was run on raw metagenomics reads using default parameters except for the following: (MOCK:max read length = 300, -insert mean = 192, -insert stddev = 98; synthetic: -max read length = 250, -insert mean = 185, -insert stddev = 90), which exclusively depend on the sequencing length of the input datasets.
Evaluation of RiboTaxa on real datasets Metagenomics and 16S gene capture by hybridization are two well-adapted techniques to target the phylogenetic 16S marker gene. For RiboTaxa evaluation, metagenomics datasets (paired-end) as well as 16S gene capture by hybridization datasets (paired-end) from different environments including soil, ocean and human gut (Table 1) were downloaded from NCBI SRA using NCBI SRA Toolkit v2. 9. 1 (https://www. ncbi. nlm. nih. gov/sra) as described earlier. The aim was to evaluate the versatility of RiboTaxa on different environments and on different techniques. All the parameters of RiboTaxa were kept as default except parameters associated with sequencing length of the input data (Ocean samples: -max read length = 101, -insert mean = 121, -insert stddev = 96; Human gut samples: -max read length = 151, -insert mean = 203, -insert stddev = 131; 16S gene capture by hybridization samples: -max read length = 151, -insert mean = 233,insert stddev = 127).
RiboTaxa results (reconstructed sequences and taxonomies) for each dataset were analysed to measure any potential phylogenetic signals of interest present in the reconstructed SSU sequences that were probably missed by initial studies. For the sequence of interest, a similarity search was conducted using blastn of NCBI BLAST + v2. 11. 0  (21)  with default parameters against the GenBank database. Phylogenetic analysis was performed using the pipeline phylogeny. fr (41). The 16S rRNA gene sequences in Fasta format were aligned with MUSCLE (42), followed by a curation step using Gblocks  (43). A phylogenetic tree was reconstructed with PhyML by using the maximum-likelihood method  (44). For each phylogenetic dataset, a percent identity matrix was also generated using Clustal Omega (45) and we used the identity thresholds defined by  (46)  to determine novel taxa, that is, 97 % for species, 94.
SSU rRNA gene reconstruction versus metagenome-based approaches To compare between SSU rRNA genes reconstruction by RiboTaxa versus metagenome-based approaches, taxonomic classification based on GTDB database (47) was also performed using Kraken2/Bracken. For this, 10 metagenomics samples from real datasets (3 from octocorals, 3 from semi-supercentenarians, 4 from N160 permafrost soil) were selected in which potentially new genus/species were identified. The aim was to see whether GTDB-based taxonomy could provide precise identification of novel species in different environments (i. e. ocean, human gut, soil) when there are few or no available genome references. The pipeline, struo  (48) , was used to build the index files of GTDB r202 for Kraken2 and Bracken. Kraken2 was then run on the metagenomics samples using default parameters to generate reports files which were then used by Bracken to calculate relative abundances.
rRNA reads extraction and nearly full-length 16S rRNA gene reconstruction using the MOCK community To extract SSU reads from metagenomics data of the MOCK community, we compared the widely used generalpurpose mapper, BBMap and the rRNA-specialized tool SortMeRNA. Even though BBMap was faster than SortMeRNA, it filtered only 6339 pairs of SSU reads (0. 55%) in contrast with SortMeRNA which filtered 10 675 pairs of SSU reads (0. 93%) from the raw metagenomics dataset. Then we compared the efficiency of four reference-based targeted assemblers (Emirge amplicon. py, Emirge. py, MATAM and MetaRib) using three input datasets: unfiltered metagenomics reads, BBMap-filtered reads or SortMeRNA-filtered reads (Figure  2 ), resulting in 12 combinations. Emirge amplicon. py reconstructed the highest number of sequences, followed by MetaRib, Emirge.
py and finally MATAM. Emirge amplicon. py output 54 sequences (375-1410 bp) from unfiltered reads, 50 sequences (497-1495 bp) from SortMeRNA-filtered reads and 42 sequences (437-1457 bp) from BBMap-filtered reads (Figure  2 ). In all cases, more 16S rRNA sequences were reconstructed from unfiltered reads followed by SortMeRNAfiltered reads and BBMap-filtered reads. The highest median length (1200 bp) resulted from sequence reconstruction by Emirge amplicon. py using SortMeRNA-filtered reads with the longest sequence being 1495 bp. Similarly, MATAM reconstructed longer sequences having a median greater than 1125 bp using the three input datasets compared to Emirge. py (median < 1100 bp). The rRNA sequences of the eight less abundant microorganisms of the MOCK (abundance lesser than 0. 02%) were not reconstructed by any combination of tools and hence could not be detected.
Taxonomic affiliation of reconstructed rRNA sequences from the MOCK community Next, we evaluated the efficiency of seven taxonomic classifiers commonly used: Centrifuge, Kraken2, mothur, RDP classifier, SPINGO, BLCA and QIIME2's Sklearn classifier. The first input data was the 53 clustered sequences obtained previously from Emirge amplicon. py and MetaRib. Precision, recall and F-measure were then calculated by comparing the 'expected' classifications of the 20 detectable microorganisms of the mock community to the classifications predicted by each taxonomic classifier using the SILVA 138. 1 database (Figure  4A ). Classification using Centrifuge, Kraken2 and mothur were limited to the genera level. Mothur and Kraken2 performed better than Centrifuge with a precision of 84. 6% and 71. 4%, respectively, compared to 69.
Compilation of RiboTaxa After evaluating all the above tools, the best performing tools were compiled to build a pipeline enabling efficient and accurate taxonomic identification of any microbial community starting from raw shotgun metagenomics data. BBTools and FastQC were integrated to perform quality control. For SSU sequence reconstruction, Emirge amplicon. py (giving the best results as de-scribed previously) was chosen. For optimal results, Emirge amplicon. py was fed with SSU reads filtered by SortMeRNA. During tools comparison, MetaRib also performed outstandingly, identifying two microorganisms of the mock community, left out by Emirge amplicon. py. Thus, to optimise SSU sequence reconstruction, MetaRib was also included. Complementarily, MetaRib was fed with unfiltered reads as this tool included a dereplication step during iteration and could handle large datasets in low computational time.
To calculate the relative abundance, unfiltered high-quality reads were mapped onto the reconstructed sequences using BBMap (BBTools). Prior to taxonomic affiliation, the reconstructed rRNA gene sequences obtained from MetaRib and EMIRGE were clustered at 97% using VSEARCH and abundances from BBmap were summed up for sequences sharing the same cluster. Finally, Sklearn classifier (QIIME 2) was used to infer taxa. Throughout the study, the default database used for Ribo-Taxa was SILVA SSU 138. 1. RiboTaxa was then compared to existing taxonomic classifiers such as Centrifuge, Kraken2, METAXA2, phyloFlash, SPINGO, BLCA and MEGAN6 to evaluate their sensitivity, specificity and computation time during taxonomic classification using then MOCK and synthetic community.
Comparison of RiboTaxa with other existing tools on the MOCK dataset Here, only the 20 detectable microorganisms were looked for in the resulting taxonomic classification. Centrifuge, Kraken2, MEGAN6, METAXA2, SPINGO and BLCA inferred taxa directly to short Illumina reads and were labelled as 'Without sequence reconstruction'. RiboTaxa and phyloFlash included a step of SSU sequence assembly before microbial profiling and were labelled as 'With sequence reconstruction' (Table  2 ). All performance measures were computed at genus and species level (except for Centrifuge and Kraken2) to compare between all classifiers. Among classifiers 'without sequence reconstruction', METAXA2 performed the best at identifying the 20 genera correctly (recall 100%) while BLCA and Centrifuge, both, yielded the lowest recall (45%) of the microbial MOCK community with 11 and 9 correct taxa, respectively.