VSEARCH: a versatile open source tool for metagenomics — TEI-tagged view

Torbjørn RognespersName, Tomáš FlouripersName, Ben NicholspersName, Christopher QuincepersName, Frédéric MahépersName (2016date). “VSEARCH: a versatile open source tool for metagenomics.” PeerJtitle·j.

Object of study — open-source, memory-unrestricted nucleotide sequence processing for metagenomicsterm, genomics, and population genomics

Research question — whether a transparent 64-bit implementation can reproduce and extend USEARCHterm functionality with competitive accuracy, speed, and robustness

Contribution — VSEARCHterm, a multithreaded sequence-analysis suite with heuristic candidate search, optimal alignment, clustering, dereplication, chimeraterm detection, FASTQterm processing, paired-read merging, and subsampling

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

VSEARCHterm: a versatile open source tool for metagenomicsterm Abstract VSEARCHterm is an open source and free of charge multithreaded 64-bit toolterm for processing and preparing metagenomicsterm, genomics and population genomics nucleotide sequence dataterm. It is designed as an alternative to the widely used USEARCHterm tool ( Edgar, 2010 ) for which the source code is not publicly available, algorithm details are only rudimentarily described, and only a memory-confined 32-bit version is freely available for academic use.

When searching nucleotide sequences, VSEARCHterm uses a fast heuristic based on words shared by the query and target sequencesterm in order to quickly identify similar sequences, a similar strategy is probably used in USEARCHterm. VSEARCHterm then performs optimal global sequence alignmentterm of the query against potential target sequences, using full dynamic programmingterm instead of the seed-and-extend heuristicterm used by USEARCHterm. Pairwise alignments are computed in parallel using vectorisationterm and multiple threadsterm.

VSEARCHterm includes most commands for analysing nucleotide sequences available in USEARCHterm version 7 and several of those available in USEARCHterm version 8, including searching (exact or based on global alignment), clustering by similarity (using length pre-sorting, abundance pre-sorting or a user-defined order), chimeraterm detection (reference-based or de novo ), dereplication (full length or prefix), pairwise alignment, reverse complementation, sorting, and subsampling. VSEARCHterm also includes commands for FASTQterm file processing, i.e., format detection, filtering, read quality statistics, and merging of paired reads. Furthermore, VSEARCHterm extends functionality with several new commands and improvements, including shuffling, rereplicationterm, masking of low-complexity sequences with the well-known DUST algorithmterm, a choice among different similarity definitions, and FASTQterm file format conversion. VSEARCHterm is here shown to be more accurate than USEARCHterm when performing searching, clustering, chimeraterm detection and subsampling, while on a par with USEARCHterm for paired-ends read merging. VSEARCHterm is slower than USEARCHterm when performing clustering and chimeraterm detection, but significantly faster when performing paired-end readsterm merging and dereplication. VSEARCHterm is available at https://github.com/torognes/vsearchterm under either the BSD 2-clause license or the GNU General Public License version 3.0.

VSEARCHterm has been shown to be a fast, accurate and full-fledged alternative to USEARCHterm. A free and open-source versatile tool for sequence analysis is now available to the metagenomicsterm community.

Introduction Rockström et al. (2009) and Steffen et al. (2015) presented biodiversity loss as a major threat for the short-term survival of humanity. Recent progress in sequencing technologies have made possible large scale studies of environmental genetic diversity, from deep sea hydrothermal vents to Antarctic lakes ( Karsenti et al., 2011 ), and from tropical forests to Siberian steppes ( Gilbert, Jansson & Knight, 2014 ). Recent clinical studies have shown the importance of the microbiomes of our bodies and daily environments for human health ( Human Microbiome Project Consortium, 2012 ). Usually focusing on universal markers (e.g., 16S rRNA, ITS, COI), these targeted metagenomicsterm studies produce many millions of sequences, and require open-source, fast and memoryterm efficient tools to facilitate their ecological interpretation.

Several pipelines have been developed for microbiome analysis, among which mothur ( Schloss et al., 2009 ), QIIME ( Caporaso et al., 2010 ), and UPARSE ( Edgar, 2013 ) are the most popular. QIIME and UPARSE are both based on USEARCHterm ( Edgar, 2010 ), a set of tools designed and implemented by Robert C. Edgar, and available at http://drive5.com/usearchterm/ . USEARCHterm offers a great number of commands and options to manipulate and analyse FASTQterm and FASTA files. However, the source code of USEARCHterm is not publicly available, algorithm details are only rudimentarily described, and only a memory-confined 32-bit version is freely available for academic use.

