Development manual for vsearch 2.33.0, built from the dev branch: it describes changes that are not released yet. The manual for the current release is published separately, next to every published version.

NAME

vsearch --search_oligodb — find every occurrence of oligos (primers, tags, barcodes) in sequences

SYNOPSIS

vsearch --search_oligodb fastxfile --db fastafile (--alnout | --blast6out | --userout) filename [options]

DESCRIPTION

The vsearch command --search_oligodb looks for short oligonucleotides (primers, tags, barcodes, adapters, splints), given in a fasta file with --db, in the sequences of a fasta or fastq file, and reports every occurrence of every oligo, with its position, on both strands. It was written for long reads that carry several copies of a target (rolling-circle amplification, concatemers, a barcode at each end), but works on sequences of any length, from a few nucleotides to hundreds of kilobases. The sequences are read one after the other and never held in memory all at once.

The oligo is aligned end to end, and the sequence is free at both ends (a semi-global alignment). An occurrence is reported when its alignment has at most --maxdiffs differences (2 by default), a difference being a substitution, an inserted nucleotide or a deleted nucleotide. Among the alignments with the fewest differences, the one with the fewest gap openings is reported. --maxgaps bounds the gap openings, and --maxgaps 0 restricts the search to substitutions.

Nucleotides are compared by overlap of their IUPAC base sets, on both sides: an oligo Y matches a C or a T in the sequence, a sequence R matches an oligo A, G, R, N… In particular, an N in the sequence matches any nucleotide. As a consequence, a run of Ns at least as long as an oligo is an occurrence without difference of that oligo, and a site with a few Ns needs fewer differences to be reported. Use --n_mismatch to count any position involving an N as a mismatch. Case is ignored, and U is read as T.

An occurrence of an oligo is one site: the alignments ending near one another (closer than half the oligo length) are collapsed into the one with the fewest differences, the leftmost in case of a tie. Occurrences of different oligos are all reported, even when they overlap.

Occurrences truncated by a sequence end (a barcode on a read that starts or ends in the middle of it, for instance) are reported, as long as at least --target_cov of the oligo (0.75 by default) lies inside the sequence. The part of the oligo off the end is not counted as differences; tlo > 1 or thi < tl shows which end is missing.

Occurrences are written in the order of the sequences, and within a sequence by position. This command is multi-threaded, and its output is the same whatever the number of threads.

Positions

For this command, the userfields qlo and qhi report the first and last positions of the occurrence in the sequence (the query), and tlo and thi the first and last positions of the oligo (the target) that are aligned. The columns 7 to 10 of --blast6out report the same four values. For an occurrence on the minus strand, qlo and qhi are positions on the plus strand, swapped (qlo > qhi), as vsearch does for all minus-strand hits. The alignment fields caln, aln, qrow and trow cover the occurrence only. (In the other commands, where the alignments are global, qlo, qhi, tlo and thi span the whole sequences.) See vsearch-userfields(7).

The field raw reports the number of differences of the alignment (as diffs does), bits is 0 and evalue is -1.

Random occurrences

The number of occurrences found by chance grows very fast with --maxdiffs. Occurrences of one primer in uniform random sequence, per 100 kb of sequence searched on both strands:

primer --maxdiffs 2 3 4
515F, 19 nt, 2 ambiguous positions 0.02 0.44 6.9
515F, 19 nt, no ambiguous position 0.00 0.13 2.6
806R, 20 nt, 3 ambiguous positions 0.02 0.55 8.4

Real sequences are not random (low-complexity regions, repeats), so these rates are a floor, not an estimate. For a 20-nucleotide primer, 3 differences is a practical ceiling on long reads, and each ambiguous position multiplies the random occurrences by 2 to 3. When possible, search for a longer oligo (a primer and its tag together, for instance), and use the diffs field to filter the output rather than a large --maxdiffs.

Truncated occurrences are the other source: the part of an oligo left inside a read can be short enough to match by chance. On a nanopore run of 1.7 million reads searched for 96 barcodes of 24 nucleotides (5 of them present) with --maxdiffs 3, truncated occurrences of the 91 absent barcodes, which can only be random, numbered 2,456,380 with --target_cov 0.5, 1,181 with 0.75, and none with 1; the present barcodes had 340,343, 38,879 and none. The whole-length occurrences of the absent barcodes numbered 26 in all three cases.

Differences with usearch

