LotuS: an efficient and user-friendly OTU processing pipeline

Background 16S ribosomal DNA (rDNA) amplicon sequencing is frequently used to analyse the structure of bacterial communities from oceans to the human microbiota. However, computational power is still a major bottleneck in the analysis of continuously enlarging metagenomic data sets. Analysis is further complicated by the technical complexity of current bioinformatics tools. Results Here we present the less operational taxonomic units scripts (LotuS), a fast and user-friendly open-source tool to calculate denoised, chimera-checked, operational taxonomic units (OTUs). These are the basis to generate taxonomic abundance tables and phylogenetic trees from multiplexed, next-generation sequencing data (454, illumina MiSeq and HiSeq). LotuS is outstanding in its execution speed, as it can process 16S rDNA data up to two orders of magnitude faster than other existing pipelines. This is partly due to an included stand-alone fast simultaneous demultiplexer and quality filter C++ program, simple demultiplexer (sdm), which comes packaged with LotuS. Additionally, we sequenced two MiSeq runs with the intent to validate future pipelines by sequencing 40 technical replicates; these are made available in this work. Conclusion We show that LotuS analyses microbial 16S data with comparable or even better results than existing pipelines, requiring a fraction of the execution time and providing state-of-the-art denoising and phylogenetic reconstruction. LotuS is available through the following URL: http://psbweb05.psb.ugent.be/lotus.