We believe that the existence of open-source solutions is beneficial for end-users and can invigorate research activities. For this reason, we have undertaken to offer a high quality open-source alternative to USEARCHterm, freely available to users without any memoryterm limitation. VSEARCHterm includes most of the USEARCHterm functions in common use, and further development may add additional features. Here we describe the details of the VSEARCHterm implementation. To assess its performance in terms of speed and quality of results, we have evaluated some of the most important functions (searching, clustering, chimeraterm detection and subsampling) and compared them to USEARCHterm. We find that VSEARCHterm delivers results that are better or on a par with USEARCHterm results.

Materials and Methods Algorithms and implementation Below is a brief description of the most important functions of VSEARCHterm and details of their implementation. VSEARCHterm command line options are shown in italics, and should be preceded by a single (-) or double dash (-TI:\,-) when used.

Reading FASTA and FASTQterm files Most VSEARCHterm commands read files in FASTA or FASTQterm format. The parser for FASTQterm files in VSEARCHterm is compliant with the standard as described by Cock et al. (2010) and correctly parses all their tests files. FASTA and FASTQterm files are automatically detected and many commands accept both as input. Files compressed with gzip or bzip2 are automatically detected and decompressed using the zlib library by Gailly & Adler (2016) or the bzip2 library by Seward (2016) , respectively. Data may also be piped into or out of VSEARCHterm, allowing for instance many separate FASTA files to be piped into VSEARCHterm for simultaneous dereplication, or allowing the creation of complex pipelines without ever having to write on slow disks.

VSEARCHterm is a 64-bit program and allows very large datasets to be processed, essentially limited only by the amount of memoryterm available. The free USEARCHterm versions are 32-bit programs that limit the available memoryterm to somewhere less than 4GB, often seriously hampering the analysis of realistic datasets.

Writing result files VSEARCHterm can output results in a variety of formats (FASTA, FASTQterm, tables, alignments, SAM) depending on the input format and command used. When outputting FASTA files, the line width may be specified using the fasta_width option, where 0 means that line wrapping should be turned off. Similar controls are offered for pairwise or multiple sequence alignments.

Searching Global pairwise sequence comparisonterm is a core functionalityterm of VSEARCHterm. Several commands compare a query sequence against a database of sequences: all-vs-all alignment ( allpairs_global ), clustering ( cluster_fast , cluster_size , cluster_smallmem ), chimeraterm detection ( uchime_denovo and uchime_ref ) and searching ( usearchterm_global ). This comparison function proceeds in two phases: an initial heuristic filtering based on shared wordsterm, followed by optimal alignment of the query with the most promising candidates.

The first phase is presumably quite similar to USEARCHterm ( Edgar, 2010 ). Heuristics are used to identify a small set of database sequencesterm that have many wordsterm in common with the query sequence. Wordsterm (or k -mers) consist of a certain number k of consecutive nucleotides of a sequence (8 by default, adjustable with the wordlengthterm·var option). All overlapping wordsterm are included. A sequence of length n then contains at most n − k + 1 unique wordsterm. VSEARCHterm counts the number of shared wordsterm between the query and each database sequenceterm. Wordsterm that appear multiple times are counted only once. To count the wordsterm in the database sequencesterm quickly, VSEARCHterm creates an index of all the 4 k possible distinct wordsterm and stores information about which database sequencesterm they appear in. For extremely frequent wordsterm, the set of database sequencesterm is represented by a bitmapterm; otherwise the set is stored as a list. A finer control of k -mer indexing is described for USEARCHterm by the pattern (binary string indicating which positions must match) and slots options. USEARCHterm has such options but seems to ignore them. Currently, VSEARCHterm ignores these two options too. The minimum number of shared wordsterm required may be specified with the minwordmatches option (10 by default), but a lower value is automatically used for short or simple query sequences with less than 10 unique wordsterm.

Comparing sequences based on statistics of shared wordsterm is a common method to quickly assess the similarity between two sequences without aligning them, which is often time-consuming. The D 2 statistic and related metrics for alignment-free sequence comparison have often been used for rapid and approximate sequence matching and their statistical properties have been well studied ( Song et al., 2014 ). The approach used here has similarities to the D 2 statistic, but multiple matches of the same wordterm are ignored.

In the second phase, searching proceeds by considering the database sequencesterm in a specific order, starting with the sequence having the largest number of wordsterm in common with the query, and proceeding with a decreasing number of shared wordsterm. If two database sequencesterm have the same number of wordsterm in common with the query, the shortest sequence is considered first. The query sequence is compared with each database sequenceterm by computing the optimal global alignmentterm. The alignment is performed using a multi-threaded and vectorised full dynamic programmingterm algorithm ( Needleman & Wunsch, 1970 ) adapted from SWIPE ( RognespersName, 2011 ). Due to the extreme memoryterm requirements of this method when aligning two long sequences, an alternative algorithmterm described by Hirschberg (1975) and Myers & Miller (1988) is used when the product of the length of the sequences is greater than 25,000,000, corresponding to aligning two 5,000 bp sequences. This alternative algorithmterm uses only a linear amount of memoryterm but is considerably slower. This second phase is probably where USEARCHterm and VSEARCHterm differ the most, as USEARCHterm by default presumably performs heuristic seed-and-extend alignment similar to BLAST ( Altschul et al., 1997 ), and only performs optimal alignment when the option fulldp (full dynamic programmingterm) is used. Computing the optimal pairwise alignment in each case gives more accurate results but is also computationally more demanding. The efficient and vectorised full dynamic programmingterm implementation in VSEARCHterm compensates that extra cost, at least for sequences that are not too long.