usearch has a search_oligodb command. vsearch’s version differs on purpose:

  • gaps are allowed (usearch aligns without gaps, at any -maxdiffs); --maxgaps 0 gives usearch’s substitution-only matching;
  • an N in the sequence always matches, where usearch drops a site that holds two Ns or more, whatever -maxdiffs;
  • occurrences truncated by a sequence end are reported;
  • occurrences are ordered by position within a sequence, not by score;
  • on the minus strand, qlo and qhi are swapped (usearch reports them in increasing order, and swaps the oligo positions instead);
  • --blast6out reports an e-value of -1 and a bit score of 0;
  • the default masking is none, and the default strand is both.

Algorithm

Each sequence is read once, and searched for all the oligos, on both strands, in the same pass. The search is bit-parallel: the state of an oligo fits in one 64-bit machine word, one bit per nucleotide of the oligo, which is why oligos are limited to 64 nucleotides. IUPAC ambiguity codes cost nothing: for each nucleotide symbol, the positions of an oligo it matches are computed once, as a bit mask.

By default, the search is the bit-vector algorithm of Myers (1999), which gives, at every position of the sequence, the fewest differences (substitutions, insertions and deletions) of an alignment of the whole oligo ending there. Its cost does not depend on --maxdiffs. With --maxgaps 0, the search is the shift-and algorithm (Baeza-Yates and Gonnet 1992), extended to mismatches by Wu and Manber (1992), with one machine word per number of mismatches, from 0 to --maxdiffs: its cost grows with --maxdiffs.

Both algorithms only report where an occurrence ends, and with how many differences. Each occurrence is then aligned again, by dynamic programming over a short window of the sequence, to recover its start and its alignment, and to choose, among the alignments with the fewest differences, the one with the fewest gap openings.

Baeza-Yates R, Gonnet GH. (1992) A new approach to text searching. Communications of the ACM 35(10):74-82 doi: 10.1145/135239.135243

Myers G. (1999) A fast bit-vector algorithm for approximate string matching based on dynamic programming. Journal of the ACM 46(3):395-415 doi: 10.1145/316542.316550

Wu S, Manber U. (1992) Fast text searching: allowing errors. Communications of the ACM 35(10):83-91 doi: 10.1145/135239.135244

OPTIONS

mandatory options

--db must be specified, along with at least one output option.

--db filename
Search for the oligos (primers, tags, barcodes…) in filename, in fasta format. Oligos can contain IUPAC ambiguity codes (see vsearch-nucleotides(7)), and must be 1 to 64 nucleotides long: an empty oligo, or an oligo longer than 64 nucleotides, is a fatal error. A UDB database is not accepted.

core options

--maxdiffs positive integer
Report an occurrence only if its alignment has at most integer differences: substitutions plus gap columns (inserted or deleted nucleotides). The default is 2. Any value up to the length of the oligo is accepted, but the number of random occurrences grows very fast with it (see DESCRIPTION). Oligo nucleotides hanging off a read end are not differences (see --target_cov). The diffs userfield reports the quantity compared against.
--maxgaps positive integer

Report an occurrence only if its alignment has at most integer gap openings: a run of consecutive gap columns counts as one, whatever its length (--maxdiffs bounds the number of gap columns). The default is unlimited. The alignment compared against --maxgaps is the one with the fewest differences, and among those, the one with the fewest gap openings; an occurrence is rejected when that alignment has too many gap openings, even if an alignment without gaps but with more substitutions would fit within --maxdiffs.

--maxgaps 0 is the exception: it switches to substitutions only. The search then looks for occurrences with up to --maxdiffs mismatches and no gap at all, as usearch does, and gapped alignments are not considered.

--n_mismatch

Count alignments of nucleotides against Ns as mismatches. By default, an alignment column holding an N is neutral: it scores zero, and it is counted as a matching column when the identity percentage is computed. With --n_mismatch, any column where at least one of the two symbols is an N (regardless of case) is scored and counted as a mismatch instead, N against N included. Both the alignment score and the post-alignment count of matches and mismatches are affected, so the identity percentage compared against --id (and reported in the output files) changes too.

The option targets N only: the other ambiguous symbols (BDHKMRSVWY) keep their default behaviour, and still count as matching any symbol they share a nucleotide with. See vsearch-pairwise_alignment_parameters(7) for the default treatment of ambiguous symbols, and for why a query can align to a long run of Ns with 100% identity.

