VSEARCH: a versatile open source tool for metagenomics

View article
PeerJ
Note that a Preprint of this article also exists, first published September 1, 2016.

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 metagenomics studies produce many millions of sequences, and require open-source, fast and memory 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 USEARCH (Edgar, 2010), a set of tools designed and implemented by Robert C. Edgar, and available at http://drive5.com/usearch/. USEARCH offers a great number of commands and options to manipulate and analyse FASTQ and FASTA files. However, the source code of USEARCH 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 USEARCH, freely available to users without any memory limitation. VSEARCH includes most of the USEARCH functions in common use, and further development may add additional features. Here we describe the details of the VSEARCH implementation. To assess its performance in terms of speed and quality of results, we have evaluated some of the most important functions (searching, clustering, chimera detection and subsampling) and compared them to USEARCH. We find that VSEARCH delivers results that are better or on a par with USEARCH results.

Materials and Methods

Algorithms and implementation

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

Reading FASTA and FASTQ files

Most VSEARCH commands read files in FASTA or FASTQ format. The parser for FASTQ files in VSEARCH is compliant with the standard as described by Cock et al. (2010) and correctly parses all their tests files. FASTA and FASTQ 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 VSEARCH, allowing for instance many separate FASTA files to be piped into VSEARCH for simultaneous dereplication, or allowing the creation of complex pipelines without ever having to write on slow disks.

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

Writing result files

VSEARCH can output results in a variety of formats (FASTA, FASTQ, 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 comparison is a core functionality of VSEARCH. Several commands compare a query sequence against a database of sequences: all-vs-all alignment (allpairs_global), clustering (cluster_fast, cluster_size, cluster_smallmem), chimera detection (uchime_denovo and uchime_ref) and searching (usearch_global). This comparison function proceeds in two phases: an initial heuristic filtering based on shared words, followed by optimal alignment of the query with the most promising candidates.

The first phase is presumably quite similar to USEARCH (Edgar, 2010). Heuristics are used to identify a small set of database sequences that have many words in common with the query sequence. Words (or k-mers) consist of a certain number k of consecutive nucleotides of a sequence (8 by default, adjustable with the wordlength option). All overlapping words are included. A sequence of length n then contains at most n − k + 1 unique words. VSEARCH counts the number of shared words between the query and each database sequence. Words that appear multiple times are counted only once. To count the words in the database sequences quickly, VSEARCH creates an index of all the 4k possible distinct words and stores information about which database sequences they appear in. For extremely frequent words, the set of database sequences is represented by a bitmap; otherwise the set is stored as a list. A finer control of k-mer indexing is described for USEARCH by the pattern (binary string indicating which positions must match) and slots options. USEARCH has such options but seems to ignore them. Currently, VSEARCH ignores these two options too. The minimum number of shared words 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 words.

Comparing sequences based on statistics of shared words is a common method to quickly assess the similarity between two sequences without aligning them, which is often time-consuming. The D2 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 D2 statistic, but multiple matches of the same word are ignored.

In the second phase, searching proceeds by considering the database sequences in a specific order, starting with the sequence having the largest number of words in common with the query, and proceeding with a decreasing number of shared words. If two database sequences have the same number of words in common with the query, the shortest sequence is considered first. The query sequence is compared with each database sequence by computing the optimal global alignment. The alignment is performed using a multi-threaded and vectorised full dynamic programming algorithm (Needleman & Wunsch, 1970) adapted from SWIPE (Rognes, 2011). Due to the extreme memory requirements of this method when aligning two long sequences, an alternative algorithm 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 algorithm uses only a linear amount of memory but is considerably slower. This second phase is probably where USEARCH and VSEARCH differ the most, as USEARCH 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 programming) 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 programming implementation in VSEARCH 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 option, the database sequence 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 USEARCH 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 similarity and presented as the search results.

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

Clustering

VSEARCH includes commands to perform de novo clustering using a greedy and heuristic centroid-based algorithm with an adjustable sequence similarity threshold 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 similarity may be found and considered. By default, the query is clustered with the centroid presenting the highest sequence similarity (distance-based greedy clustering, DGC), or, if the sizeorder option is turned on, the centroid with the highest abundance (abundance-based greedy clustering, AGC) (He et al., 2015; Westcott & Schloss, 2015; Schloss, 2016). VSEARCH 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, VSEARCH 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 rereplication