If the resulting alignment indicates a similarity equal to or greater than the value specified with the id optionterm, the database sequenceterm is accepted. If the similarity is too low, it is rejected. Several other options may also be used to determine how similarity is computed ( iddef , as USEARCHterm used to offer up to version 6), and which sequences should be accepted and rejected, either before (e.g., self , minqsize ) or after alignment (e.g., maxgaps , maxsubs ). The search is terminated when either a certain number of sequences have been accepted (1 by default, adjustable with the maxaccepts option), or a certain number of sequences have been rejected (32 by default, adjustable with the maxrejects option). The accepted sequences are sorted by sequence similarityterm·var and presented as the search results.

VSEARCHterm also includes a search_exact command that only identifies exact matches to the query. It uses a hash tableterm in a way similar to the full-length dereplicationterm command described below.

Clustering VSEARCHterm includes commands to perform de novo clustering using a greedy and heuristic centroid-based algorithmterm with an adjustable sequence similarity thresholdterm specified with the id option (e.g., 0.97). The input sequences are either processed in the user supplied order ( cluster_smallmem ) or pre-sorted based on length ( cluster_fast ) or abundance (the new cluster_size option). Each input sequence is then used as a query in a search against an initially empty database of centroid sequences. The query sequence is clustered with the first centroid sequence found with similarity equal to or above the threshold. The search is performed using the heuristic approach described above which generally finds the most similar sequences first. If no matches are found, the query sequence becomes the centroid of a new cluster and is added to the database. If maxaccepts is higher than 1, several centroids with sufficient sequence similarityterm·var may be found and considered. By default, the query is clustered with the centroid presenting the highest sequence similarityterm·var (distance-based greedy clusteringterm, DGC), or, if the sizeorder option is turned on, the centroid with the highest abundance (abundance-based greedy clusteringterm, AGC) ( He et al., 2015 ; Westcott & Schloss, 2015 ; Schloss, 2016 ). VSEARCHterm performs multi-threaded clustering by searching the database of centroid sequences with several query sequences in parallel. If there are any non-matching query sequences giving rise to new centroids, the required internal comparisons between the query sequences are subsequently performed to achieve correct results. For each cluster, VSEARCHterm can create a simple multiple sequence alignment using the center star method ( Gusfield, 1993 ) with the centroid as the center sequence, and then compute a consensus sequence and a sequence profile.