--strand plus|both
Search the oligos on both strands (default), or on the plus strand only, i.e. only as they are written in the --db file. Keywords are case-insensitive. Unlike the other search commands, the default here is both.
--target_cov real

Report an occurrence truncated by a read end only if at least real of the oligo lies inside the read. The default is 0.75. The oligo nucleotides hanging off the read end cost nothing (they are not counted in diffs), but the shorter the part left inside the read, the more often it matches by chance: at 0.5, random partial matches outnumber real truncated oligos (see DESCRIPTION). Use --target_cov 1 to report whole occurrences only.

Unlike the other commands, where the target coverage counts only the target nucleotides aligned with a query nucleotide, the fraction compared against here is the span of the oligo inside the read, deleted nucleotides included: a deletion is a difference, bounded by --maxdiffs, not a truncation. The tcov userfield keeps its usual definition, so it is below 100% for an occurrence with a deletion even when the whole oligo lies inside the read; tlo > 1 or thi < tl tells a truncated occurrence apart.

--threads positive integer
Set the number of computation threads to use, from 1 to 1024. The number of threads should not exceed the number of available CPU cores. The value 0 is also accepted and, like the default, uses all available cores; decimal values are truncated to their integer part. On Linux, “available” accounts for the CPU affinity mask and the cgroup CPU quota of the running process, so a job confined by taskset(1), Slurm, Docker or Kubernetes launches one thread per core it was actually granted, rather than one per core the machine has.

output options

--alnout filename
Write the alignment of each occurrence to filename in a human-readable format, with its positions in the query and in the oligo, to check occurrences by eye (while choosing a value for --maxdiffs, for instance). Use --rowlen to set the alignment line width.
--blast6out filename
Write the occurrences to filename using a BLAST-like tab-separated format, one line per occurrence, with twelve fields: query label, oligo label, percentage identity, alignment length, mismatches, gap openings, first and last query positions of the occurrence (swapped on the minus strand), first and last oligo positions aligned, expectation value (always -1), and bit score (always 0). The positions are those of the qlo, qhi, tlo and thi userfields.
--userout filename
Write user-defined tab-separated output to filename. Select and order the fields with --userfields.

secondary options

--bzip2_decompress
Specify that the input pipe is streaming data compressed using Huffman coding. See bzip2(1) for more details. This option is required when compressed data arrives on standard input through a pipe (‘-’), where the format cannot be detected without consuming the stream. It is not needed when reading from a regular file compressed with bzip2, nor when such a file is redirected to standard input: compression is then detected automatically, and a contradicting option is ignored (with a warning when the input is standard input). Pipes other than standard input, such as shell process substitutions and named FIFOs, are always read as uncompressed data; compressed data must arrive on standard input or as a named file.
--gzip_decompress
Specify that the input pipe is streaming data compressed using Lempel-Ziv coding. See gzip(1) for more details. This option is required when compressed data arrives on standard input through a pipe (‘-’), where the format cannot be detected without consuming the stream. It is not needed when reading from a regular file compressed with gzip, nor when such a file is redirected to standard input: compression is then detected automatically, and a contradicting option is ignored (with a warning when the input is standard input). Pipes other than standard input, such as shell process substitutions and named FIFOs, are always read as uncompressed data; compressed data must arrive on standard input or as a named file.
--hardmask

Replace masked nucleotides with Ns, rather than lowercasing them.

This is also the only masking option that reaches the alignment. Soft and dust masking (see --qmask and --dbmask) only keep masked words out of the k-mer pre-filter that selects candidates; the masked region is still aligned and scored like any other. An N, on the other hand, scores zero and counts as a matching column, so hard masking makes a masked region match whatever it is aligned against, and the reported identity can only go up. Add --n_mismatch to count those columns as mismatches instead.