Background
Next generation sequencing platforms are reducing the cost of collecting metagenomic data from large environmental and clinical microbial ecosystems. With 16S rDNA amplicon sequencing becoming a mainstream approach in these research areas, there is a need to optimize computer resources to handle this data.
Although online services exist to process 16S rDNA data such as the Ribosomal Database Project (RDP) pipeline [1] or the PyroTagger pipeline [2], a single HiSeq run can yield up to 6 × 10 9 sequences 1 , which challenges uploading capabilities. Large initiatives, like the Human Microbiome Project (HMP), have used the 16S pipelines Quantitative Insights Into Microbial Ecology (QIIME) [3] and mothur [4], applications that can be installed and run locally. mothur follows the philosophy of incorporating all tools in one software package, while QIIME (partly) relies on 3rd party software. Both pipelines present a complete package with tools to interpret the community composition, but more work-intensive components like denoising and sequence clustering are designed for cluster environments.
We developed less OTU scripts (LotuS) as an opensource pipeline in Perl and C++, which clusters operational taxonomic units (OTUs), generates taxonomic-level abundance matrices and a phylogenetic tree of the OTUs directly from non-demultiplexed sequencing files ( Figure 1). Currently, 454, HiSeq and MiSeq technologies are supported. Two popular proprietary programs, UPARSE [5] and RDP classifier [6], were incorporated to cluster OTUs, remove sequencing noise and chimeric sequences and classify OTUs taxonomically. Alternatively, OTU seed sequences can be classified by BLAST lowest common ancestor (LCA) comparison to either greengenes [7] or SILVA 16S rDNA databases [8]. LotuS installation does not require root access or changes in system paths and is automatically performed with a script that can additionally also update LotuS to newer versions. This script downloads, installs and configures all proprietary software and databases as well as configuring LotuS to these, with the exception of UPARSE, which requires user registration. Tutorials about LotuS usage and subsequent R numerical analysis are provided on the LotuS webpage (http://psbweb05.psb.ugent.be/ lotus/).

Sequence filtering
Sequences are filtered with a novel C++ program bundled with the pipeline, sdm, that has evolved over the course of this project from a simple demultiplexer to a general purpose sequence file formatting, quality filtering/adjustment and OTU sequence ("Seed") picking tool, optimized for speed and a high recovery rate of sequences. Further benefits of this program include the simultaneous sorting of input sequences into "high-" and "mid-" quality sequences, dependent on the overall quality and length. These two sequence bins will be important in the further LotuS workflow: the high quality sequences are used in the sensitive OTU clustering process, reducing the number of spurious OTU's caused by sequencing errors and not biological diversity. Both, high-and mid-quality sequences are used to count the occurrence of OTU's in single samples, with the intent that mid-quality sequences, which are mapping to an established OTU, will not confound overall diversity measures but add to the count of given biological entities present in a sample.
Default sdm options for 454 and MiSeq sequences are provided with LotuS; these can be modified to filter input sequences after average quality, accumulated error over the sequence, quality in a freely definable window and remove 5′ low-quality bases filtered for these criteria. Furthermore, sequences are filtered for min/max length, ambiguous nucleotides, max barcode and primer errors, polynucleotide runs, and trimmed for adapter sequences, if present. sdm/LotuS accepts fasta + quality, fastq and gzipped versions of these as input. Furthermore, sdm can be used on non-demultiplexed sequences to do a quality filtering of sequences, e.g. prior to assembly [9,10].

OTU clustering
Filtered sequences are clustered into OTUs with UPARSE. UPARSE is implemented as described by [5], with the exception that the "-cluster_OTUs" command is executed with the additional parameters -uparse_maxhot 62 -uparse_maxdrop 12 that increase the number of potential best hits that are explicitly aligned. This makes the overall clustering slightly slower, but in our experience, it results into more consistent OTUs on MiSeq datasets in case studies, thereby reducing the total number of OTUs. Also, we noticed that default UPARSE can propose a minor fraction of OTUs that are overlapping with other OTUs (mapped within the desired OTU sequence similarity, e.g. 97%, to existing OTUs). OTU abundance is estimated by mapping mid-and high-quality filtered input sequences onto the newly created OTUs.

Seed extension
The OTU sequence, here termed OTU seed sequence, should fulfil the following criteria to represent the OTU for sequence matching, taxonomic annotation and tree reconstruction: it should be as long as possible, be close to the median of all amplicons clustered to the OTU (representing the centre of the OTU) and contain the least amount of sequencing errors. One common practice is to use a consensus sequence of all amplicons clustered to an OTU, also used by UPARSE. However, a consensus sequence could be the average of two or more strains that constitute the OTU and the UPARSE denoising algorithm is critically dependent on sequences pruned to a reduced length (typically 250bp for 454 sequencing). To resolve this, sdm searches within all high quality input sequences for a sequence matching the above criteria, in a process we call "Seed extension". In brief, all input sequences are aligned with the consensus OTU sequence using usearch [11]. From these, sdm selects iteratively the hit being closest to the OTU median, having the highest overall mean accumulated error and the longest overall sequence length, with low-quality 5′ nucleotides being removed. In the case of paired MiSeq or HiSeq sequences, the highest-quality pair is selected and these are then merged using flash [12].

Taxonomic annotation of OTUs
By default, the OTU taxonomy is derived from RDP naïve bayes classifier annotations [6]. The minimum acceptance confidence for RDP is by default set to 0.8 but can be modified via "-rdp_thr". Alternatively, extended OTU seed sequences are aligned against a reference 16S rDNA database with BLAST + [13]; currently, we support greengenes [7] and SILVA [14] 16S rDNA databases, but this is easily extendable to other ribosomal databases. A lowest common ancestor algorithm is used to assign a taxonomy modified to only consider hits within 1.5% sequence identity to the overall best hit to the current OTU seed, with a limit of 200 total hits. This threshold is chosen to allow for sequencing errors to still include relevant hits, while limiting the space of possible hits within a close enough range to the current best hit; it was validated on 200 bp simulated reads (see Additional file 1) and is a good tradeoff between precision and specificity of taxonomic assignments (Additional file 2). At each taxonomic level, the taxonomic consistency of each hit is evaluated. By default, a taxonomic assignment is accepted if >90% of references are the same taxa, modifiable with the LotuS parameter "-LCA_frac". If a reference has no taxonomic information for a given level (and taxonomic assignments were consistent thus far), it is discarded. In an extreme example, this can lead to species-level taxonomic assignments even if only one reference is assigned to species level but all other references have no taxonomic information.
Furthermore, the sequence similarity of the best hit is used to delimit the taxonomic depth of an assignment. Even if the best hit has a known species name, if the identity of the OTU to the reference is 96%, the taxonomy will only be used up to genus level, as a best reference at a 96% threshold could indicate that this OTU is not represented by a reference species in the database. The default parameters are to limit species at 97%, genus at 95%, family at 93%, class at 91%, order at 88% and phylum at 78% sequence similarity, though users can change these parameters as the established 16S rRNA gene taxonomy does strictly speaking not follow consistent cutoffs for different taxonomic levels [15,16].

Multiple alignment and phylogenetic reconstruction
The calculation of a phylogenetic tree of the extended OTU seed sequences is an optional step. The phylogenetic tree can be used to fulfil requirements for calculating diversity indices such as Faith's phylogenetic diversity [17] or between sample UniFrac distances [18]. For this step, the sequences are aligned with Clustal Ω [19] with default parameters for nucleotide alignments. From the aligned sequences, a phylogenetic tree is reconstructed using the gamma model of sequence evolution (options "-nt -gamma -no2nd -fastest -spr 4") in FastTree2 [20], as recommended by its author (http://www.microbesonline.org/fasttree/) and is saved in Phylip format.

LotuS Output
LotuS saves the output in the specified output folder in a structure that can be directly integrated into specialized analytical packages of statistical software packages. OTU abundance matrix, the phylogeny for each OTU, a phylogenetic tree and a .biom formatted OTU matrix are stored in this folder. Three subfolders contain a) the run logs and processing reports, b) copies of configuration files and c) higher level taxa abundance matrices.