Dereplication and rereplicationterm Full-length dereplicationterm ( derep_fulllength ) is performed using a hash tableterm with an open addressing and linear probing strategy based on the Google CityHash hash functions (written by Geoff Pike and Jyrki Alakuijala, and available at https://github.com/google/cityhash ). The hash tableterm is initially empty. For each input sequence, the hash is computed and a lookup in the hash tableterm is performed. If an identical sequence is found, the input sequence is clustered with the matching sequence; otherwise the input sequence is inserted into the hash tableterm.

Prefix dereplicationterm ( derep_prefix ) is also implemented. As with full-length dereplicationterm, identical sequences are clustered. In addition, sequences that are identical to prefixes of other sequences will also be clustered together. If a sequence is identical to the prefix of multiple sequences, it is generally not defined how prefix clustering should behave. VSEARCHterm resolves this ambiguity by clustering the sequence with the shortest of the candidate sequences. If they are equally long, priority will be given to the most abundant, the one with the lexicographically smaller identifier or the one with the earliest original position, in that order.

To perform prefix dereplicationterm, VSEARCHterm first creates an initially empty hash tableterm. It then sorts the input sequences by length and identifies the length s of the shortest sequence in the dataset. Each input sequence is then processed as follows, starting with the shortest: If an exact match to the full input sequence is found in the hash tableterm, the input sequence is clustered with the matching hash tableterm sequence. If no match to the full input sequence is found, the prefixes of the input sequence are considered, starting with the longest prefix and proceeding with shorter prefixes in order, down to prefixes of length s . If a match is now found in the hash tableterm, the sequences are clustered, the matching sequence is deleted from the hash tableterm and the full input sequence is inserted into the hash tableterm instead. If no match is found for any prefix, the full sequence is inserted into the hash tableterm. In the end, the remaining sequences in the hash tableterm will be output with accumulated abundances for all sequences in each cluster.

In order to identify matches in the hash tableterm during prefix dereplicationterm, a hash is computed for each full-length input sequence and all its prefixes. The hash function used is the 64-bit Fowler–Noll–Vo 1a hash functionterm ( Fowler, Noll & Vo, 1991 ), which is simple and quick to compute for such a series of sequences by adding one nucleotide at a time.

The sequences resulting from dereplication and many other commands may be relabeled with a given prefix followed by a sequentially increasing number. VSEARCHterm exclusively also offers the possibility of relabelling each sequence with the SHA-1term ( Eastlake & Jones, 2001 ) or MD5term ( Rivest, 1992 ) message digest (hash) of the sequence. These are strings that are highly likely to be unique for each sequence. Before the digest is computed, the sequence is normalized by converting U’s to T’s and converting all symbols to upper case. VSEARCHterm includes public domain code for the MD5term algorithm written by Alexander Peslyak, and for SHA1term by Steve Reid and others.

VSEARCHterm also includes a new command ( rereplicate ) to perform rereplicationterm that can be used to recreate datasets as they were before full-length dereplicationterm, but of course original labels cannot be recreated.

Chimeraterm detection Chimerasterm are detected either de novo ( uchime_denovo command) or with a reference database ( uchime_ref command) using the UCHIME algorithmterm described by Edgar et al. (2011) . VSEARCHterm will divide each query sequence into four segments and look for similarity of each segment to sequences in the set of potential parents using the heuristic search function described earlier. It will consider the four best candidates for each segment using maxaccepts 4 and maxrejects 16, and an id threshold of 0.55. VSEARCHterm optionally outputs borderline sequences, that is, sequences having a high enough score (as specified with the minh option) but with too small a divergence from the closest parent (as specified with the mindiv option). Multi-threading is supported for reference-based chimeraterm detection.

Low-complexity sequence masking VSEARCHterm includes a highly optimized and parallelized implementation of the Dust algorithmterm by Tatusov and Lipman for masking of simple repeats and low-complexity nucleotide sequencesterm. It is considerably faster than the implementation of the same algorithm in USEARCHterm. Their code available at ftp://ftp.ncbi.nlm.nih.gov/pub/tatusov/dust/version1/src/ is in the public domain. VSEARCHterm uses this algorithm by default, while USEARCHterm by default uses an undocumented rapid masking algorithm called fastnucleo . VSEARCHterm performs soft-maskingterm automatically for the pairwise alignment, search, clustering and chimeraterm detection commands. This behaviour can be controlled with the hardmask option to replace masked symbols with N’s instead of lower-casing them, and the dbmask and qmask options, which selects the masking algorithm (none, dust or soft) used for the database and query sequences, respectively. Masking may also be performed explicitly on an input file using the fastx_mask and maskfasta commands.

FASTQterm file processing VSEARCHterm includes commands to detect the FASTQterm file version and the range of quality scores used ( fastqterm_chars ), as well as two commands for computing sequence quality statistics ( fastqterm_stats and fastqterm_eestats ). It can also truncate and filter sequences in FASTQterm files based on various criteria ( fastqterm_filter ). A new command is added to convert between different FASTQterm file versions and quality encodings ( fastqterm_convert ), e.g., from the old Phred+64 encoded Illumina FASTQterm files to the newer Phred+33 format.

Merging of paired-end readsterm Merging of paired-end readsterm is supported by VSEARCHterm using the fastqterm_mergepairs command. The method used has some similarity to PEAR ( Zhang et al., 2014 ) and recognises options similar to USEARCHterm. The algorithm computes the optimal ungapped alignmentterm of the overlapping region of the forward sequence and the reverse-complemented reverse sequence. The alignment requires a minimum overlap length (specified with the fastqterm_minovlen option, default 10), a maximum number of mismatches ( fastqterm_maxdiffs option, default 5), and a minimum and maximum length of the merged sequence ( fastqterm_minmergelen option, default 1, and fastqterm_maxmergelen option, default infinite). Staggered read pairs, i.e., read pairs where the 3′ end of the reverse read has an overhang to the left of the 5′ end of the forward read, are not allowed by default, but may be turned on by the fastqterm_allowmergestagger option. VSEARCHterm uses a match score (alpha) of +4 and a mismatch score (beta) of −5 for perfect quality residues. These scores are weighted by the probability that these two residues really match or mismatch, respectively, taking quality scores into account. These probabilities are computed in a way similar to PEAR score method 2 described in ‘Algorithms and implementation’ of the PEAR paper ( Zhang et al., 2014 ), but VSEARCHterm assumes all nucleotide background frequencies are 0.25. When merging sequences, VSEARCHterm computes posterior quality scoresterm for the overlapping regions as described by Edgar & Flyvbjerg (2015) . For speed, scores and probabilities are pre-computed for all possible quality scores.

Sorting and shuffling VSEARCHterm can sort FASTA files by decreasing sequence lengthterm·var ( sortbylength ) or abundance ( sortbysize ). VSEARCHterm can also perform shuffling of FASTA files in random order ( shuffle ). A seed value for the pseudo random number generatorterm may be provided by the randseed option to obtain replicable results.

Subsampling Sequences in FASTA and FASTQterm files can be subsampledterm ( fastx_subsample ) by randomly extracting a certain number ( sample_size ) or percentage ( sample_pct ) of the input sequences. Abundances may be taken into account, giving results as if the input sequences were rereplicated, subsampledterm and then dereplicated.

Results and Discussion Supported commands and options VSEARCHterm implements the following commands available in USEARCHterm version 7: allpairs_global , cluster_fast , cluster_smallmem , derep_fulllength , derep_prefix , fastqterm_chars , fastqterm_filter , fastqterm_mergepairs , fastqterm_stats , fastx_mask , maskfasta , sortbylength , sortbysize , uchime_denovo , uchime_ref and usearchterm_global . In addition, the following commands available in USEARCHterm version 8 have been implemented: fastqterm_eestats , fastx_revcomp , fastx_subsample and search_exact . VSEARCHterm additionally includes a few new commands that do not exist in USEARCHterm: cluster_size , fastqterm_convert , rereplicate and shuffle .

Some USEARCHterm version 7 commands have not yet been implemented in VSEARCHterm. We have not prioritized commands related to amino acid sequences ( findorfs ), local alignment ( allpairs_local , pairs_local , search_local , ublast ), brute-force search ( search_global , pairs_global ), UDB databases ( makeudb_ublast , makeudb_usearchterm , udb2fasta , udbinfo , udbstats ), and the UPARSE pipeline ( cluster_otus , uparse_ref ).

Almost all USEARCHterm 7 options are supported, except for those related to non-standard database indexing ( alpha , dbaccelpct , dbstep , pattern , slots ) as well as local alignments and alignment heuristics ( band , hspw , lext , lopen , matrix , minhsp , xdrop_g , xdrop_nw , xdrop_u ).

The same command and option names as in USEARCHterm version 7 have generally been used in order to make VSEARCHterm an almost drop-in replacementterm. In fact, in QIIME many commands will run fine if an alias or link from usearchterm to vsearchterm is made. Detailed documentation of VSEARCHterm is available as a man page. We will consider adding further commands and options to VSEARCHterm in the future.

Performance assessment The performance of the most important functions of VSEARCHterm version 2.0.3 (64-bit) was evaluated and compared to the 32-bit versions of USEARCHterm version 7.0.1090 and 8.1.1861. All datasets used were small enough to fit comfortably in the memoryterm allocated to a 32-bit process. Chimeraterm detection was also compared to UCHIME version 4.2. All tests were run on GNU/Linux CentOS 6.7 compute nodes with 16 physical cores (Intel(R) Xeon(R) CPU E5-2670 0 @ 2.60 GHz) and 64GB RAM. Programs were run with 8 threads, if possible. All times indicated are wall-clock times. All scripts and data necessary to perform the evaluations are available in the GitHub repository at https://github.com/torognes/vsearch-eval/ to enable independent replication.

Searching Evaluation of search accuracy was carried out as described in the USEARCHterm paper ( Edgar, 2010 ), its supplementary, and on the website ( http://drive5.com/usearchterm/benchmark_rfamterm.html ), by assessing the ability of the programs to identify RNA sequences belonging to the same family in RFAMterm ( Burge et al., 2013 ). The 383,004 sequences in Rfamterm version 11 were randomly shuffled and then the first sequence from each of the 2,085 (out of 2,208) families that contained at least 2 members was selected as a representative and used as a query against the remaining 380,919 sequences. The programs were run with options id 0.0, minseqlength 1, maxaccepts 1, maxrejects 32, and strand plus. If the matching sequence found belonged to the same family, it was considered a true positive, otherwise it was considered as a false positive. We combined the results from 20 shufflings and plotted the results in the ROC-like curve shown in Fig. 1 . For a false discovery rateterm between 0.010 and 0.015, VSEARCHterm is more accurate than USEARCHterm’s latest version. For lower values, the three programs have similar accuracies. At higher false discovery ratesterm, USEARCHterm version 8 has an advantage.

Search accuracy on the RFAMterm v11 dataset. USEARCHterm version 7 (blue), USEARCHterm version 8 (orange) and VSEARCHterm (black) was run using the usearchterm_global command on subsets of the RFAMterm dataset to identify members of the same families. The plot shows the true positive rate (also known as the recallterm or sensitivity) as a function of the false discovery rateterm at varying sequence similarityterm·var levels. This curve is based on data from 20 shufflings of the dataset.

The time to search the Rfamterm database as described above was measured. To avoid extremely short running times, 1,000 replicates of the datasets were used. USEARCHterm version 7 required on average 5 min 29 s for the search, USEARCHterm version 8 took 5 min 57 s, while VSEARCHterm took 5 min 26 s.

Clustering Westcott & Schloss (2015) have already carried out an evaluation of the clustering performance of VSEARCHterm. They tested the ability of several tools to assign OTUs for 16S rRNA sequences and “demonstrated that for the greedy algorithms VSEARCHterm produced assignments that were comparable to those produced by USEARCHterm making VSEARCHterm a viable free and open source alternative to USEARCHterm”. Schloss (2016) also evaluated de novo clustering by VSEARCHterm.

We independently evaluated the clustering accuracy of USEARCHterm and VSEARCHterm as described for Swarm ( MahépersName et al., 2014 ) using two mock datasets, one with an even and one with uneven composition of 57 archaea and bacteria. The datasets were first dereplicated. Then the taxonomy of the unique sequences was assigned by a search against the set of rRNA reference sequences representing the species in the mock datasets, carried out with the usearchterm_global command of USEARCHterm. The sequences were shuffled randomly 10 times and clustering was performed at 20 different similarity levels ranging from 80% to 99% in steps of 1%. Clustering was carried out in two ways, first using the cluster_fast command that pre-sorts the sequences by length, and then using the cluster_smallmem command that simply processes the sequences in the user-supplied order (in that case, we supplied sequences ordered by decreasing abundance using the sortbysize command). We then compared the clusters obtained to the assigned species and computed the recallterm, precisionterm and the adjusted Rand indexterm of the classifications. Recallterm measures to what extent amplicons assigned to the same speciesterm are grouped together in the same OTUterm (i.e., not over-splitting). Precisionterm measures to what extent amplicons in an OTU are assigned to the same speciesterm (i.e., not over-grouping). The adjusted Rand indexterm ( Rand, 1971 ; Hubert & Arabie, 1985 ) summarises both precisionterm and recallterm using a measure of agreement between two clusterings while adjusting for expected valuesterm by chance. The average values over the all shufflings are presented in Figs. 2 and 3 for the even and uneven datasets, respectively. For abundance-sorted sequences, the difference between VSEARCHterm and USEARCHterm version 8 is negligible. The difference is larger for length-sorted sequences. When using length sorting, USEARCHterm 8 (as well as version 7 on the even dataset) shows better precisionterm than VSEARCHterm for similarity levels below 93%. However, OTU delineation is usually performed at higher similarity values, typically 97%. In the case of our benchmark, overall accuracy as measured by the adjusted Rand indexterm is maximised at 95–97% similarity. This is precisely the region where for length sorting at least VSEARCHterm outperforms USEARCHterm.

Clustering accuracy on the even dataset. USEARCHterm version 7 (blue) and 8 (orange) and VSEARCHterm (black) was run using abundance sorting ( cluster_smallmem ) (A, C, E) and length sorting ( cluster_fast ) (B, D, F) on the even dataset. The performance is indicated with the adjusted Rand indexterm (A, B), recallterm (C, D) and precisionterm (E, F) metrics.

Clustering accuracy on the uneven dataset. USEARCHterm version 7 (blue) and 8 (orange) and VSEARCHterm (black) was run using abundance sorting ( cluster_smallmem ) (A, C, E) and length sorting ( cluster_fast ) (B, D, F) on the uneven dataset. The performance is indicated with the adjusted Rand indexterm (A, B), recallterm (C, D) and precisionterm (E, F) metrics.

The time usedterm for clustering is shown in Fig. 4 . The time usedterm depended on the dataset, algorithm and clustering thresholdterm. The USEARCHterm programs were in general 2–3 times faster than VSEARCHterm. In general the difference in speed was smaller for higher thresholds, especially at 99% similarity. Clustering with USEARCHterm using the fulldp option is an order of magnitude slower than VSEARCHterm, but contrary to our expectations it does not seem to improve accuracy.

Clustering speed. Median wall time in seconds to cluster the even (A, B) and uneven (C, D) datasets using USEARCHterm version 7 (blue) and 8 (orange) and VSEARCHterm (black) using abundance sorting ( cluster_smallmem ) (A, C) and length sorting ( cluster_fast ) (B, D).

Dereplication Measurements of dereplication speed were performed on the even and uneven datasets described earlier as well as on the BioMarKs dataset ( Logares et al., 2014 ). For full-length dereplicationterm ( derep_fulllength ) VSEARCHterm was about 40–50% faster than USEARCHterm version 7 and 50–70% faster than version 8 on all three datasets. All programs were approximately equally fast on prefix dereplicationterm ( derep_prefix ) of the even and uneven datasets. However, prefix dereplicationterm of the BioMarKs dataset was extremely slow with USEARCHterm. USEARCHterm version 7 used more than 4 min and version 8 more than 27 min, while VSEARCHterm used less than 4 s. The prefix dereplicationterm algorithm used in USEARCHterm appears ineffective when dealing with short sequences. Removing the 811 sequences shorter than 200 bp out of the 312,503 sequences of the BioMarKs dataset reduces the running time of USEARCHterm version 7 and 8 down to just 5 and 6 s, respectively.

Chimeraterm detection We evaluated the chimeraterm detection accuracy of VSEARCHterm and USEARCHterm in two ways, first using a method similar to that performed for UCHIME, and then using a new chimeraterm simulation procedure based on sequences from Greengenesterm ( DeSantis et al., 2006 ) and SILVAterm ( Quast et al., 2013 ) sequences.

Chimeraterm detection performance with the SIMM datasetterm. UCHIME (UC), USEARCHterm version 7 (U7) and 8 (U8), and VSEARCHterm (V) was run using the uchime_ref algorithm on the SIMM datasetterm that was originally also used to evaluate the UCHIME algorithmterm. Divergence is the percentage of similarity to the original sequences. Noise is either zero (–) or the percentage of indels (i1–i5) or substitutions (m1-5) added. The number of chimerasterm detected out of 100 of each type is shown. The best results in each category are shaded.

First we repeated the evaluation of the uchime_ref command described in the UCHIME paper ( Edgar et al., 2011 ) using the SIMM datasetterm downloaded from http://drive5.com/uchime/uchime_download.html . The dataset consists of 900 simulated chimerasterm that are approximately 250 bp long. The chimerasterm were generated from 2, 3 or 4 segments selected randomly from 86 original sequences and have similarities in the ranges 90–95%, 95–97% and 97–99% to the original sequences. They were either used unmodified or with 1–5% indels or 1–5% substitutions. We assessed the performance of (i) the original open-source UCHIME version 4.2 program, (ii) USEARCHterm version 7, (iii) USEARCHterm version 8, and (iv) VSEARCHterm. The results are shown in Table 1 and indicate that VSEARCHterm is superior to the other tools in almost all cases, and in particular when indels were added. The original UCHIME program was found to be quite, but also considerably slower than all the other tools. USEARCHterm was better than VSEARCHterm in only 4 out of 99 cases.

Next, we tested reference-based ( uchime_ref ) and de novo ( uchime_denovo ) chimeraterm detection using sequences from the 2011 version of Greengenesterm downloaded from http://greengenesterm.lbl.gov/Download/Sequence_Data/Fasta_data_files/ and from version 106 (May 2011) of the SILVAterm database downloaded from https://www.arb-silva.de/no_cache/download/archive/release_106/Exports/ . Sequences from the 16S rRNA V4 region was computationally extracted using the 515F (5′-GTGNCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) primers, and 8,000 reads were randomly selected from each database. PCR was simulated using a new simulation algorithm known as Simera ( NicholspersName & QuincepersName, 2016 ) (available at https://github.com/bnichols1979/Simera ) that includes amplification and creation of PCR artefacts like chimerasterm. We sampled 30,000 reads (-s 30,000) and generated 20,000 potential chimerasterm (-c 20,000). Defaults were used for other options to Simera. The output sequences were then fed into an Illumina MiSeq noise simulator ( Schirmer et al., 2015 ) ending up with 14,966 reads based on Greengenesterm and 14,952 reads based on SILVAterm, of which 1,262 and 1,640 reads contain chimeric sequences, respectively. Next, the sequences were either clustered using the cluster_fast command at 97% identity or dereplicated. VSEARCHterm and USEARCHterm version 7 and 8 were run using the uchime_denovo command and then using the uchime_ref command with the Gold database downloaded from http://drive5.com/uchime/uchime_download.html as the reference database. To assess the performance, the results were sorted based on the chimera scoreterm·var, and then the ability to classify individual sequences correctly into chimeric and non-chimeric was plotted as ROC curves. The curves reflect the accuracy of classifying individual reads, not clusters, as abundances were taken into account. The plots in Figs. 5 and 6 show that de novo chimera detectionterm performs better than reference-based detectionterm, with the SILVAterm dataset in particular, but it does of course depend on the reference database used. VSEARCHterm performs better than both versions of USEARCHterm for de novo chimera detectionterm. For reference-based detectionterm VSEARCHterm also performs better for the Greengenesterm dataset, while none of the programs work well with the SILVAterm dataset. Clustering at 97% appears to be more appropriate than dereplication. In this test, the USEARCHterm programs were about twice as fast as VSEARCHterm for de novo detection, while they were about 10–30% faster than VSEARCHterm for reference-based detectionterm.

Chimeraterm detection performance with the Greengenesterm dataset shown with ROC curves. USEARCHterm version 7 (blue) and 8 (orange) and VSEARCHterm (black) was run using the uchime_denovo (A, B) and the uchime_ref (C, D) commands on simulated Illumina data based on the Greengenesterm database that has either been clustered with a 97% identity threshold (using the cluster_fast command in VSEARCHterm) (A, C) or dereplicated (using the derep_fulllength command in VSEARCHterm) (B, D).

Chimeraterm detection performance on the SILVAterm dataset shown with ROC curves. USEARCHterm version 7 (blue) and 8 (orange) and VSEARCHterm (black) was run using the uchime_denovo (A, B) and the uchime_ref (C, D) commands on simulated Illumina data based on the SILVAterm database that has either been clustered with a 97% identity threshold (using the cluster_fast command in USEARCHterm) (A, C) or dereplicated (using the derep_fulllength command in VSEARCHterm) (B, D).

Merging of paired-end readsterm Evaluation of paired-end readsterm merging performance was carried out in a manner similar to that described for the evaluation of PEAR ( Zhang et al., 2014 ). We used whole genome sequencing data from Staphylococcus aureus subspecies aureus strain USA 300 TCH 1516 sequenced by MacCallum et al. (2009) and retrieved from the GAGE-B repository ( http://ccb.jhu.edu/gage_b/ ). The S.aureus reads were 101 bp long from on average 180 bp long fragments, giving a 45X coverage of the genome. We also used Methylococcus capsulatus strain Bath 16S rRNA V3 region amplicon reads sequenced by Masella et al. (2012) . These reads were 108 bp long and the pairs should have an overlap of exactly 18 bp. Merging options were set to allow a minimum overlap of 10 bp and a maximum of 5 mismatches (USEARCHterm 7 and 8 have different default values for those), while other options were left at defaults. All programs were run with 8 threads. Merged sequences that could be perfectly aligned to their respective reference sequences (either the entire genome or the specific rRNA region) using BWA MEM ( Li & Durbin, 2009 ) were considered correctly mergedterm. The results are shown in Table 2 . The numbers indicate that USEARCHterm version 7 merges the most reads for both bacteria, but also has the lowest percentage of correctly mergedterm pairs of those merged. USEARCHterm version 8 merges the fewest reads, but has the highest percentage of correctly mergedterm reads of those merged. VSEARCHterm is in the middle by merging more reads than USEARCHterm 8 with only a small decrease in the percentage of correct mergesterm·var. VSEARCHterm is about twice as fast as USEARCHterm 8 and 4–5 times faster than USEARCHterm version 7.

Paired-end readsterm merging performance. The number of sequence pairs, merged pairsterm·var, and correctly mergedterm pairs are shown for each bacterium and program. The percentage of reads merged, as well as the percentage of correctly mergedterm reads both of the merged reads and of all reads are also shown. Times are in seconds using 8 threads.

Subsampling performance. The observed distribution of the maximum amplicon abundance in 10,000 random subsamplings of 5% of the TARA V9 dataset using VSEARCHterm (A, black) and USEARCHterm version 8 (B, orange) is shown. The expected mean abundance is 782,133.5 (blue dashed line).

Subsampling We evaluated the subsampling commands of USEARCHterm version 8 and VSEARCHterm to check if the results obtained correspond to those expected. We performed 10,000 random subsamplings of 5% of the 9.5 million unique sequences in the TARA V9 dataset ( Karsenti et al., 2011 ). To make this possible with the 32-bit USEARCHterm, we first downsampled the dataset once to 10% using VSEARCHterm and then randomly subsampledterm it again at 50% with either USEARCHterm or VSEARCHterm. Plots of the distribution of the abundance of the most abundant sequence in each subsampling are shown in Fig. 7 . The highest amplicon abundance in the original dataset is 15,638,316. After the initial 10% subsampling, the highest abundance was 1,564,267. After the second subsampling, the top abundances should therefore have a distribution centred on a value of 782,133.5. As can be seen from the figure, the USEARCHterm distribution has a mean that is about 2,000 too small, while the VSEARCHterm distribution is correctly centred on the expected valueterm. Subsampling experiments were also performed at 2.5%, 1.5% and 0.5% with similar results, although the errors were of decreasing size. USEARCHterm seems to under-sample abundant ampliconsterm and to over-sample rare ampliconsterm.

Conclusions VSEARCHterm supports almost all of the commands and options for nucleotide sequence analysis in USEARCHterm version 7 as well as several new features. It has a 64-bit design and handles large datasets virtually only limited by the amount of available memoryterm. We have demonstrated that VSEARCHterm is in general more accurate than USEARCHterm when performing searching, clustering, chimeraterm detection and subsampling. The accuracy is on a par with USEARCHterm for paired-end readsterm merging. VSEARCHterm is faster than USEARCHterm when performing dereplication and merging of paired-end readsterm, but slower for clustering and chimeraterm detection. We will continue to improve the accuracy, speed and robustness of VSEARCHterm in the future, as well as adding new features.