--log filename
Write messages to filename. Messages include program version, start and finish times, elapsed time, amount of memory available, maximum amount of memory consumed, number of cores and command line options, and if need be, command-specific informational messages, warnings, and errors. Messages are written in addition to the standard error, not instead of it; use --quiet to silence the standard error copy (see the --quiet entry for what it spares).
--no_progress
Suppress the gradually increasing progress indicator normally written to the standard error stderr(3).
--notrunclabels
Retain whole sequence headers in output files. By default, vsearch truncates sequence headers at first space or tabulation. This option suppresses truncation.
--qmask none|dust|soft
Mask regions in the query sequences using the dust method or the soft method, or none to suppress masking. Values are case-insensitive. The default is none, unlike the other search commands: the search ignores case, so soft masking has no effect on it, and only matters with --hardmask, which turns the masked nucleotides into Ns. Beware that an N matches any oligo nucleotide (see --n_mismatch), so hard masking creates occurrences in the masked regions rather than hiding them.
--quiet
Suppress messages to the standard output stdout(3) and standard error stderr(3), except for warnings and error messages.
--rowlen positive integer
Set the width of alignment lines in --alnout output. The default value is 64. Set to 0 to disable line wrapping.
--sizein
Use the abundance annotations present in sequence headers when reading fasta or fastq file. Search for the pattern [>@;]size=integer[;]. Entries without abundance annotations are silently assumed to be of size=1.
--userfields string
Select and order the fields written to --userout output. Fields are separated by + (e.g. query+target+id). See vsearch-userfields(7) for a complete description of all available fields.

EXAMPLES

Find every occurrence of a set of primers and barcodes in nanopore reads, with up to three differences, and write one line per occurrence:

vsearch \
    --search_oligodb reads.fastq.gz \
    --db oligos.fasta \
    --maxdiffs 3 \
    --userout hits.tsv \
    --userfields query+target+qstrand+qlo+qhi+tlo+thi+diffs

The same search, without gaps, as usearch would do it, and with the alignments written out for checking:

vsearch \
    --search_oligodb reads.fastq.gz \
    --db oligos.fasta \
    --maxdiffs 3 \
    --maxgaps 0 \
    --alnout hits.aln \
    --blast6out hits.tsv

Convert the occurrences to BED (0-based start, end excluded, both positions on the plus strand):

vsearch \
    --search_oligodb reads.fastq \
    --db oligos.fasta \
    --userout hits.tsv \
    --userfields query+qlo+qhi+target+diffs+qstrand

awk 'BEGIN {OFS = "\t"}
     {lo = ($2 < $3) ? $2 : $3 ; hi = ($2 < $3) ? $3 : $2 ;
      split($1, label, " ") ;
      print label[1], lo - 1, hi, $4, $5, $6}' hits.tsv > hits.bed

SEE ALSO

vsearch-cut(1), vsearch-search_exact(1), vsearch-usearch_global(1), vsearch-cigar(5), vsearch-fasta(5), vsearch-fastq(5), vsearch-nucleotides(7), vsearch-usearch(7), vsearch-userfields(7)

CITATION

Rognes T, Flouri T, Nichols B, Quince C, Mahé F. (2016) VSEARCH: a versatile open source tool for metagenomics. PeerJ 4:e2584 doi: 10.7717/peerj.2584

REPORTING BUGS

Submit suggestions and bug-reports at https://github.com/torognes/vsearch/issues, send a pull request on https://github.com/torognes/vsearch, or compose a friendly or curmudgeont e-mail to Torbjørn Rognes (torognes@ifi.uio.no).

AVAILABILITY

Source code and binaries are available at https://github.com/torognes/vsearch.

These manual pages are also published online at https://torognes.github.io/vsearch/. That address always serves the manual of the current release; the manuals of earlier releases are listed at https://torognes.github.io/vsearch/versions/.

COPYRIGHT

Copyright (C) 2014-2026, Torbjørn Rognes, Frédéric Mahé and Tomás Flouri

All rights reserved.

Contact: Torbjørn Rognes torognes@ifi.uio.no, Department of Informatics, University of Oslo, PO Box 1080 Blindern, NO-0316 Oslo, Norway

This software is dual-licensed and available under a choice of one of two licenses, either under the terms of the GNU General Public License version 3 or the BSD 2-Clause License.

GNU General Public License version 3

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

The BSD 2-Clause License

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

  1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

  2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

ACKNOWLEDGMENTS

We would like to thank the authors of the following projects for making their source code available:

  • vsearch includes code from Google’s CityHash project by Geoff Pike and Jyrki Alakuijala, providing some excellent hash functions available under a MIT license.
  • vsearch includes code derived from Tatusov and Lipman’s DUST program that is in the public domain.
  • vsearch includes public domain code written by Alexander Peslyak for the MD5 message digest algorithm.
  • vsearch includes public domain code written by Steve Reid and others for the SHA1 message digest algorithm.
  • vsearch binaries may include code from the zlib library, copyright Jean-Loup Gailly and Mark Adler.
  • vsearch binaries may include code from the bzip2 library, copyright Julian R. Seward.
This site uses Just the Docs, a documentation theme for Jekyll.