Results and discussion
We used a simulation of 1000 greengenes [7] 16S sequences, truncated and randomly mutated (see Additional file 1), to validate the LCA algorithm and research the influence of read length on taxonomic classification. This showed that longer reads are assigned with a higher confidence in RDP (Additional file 3a) and the fraction of 16S reads that remains unclassified using our LCA algorithm significantly decreases with longer read length, when using a reference database from which the 1,000 simulated reads were removed (Additional file 3b). With reads ≥250 bp, our simulation converged to 100% precision and specificity (Additional file 4a,b). To simulate a rare situation where no close relatives are present in the reference database, we used the LCA algorithm only on database hits that had <97% identity to the target read. Here the fraction of taxonomic assignments is in general lower and no species-level assignments were made, as expected given default parameters (Additional file 3c). Precision and specificity of assignments were lowered, and here, the longer read lengths were especially important (Additional file 4c,d).
Read quality is decreasing with increasing read length in 454 and illumina sequencing [5], and UPARSE improves OTU clustering by using only the high-quality 5′ DNA. The here proposed OTU seed extension is important for taxonomic classification and tree building, especially for badly characterized species. It takes advantage of the improved, fast OTU clustering, while using long, high-quality reads for taxonomic annotations and multiple sequence alignments.
We tested the validity and performance of LotuS on cecal gut samples from five different mice strains [21]. It consists of two 454 GS FLX runs and a total of 393,070 reads, the expected read length is 400-500 nucleotides. We used five methods (described in Additional file 1) to calculate OTU abundance matrices: LotuS using RDP taxonomy (LR), LotuS using BLAST taxonomy (LB), QIIME de novo OTU creation (QDN), QIIME with reference-based OTUs (QRE) and mothur (MOT   clippings and filtering is used. In mothur, this was 248 bp, because only high-quality and informative sequence parts were retained, similar to the 250 bp cutoff used by LotuS to cluster OTUs. OTU and taxon-read abundance matrices from these three pipelines were processed in the R statistical computer language. In general, the sample composition is very consistent between the five methods ( Figure 2), when taking into consideration differences in taxonomy due to the database/assignment algorithm, that were different in all three pipelines. To measure how well the sample relationships are maintained if different pipeline versions are used, the correlations between Bray-Curtis (BC) and weighted UniFrac (wUF) inter-sample distances were compared in an approach similar to the Mantel test [22] (Additional file 5: Table S1 and Additional file 6: Table S2). The average correlation was 0.993 ± 0.005, 0.988 ± 0.001 and 0.958 ± 0.02, on OTU BC, wUF and genus BC distances, respectively.. Both LotuS-derived wUF were closer to those of mothur than QRE or QDN; they were also closer to QRE than QDN, where the guide tree is based on the greengenes phylogenetic tree, indicating that our de novo constructed guide tree is similar to the greengenes one. OTU diversity varied strongly between methods [5]: QDN and QRE clustered most OTUs (6,148 and 1,467), followed by mothur (913) and LotuS (369) ( Table 1). When rarefying samples to 2,000 reads, the average OTU number across all samples was similar between LotuS and mothur, with 119 and 118 OTUs/sample. QIIME-derived data had higher sample diversity, with 243 and 201 rarefied OTUs in QDN and QRE, respectively. The total number of reads in the abundance matrices differs only slightly between pipelines (345,801 ± 9,622).
The number of unassigned taxa at a given taxonomic level is usually lowest in QRE and highest in QDN. QRE does not use a LCA algorithm to exclude cross hits and this is obvious at genus level where only 54.7% of taxa are undefined, which is lower than LB (65.9% of genera unassigned, Additional file 2). QRE is by default excluding sequences that do not match a reference with <97% and the unassigned taxa are thus from the greengenes taxonomic assignment. From the known genera, 55% are shared between all methods, and this applies for 71% of families (Additional file 7). Execution times are the most outstanding difference between tested pipelines ( Table 2). The same dataset was processed with all pipelines using the same 2.7 GHz Intel Core i7 MacBook (1 core). Default LotuS was 25 to 200 times faster than the other pipelines and demultiplexing and quality filtering with sdm was 4-6 times faster than these processes in other pipelines. However, when matching reads to a reference database, LotuS is limited by the speed of BLAST+.
One of the main problems in testing the validity of metagenomic processing is that the real taxonomic composition of a metagenomic sample is unknown or only available for artificial datasets [23], and thus, the absolute error due to bioinformatics processing is hard to estimate. Here we sequenced 40 technical replicates in two separate MiSeq runs to circumvent this problem, as we can determine a relative error between replicates that could be attributed to a pipeline. In this dataset, LotuS recovered 28,789,221 reads from the two MiSeq runs, 19% or 10% more reads than QRE or QDN, respectively. The average sample richness at 15,000 rarefied sequences was highest in QDN (1,333 ± 393.5) followed by LotuS (720.5 ± 186) and closed-reference OTU picking QRE (469.6 ± 146). This could be related to QRE being too stringent as not all OTUs might be present in reference database and QDN being too lenient. Comparing the reproducibility of technical replicates between the two MiSeq runs, we used Canberra, Bray-Curtis and Jensen-Shannon distances between technical replicates. The mean of these was measured for each pipeline variant and on OTU, genus and family level. Overall, no consistent trend was detected for one pipeline recovering more reproducible compositions that any other  ( Table 3). Comparing the reproducibility of richness, LotuS OTU matrices (identical between LR and LB) had an average richness difference of 5.7 ± 5.9%, lower compared to QRE (7.5 ± 6%, p = 0.0065 compared to LotuS) and QDN (6.4 ± 6.3%, p = 0.17 compared to LotuS), thus richness was more stable with LotuS among technical replicates. Runtime was again fastest in LotuS (Table 2), though differences in execution time were not as extreme as observed for the 454 dataset with LR being~8 times faster than Qiime.

Conclusions
The novel LotuS pipeline is able to handle small to very large 16S datasets on a personal computer and effortlessly integrate multiple sequencing runs. Computational efficiency is very high due to a selection of state-of-theart proprietary software like UPARSE for denoising and sequence clustering and sdm for demultiplexing and sequence filtering. Comparison to other pipelines suggests a high similarity in higher taxonomic composition to existing tools, but on OTU level, the de-novo-called OTU shows an increased richness compared to closed-reference OTU calling and less richness than non-denoised de novo OTU calling, as expected. This pipeline has the advantage of state-of-the-art, flowgram-independent denoising and long, high-quality OTU sequences from the OTU seed extension step used for phylogenetic tree construction and taxonomic annotation.