Full-length dereplication (derep_fulllength) is performed using a hash table 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 table is initially empty. For each input sequence, the hash is computed and a lookup in the hash table 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 table.

Prefix dereplication (derep_prefix) is also implemented. As with full-length dereplication, 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. VSEARCH 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 dereplication, VSEARCH first creates an initially empty hash table. 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 table, the input sequence is clustered with the matching hash table 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 table, the sequences are clustered, the matching sequence is deleted from the hash table and the full input sequence is inserted into the hash table instead. If no match is found for any prefix, the full sequence is inserted into the hash table. In the end, the remaining sequences in the hash table will be output with accumulated abundances for all sequences in each cluster.

In order to identify matches in the hash table during prefix dereplication, 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 function (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. VSEARCH exclusively also offers the possibility of relabelling each sequence with the SHA-1 (Eastlake & Jones, 2001) or MD5 (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. VSEARCH includes public domain code for the MD5 algorithm written by Alexander Peslyak, and for SHA1 by Steve Reid and others.

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

Chimera detection

Chimeras are detected either de novo (uchime_denovo command) or with a reference database (uchime_ref command) using the UCHIME algorithm described by Edgar et al. (2011). VSEARCH 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. VSEARCH 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 chimera detection.

Low-complexity sequence masking

VSEARCH includes a highly optimized and parallelized implementation of the Dust algorithm by Tatusov and Lipman for masking of simple repeats and low-complexity nucleotide sequences. It is considerably faster than the implementation of the same algorithm in USEARCH. Their code available at ftp://ftp.ncbi.nlm.nih.gov/pub/tatusov/dust/version1/src/ is in the public domain. VSEARCH uses this algorithm by default, while USEARCH by default uses an undocumented rapid masking algorithm called fastnucleo. VSEARCH performs soft-masking automatically for the pairwise alignment, search, clustering and chimera 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.

FASTQ file processing

VSEARCH includes commands to detect the FASTQ file version and the range of quality scores used (fastq_chars), as well as two commands for computing sequence quality statistics (fastq_stats and fastq_eestats). It can also truncate and filter sequences in FASTQ files based on various criteria (fastq_filter). A new command is added to convert between different FASTQ file versions and quality encodings (fastq_convert), e.g., from the old Phred+64 encoded Illumina FASTQ files to the newer Phred+33 format.

Merging of paired-end reads

Merging of paired-end reads is supported by VSEARCH using the fastq_mergepairs command. The method used has some similarity to PEAR (Zhang et al., 2014) and recognises options similar to USEARCH. The algorithm computes the optimal ungapped alignment of the overlapping region of the forward sequence and the reverse-complemented reverse sequence. The alignment requires a minimum overlap length (specified with the fastq_minovlen option, default 10), a maximum number of mismatches (fastq_maxdiffs option, default 5), and a minimum and maximum length of the merged sequence (fastq_minmergelen option, default 1, and fastq_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 fastq_allowmergestagger option. VSEARCH 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 VSEARCH assumes all nucleotide background frequencies are 0.25. When merging sequences, VSEARCH computes posterior quality scores 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

VSEARCH can sort FASTA files by decreasing sequence length (sortbylength) or abundance (sortbysize). VSEARCH can also perform shuffling of FASTA files in random order (shuffle). A seed value for the pseudo random number generator may be provided by the randseed option to obtain replicable results.

Subsampling

Sequences in FASTA and FASTQ files can be subsampled (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, subsampled and then dereplicated.

Results and Discussion

Supported commands and options

VSEARCH implements the following commands available in USEARCH version 7: allpairs_global, cluster_fast, cluster_smallmem, derep_fulllength, derep_prefix, fastq_chars, fastq_filter, fastq_mergepairs, fastq_stats, fastx_mask, maskfasta, sortbylength, sortbysize, uchime_denovo, uchime_ref and usearch_global. In addition, the following commands available in USEARCH version 8 have been implemented: fastq_eestats, fastx_revcomp, fastx_subsample and search_exact. VSEARCH additionally includes a few new commands that do not exist in USEARCH: cluster_size, fastq_convert, rereplicate and shuffle.

Some USEARCH version 7 commands have not yet been implemented in VSEARCH. 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_usearch, udb2fasta, udbinfo, udbstats), and the UPARSE pipeline (cluster_otus, uparse_ref).

Almost all USEARCH 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 USEARCH version 7 have generally been used in order to make VSEARCH an almost drop-in replacement. In fact, in QIIME many commands will run fine if an alias or link from usearch to vsearch is made. Detailed documentation of VSEARCH is available as a man page. We will consider adding further commands and options to VSEARCH in the future.

Performance assessment

The performance of the most important functions of VSEARCH version 2.0.3 (64-bit) was evaluated and compared to the 32-bit versions of USEARCH version 7.0.1090 and 8.1.1861. All datasets used were small enough to fit comfortably in the memory allocated to a 32-bit process. Chimera 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 USEARCH paper (Edgar, 2010), its supplementary, and on the website (http://drive5.com/usearch/benchmark_rfam.html), by assessing the ability of the programs to identify RNA sequences belonging to the same family in RFAM (Burge et al., 2013). The 383,004 sequences in Rfam 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 rate between 0.010 and 0.015, VSEARCH is more accurate than USEARCH’s latest version. For lower values, the three programs have similar accuracies. At higher false discovery rates, USEARCH version 8 has an advantage.

Search accuracy on the RFAM v11 dataset.
Figure 1: Search accuracy on the RFAM v11 dataset.
USEARCH version 7 (blue), USEARCH version 8 (orange) and VSEARCH (black) was run using the usearch_global command on subsets of the RFAM dataset to identify members of the same families. The plot shows the true positive rate (also known as the recall or sensitivity) as a function of the false discovery rate at varying sequence similarity levels. This curve is based on data from 20 shufflings of the dataset.

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

Clustering

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

We independently evaluated the clustering accuracy of USEARCH and VSEARCH as described for Swarm (Mahé 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 usearch_global command of USEARCH. 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 recall, precision and the adjusted Rand index of the classifications. Recall measures to what extent amplicons assigned to the same species are grouped together in the same OTU (i.e., not over-splitting). Precision measures to what extent amplicons in an OTU are assigned to the same species (i.e., not over-grouping). The adjusted Rand index (Rand, 1971; Hubert & Arabie, 1985) summarises both precision and recall using a measure of agreement between two clusterings while adjusting for expected values 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 VSEARCH and USEARCH version 8 is negligible. The difference is larger for length-sorted sequences. When using length sorting, USEARCH 8 (as well as version 7 on the even dataset) shows better precision than VSEARCH 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 index is maximised at 95–97% similarity. This is precisely the region where for length sorting at least VSEARCH outperforms USEARCH.

Clustering accuracy on the even dataset.
Figure 2: Clustering accuracy on the even dataset.
USEARCH version 7 (blue) and 8 (orange) and VSEARCH (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 index (A, B), recall (C, D) and precision (E, F) metrics.
Clustering accuracy on the uneven dataset.
Figure 3: Clustering accuracy on the uneven dataset.
USEARCH version 7 (blue) and 8 (orange) and VSEARCH (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 index (A, B), recall (C, D) and precision (E, F) metrics.

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

Clustering speed.
Figure 4: Clustering speed.
Median wall time in seconds to cluster the even (A, B) and uneven (C, D) datasets using USEARCH version 7 (blue) and 8 (orange) and VSEARCH (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 dereplication (derep_fulllength) VSEARCH was about 40–50% faster than USEARCH version 7 and 50–70% faster than version 8 on all three datasets. All programs were approximately equally fast on prefix dereplication (derep_prefix) of the even and uneven datasets. However, prefix dereplication of the BioMarKs dataset was extremely slow with USEARCH. USEARCH version 7 used more than 4 min and version 8 more than 27 min, while VSEARCH used less than 4 s. The prefix dereplication algorithm used in USEARCH 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 USEARCH version 7 and 8 down to just 5 and 6 s, respectively.

Chimera detection

We evaluated the chimera detection accuracy of VSEARCH and USEARCH in two ways, first using a method similar to that performed for UCHIME, and then using a new chimera simulation procedure based on sequences from Greengenes (DeSantis et al., 2006) and SILVA (Quast et al., 2013) sequences.

Table 1:
Chimera detection performance with the SIMM dataset.
UCHIME (UC), USEARCH version 7 (U7) and 8 (U8), and VSEARCH (V) was run using the uchime_ref algorithm on the SIMM dataset that was originally also used to evaluate the UCHIME algorithm. 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 chimeras detected out of 100 of each type is shown. The best results in each category are shaded.
2 segments 3 segments 4 segments
Divergence Noise UC U7 U8 V UC U7 U8 V UC U7 U8 V
97–99% 89 88 88 89 56 52 52 55 38 33 34 35
i1 79 79 77 85 46 44 43 53 32 27 24 34
i2 64 57 56 77 33 32 31 56 24 20 18 33
i3 48 45 36 72 37 35 29 45 16 17 16 21
i4 29 24 23 65 18 11 13 40 9 9 8 25
i5 27 22 16 53 15 12 12 39 7 8 6 17
m1 83 83 83 81 53 48 48 53 33 29 29 30
m2 73 71 71 72 49 44 44 50 28 22 22 27
m3 66 66 66 68 40 40 39 44 21 20 21 21
m4 55 54 53 57 28 24 23 28 21 18 18 19
m5 44 44 42 48 20 19 18 28 16 14 12 12
95–97% 100 100 100 100 80 77 76 79 64 60 59 63
i1 100 98 98 100 77 75 72 75 54 55 53 61
i2 96 94 93 99 60 55 55 71 48 44 44 60
i3 86 82 82 95 61 50 52 70 38 36 31 53
i4 75 66 64 95 48 41 39 64 29 29 22 47
i5 64 58 53 86 37 32 25 60 24 19 19 46
m1 99 99 99 99 76 73 73 76 60 57 57 60
m2 98 97 97 97 71 69 69 71 50 48 46 48
m3 93 94 94 96 63 61 61 64 41 41 41 42
m4 92 92 90 93 56 55 54 57 39 39 37 41
m5 86 86 85 86 53 51 51 56 35 35 34 34
90–95% 100 100 100 100 93 93 93 93 88 88 88 86
i1 100 100 100 100 88 88 87 91 86 86 87 88
i2 99 97 99 99 83 79 78 88 74 72 72 84
i3 100 100 100 100 79 76 75 88 74 69 70 82
i4 99 94 96 99 80 71 72 84 66 62 61 79
i5 95 84 86 99 74 65 65 88 55 48 48 71
m1 100 100 100 100 89 89 89 92 87 87 86 85
m2 100 100 100 100 87 87 87 89 78 78 78 79
m3 100 99 99 100 86 86 86 89 76 76 78 80
m4 100 100 100 100 82 82 84 83 73 73 72 78
m5 99 98 98 99 82 81 82 84 75 73 75 79
DOI: 10.7717/peerj.2584/table-1

First we repeated the evaluation of the uchime_ref command described in the UCHIME paper (Edgar et al., 2011) using the SIMM dataset downloaded from http://drive5.com/uchime/uchime_download.html. The dataset consists of 900 simulated chimeras that are approximately 250 bp long. The chimeras 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) USEARCH version 7, (iii) USEARCH version 8, and (iv) VSEARCH. The results are shown in Table 1 and indicate that VSEARCH 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. USEARCH was better than VSEARCH in only 4 out of 99 cases.

Next, we tested reference-based (uchime_ref) and de novo (uchime_denovo) chimera detection using sequences from the 2011 version of Greengenes downloaded from http://greengenes.lbl.gov/Download/Sequence_Data/Fasta_data_files/ and from version 106 (May 2011) of the SILVA 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 (Nichols & Quince, 2016) (available at https://github.com/bnichols1979/Simera) that includes amplification and creation of PCR artefacts like chimeras. We sampled 30,000 reads (-s 30,000) and generated 20,000 potential chimeras (-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 Greengenes and 14,952 reads based on SILVA, 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. VSEARCH and USEARCH 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 score, 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 detection performs better than reference-based detection, with the SILVA dataset in particular, but it does of course depend on the reference database used. VSEARCH performs better than both versions of USEARCH for de novo chimera detection. For reference-based detection VSEARCH also performs better for the Greengenes dataset, while none of the programs work well with the SILVA dataset. Clustering at 97% appears to be more appropriate than dereplication. In this test, the USEARCH programs were about twice as fast as VSEARCH for de novo detection, while they were about 10–30% faster than VSEARCH for reference-based detection.

Chimera detection performance with the Greengenes dataset shown with ROC curves.
Figure 5: Chimera detection performance with the Greengenes dataset shown with ROC curves.
USEARCH version 7 (blue) and 8 (orange) and VSEARCH (black) was run using the uchime_denovo (A, B) and the uchime_ref (C, D) commands on simulated Illumina data based on the Greengenes database that has either been clustered with a 97% identity threshold (using the cluster_fast command in VSEARCH) (A, C) or dereplicated (using the derep_fulllength command in VSEARCH) (B, D).
Chimera detection performance on the SILVA dataset shown with ROC curves.
Figure 6: Chimera detection performance on the SILVA dataset shown with ROC curves.
USEARCH version 7 (blue) and 8 (orange) and VSEARCH (black) was run using the uchime_denovo (A, B) and the uchime_ref (C, D) commands on simulated Illumina data based on the SILVA database that has either been clustered with a 97% identity threshold (using the cluster_fast command in USEARCH) (A, C) or dereplicated (using the derep_fulllength command in VSEARCH) (B, D).

Merging of paired-end reads

Evaluation of paired-end reads 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 (USEARCH 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 merged. The results are shown in Table 2. The numbers indicate that USEARCH version 7 merges the most reads for both bacteria, but also has the lowest percentage of correctly merged pairs of those merged. USEARCH version 8 merges the fewest reads, but has the highest percentage of correctly merged reads of those merged. VSEARCH is in the middle by merging more reads than USEARCH 8 with only a small decrease in the percentage of correct merges. VSEARCH is about twice as fast as USEARCH 8 and 4–5 times faster than USEARCH version 7.

Table 2:
Paired-end reads merging performance.
The number of sequence pairs, merged pairs, and correctly merged pairs are shown for each bacterium and program. The percentage of reads merged, as well as the percentage of correctly merged reads both of the merged reads and of all reads are also shown. Times are in seconds using 8 threads.
Bacterium Program Pairs Merged Correct %Merged %Cor/Mer %Cor/All Time (s)
Staphylococcus aureus USEARCH 7 647,052 273,438 270,849 42.26 99.05 41.86 11.65
USEARCH 8 647,052 203,729 202,003 31.49 99.15 31.22 4.69
VSEARCH 647,052 214,988 213,103 33.23 99.12 32.93 2.15
Methylococcus capsulatus strain Bath USEARCH 7 673,845 643,903 642,720 95.56 99.82 95.38 14.43
USEARCH 8 673,845 554,099 553,747 82.23 99.94 82.18 6.27
VSEARCH 673,845 581,752 581,346 86.33 99.93 86.27 3.61
DOI: 10.7717/peerj.2584/table-2
Subsampling performance.

Figure 7: Subsampling performance.

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

Subsampling

We evaluated the subsampling commands of USEARCH version 8 and VSEARCH 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 USEARCH, we first downsampled the dataset once to 10% using VSEARCH and then randomly subsampled it again at 50% with either USEARCH or VSEARCH. 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 USEARCH distribution has a mean that is about 2,000 too small, while the VSEARCH distribution is correctly centred on the expected value. Subsampling experiments were also performed at 2.5%, 1.5% and 0.5% with similar results, although the errors were of decreasing size. USEARCH seems to under-sample abundant amplicons and to over-sample rare amplicons.

Conclusions

VSEARCH supports almost all of the commands and options for nucleotide sequence analysis in USEARCH 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 memory. We have demonstrated that VSEARCH is in general more accurate than USEARCH when performing searching, clustering, chimera detection and subsampling. The accuracy is on a par with USEARCH for paired-end reads merging. VSEARCH is faster than USEARCH when performing dereplication and merging of paired-end reads, but slower for clustering and chimera detection. We will continue to improve the accuracy, speed and robustness of VSEARCH in the future, as well as adding new features.

6241 Citations   Views   Downloads