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

NAME

vsearch --fastq_mergepairs — merge paired-end reads into one sequence

SYNOPSIS

vsearch --fastq_mergepairs fwdfile --reverse revfile (--fastaout | --fastqout | other output options) filename [options]

DESCRIPTION

The vsearch command --fastq_mergepairs merges paired-end sequence reads into a single sequence by aligning the forward and reverse reads and combining their overlapping regions. The forward reads are specified as the argument to this option; the reverse reads are specified with --reverse. Reads are matched by position: the first forward read is paired with the first reverse read, the second with the second, and so on. Labels are not used for matching; the run stops with a fatal error if the two files contain different numbers of reads.

The reverse read is reverse-complemented before alignment. Merging requires an overlap between the two reads of at least --fastq_minovlen bases (default 10, minimum 5). Read pairs with too many mismatches in the overlap — more than --fastq_maxdiffs (default 10) or more than --fastq_maxdiffpct percent (default 100.0%) — are discarded. Additional heuristics prevent merging of read pairs that cannot be aligned reliably.

In the merged region, quality scores from the two reads are combined using the Phred score formula. Outside the overlap, the quality scores from the contributing read are used directly. The quality score of N bases is replaced by the minimum score (Q0) in every output, including the not-merged output files. Output quality scores can be clamped with --fastq_qmaxout and --fastq_qminout (these apply only to the merged region). Unlike the commands that pass an input quality through, --fastq_mergepairs keeps the pre-2.32.0 --fastq_qmaxout default of 41, because the score it clamps is computed rather than read: two agreeing Q40 bases have a posterior quality of Q85, and reporting it would change the merged output of every run. Pass --fastq_qmaxout 93 for the unclamped posterior. The merged scores are written with the input offset, --fastq_ascii: this is the only command whose output encoding tracks its input encoding, and it does not accept --fastq_asciiout.

Staggered pairs — where the 3’ end of the reverse read extends past the 5’ end of the forward read — are discarded by default. Use --fastq_allowmergestagger to allow them; the overhanging portions of both reads are excluded from the merged sequence, which then covers the overlap only.

Reads can be pre-filtered with --fastq_truncqual, --fastq_minlen, --fastq_maxlen (the length bounds apply after truncation), and --fastq_maxns; the expected-error filter --fastq_maxee applies to the merged sequence. Bounds on the merged sequence length are set with --fastq_minmergelen and --fastq_maxmergelen.

To illustrate a merge with a 6-base overlap:

Forward (5'→3'):   AAAATTTTTT
Reverse (5'→3'):   GGGGAAAAAA

Aligned:
  Forward:      AAAATTTTTT
  Rev-comp:         TTTTTTCCCC   (reverse complement of the reverse read)

Overlap:            TTTTTT      (6 bases; scores combined)
Merged:         AAAATTTTTTCCCC

At the end of the run, vsearch prints a report — to standard error (unless --quiet is given), or to the log file when --log is used — giving the number of merged pairs and, for the pairs that could not be merged, a breakdown by reason. Most reasons correspond directly to a user-adjustable threshold: reads too short (after truncation) (--fastq_minlen), reads too long (after truncation) (--fastq_maxlen), too many N's (--fastq_maxns), too many differences (--fastq_maxdiffs), too high percentage of differences (--fastq_maxdiffpct), overlap too short (--fastq_minovlen), expected error too high (--fastq_maxee), merged fragment too short (--fastq_minmergelen), merged fragment too long (--fastq_maxmergelen), and staggered read pairs (allowed with --fastq_allowmergestagger). The remaining reasons reflect the alignment heuristics and are not directly controllable:

too few kmers found on same diagonal: Too few k-mers shared by the two reads fall on a common alignment diagonal, so no candidate overlap could be located and no alignment was attempted.
multiple potential alignments: More than one overlap of comparable quality was found — for instance when the overlap region contains a tandem repeat — so the correct alignment is ambiguous and the pair is left unmerged.
alignment score too low, or score drop too high: An overlap was found and aligned, but its score remained below the internal threshold, or the score dropped too sharply within the overlap. This usually points to clustered mismatches or an indel in the overlap region.

ALGORITHM

Merging happens in three stages: candidate overlaps are located with shared words, each candidate is scored by an ungapped alignment whose scores are derived from the base qualities, and the best candidate is accepted only if it is both good enough and unambiguous. The thresholds below are internal and have no options of their own.

Locating candidates. The forward read is indexed by its 5-mers (a shorter word than the search commands use), and every possible overlap length is rated by the number of 5-mers the two reads share at that offset. Only offsets sharing at least four 5-mers are examined further. A pair where no offset reaches that count is discarded as too few kmers found on same diagonal, without any alignment being attempted. When --fastq_minovlen is below 9 the requirement is relaxed to --fastq_minovlen minus 4.

Scoring a candidate. Each surviving offset is aligned without gaps, the forward read against the reverse complement of the reverse read, and scored in bits as a log-odds ratio. Writing p for the probability that two truly identical bases are observed as a match, given the error probabilities p_f and p_r of the forward and reverse base:

p = 1 - p_f - p_r + 4 x p_f x p_r / 3

a matching column scores log2(p / 0.25) and a mismatching column log2((1 - p) / 0.75), the latter capped at -4 bits so that one disagreement between two poor bases cannot dominate. Two good bases that agree are therefore worth almost 2 bits each, while a disagreement between two good bases costs many. The formulas are those of Edgar & Flyvbjerg (2015), which also give the posterior qualities written for the merged region.

Accepting a candidate. The running score and its maximum so far are tracked as the overlap is walked, and a candidate whose score falls more than 16 bits below its own maximum is dropped: that is how a run of clustered mismatches is rejected even when the two ends align well. A candidate is acceptable once its final score reaches 16 bits, or 1.6 times --fastq_minovlen when that is below 9 — about nine matching bases of good quality, since eight Q40 matches come to 15.998 bits. If more than one candidate is acceptable the overlap is ambiguous and the pair is discarded as multiple potential alignments, which is what a tandem repeat in the overlap produces. Otherwise the highest-scoring candidate is kept, and a pair whose best candidate never reached the threshold is discarded as alignment score too low, or score drop too high.

Because the alignment is ungapped, an indel in the overlap is not modelled as such: it shows up as a run of mismatches and is rejected by the score-drop rule.

OPTIONS

mandatory options

--reverse filename: Specify the FASTQ file containing the reverse reads.

At least one output option is required: one of the merged-read outputs below, one of the not-merged outputs (--fastaout_notmerged_fwd, --fastaout_notmerged_rev, --fastqout_notmerged_fwd, --fastqout_notmerged_rev), or --eetabbedout (see the following sections).

--fastaout filename: Write merged sequences to filename, in fasta format.
--fastqout filename: Write merged sequences to filename, in fastq format (see vsearch-fastq(5)).

core options

--fastq_ascii 33|64: Specify the offset used as the basis for the fastq quality score when reading fastq files. For example, an offset of 33 means that a quality value of 41 is represented by the 74th ASCII symbol (33 + 41 = 74), which is ‘J’. See ascii(7) for a view of the ASCII character set. The offset value is either 33 or 64, default is 33.

The offset matters even to a command that never decodes a quality score: it is what the reader compares the observed quality symbols against before warning that the file may use the other encoding, it sets the default of --fastq_qmax (the highest score the offset can represent), and the sum rules on --fastq_qmin and --fastq_qmax are stated in terms of it.

--fastq_minovlen positive integer: Set the minimum length of the overlap region between the forward and reverse reads. Must be at least 5. The default is 10.
--fastq_maxdiffs positive integer: Set the maximum number of mismatches allowed in the overlap region. This option has a strong influence on the merging success rate. The default is 10.
--fastq_maxdiffpct real: Set the maximum percentage of mismatches allowed in the overlap region (0.0 to 100.0). Additional heuristics in the merging algorithm may still discard pairs with a high mismatch rate. The default is 100.0. usearch bounds the same quantity from the other side, as a minimum identity: --fastq_maxdiffpct x is equivalent to usearch’s --fastq_pctid 100 - x.
--fastq_allowmergestagger: Allow merging of staggered read pairs. Staggered pairs arise when a very short fragment is sequenced: the 3’ end of the reverse read extends beyond the 5’ end of the forward read. The overhanging portion of the reverse read is discarded. By default, staggered pairs are not merged (see --fastq_nostagger).
--fastq_nostagger: Discard staggered read pairs. This is the default behaviour. See --fastq_allowmergestagger to change it.
--fastq_minmergelen positive integer: Discard merged sequences shorter than positive integer bases. The default is 0 (no minimum).
--fastq_maxmergelen positive integer: Discard merged sequences longer than positive integer bases. The default is 1,000,000.
--fastq_qmax integer: Specify the maximal quality score accepted when reading fastq sequences. Stop with an error message if a quality score higher than the specified value is read. The offset (see --fastq_ascii) plus the maximal score may not exceed 126, the last printable ASCII character: scores up to 93 with offset 33, up to 62 with offset 64. The default is the highest score the offset can represent (93 with offset 33, 62 with offset 64), so no quality score is rejected unless this option is lowered. Before version 2.32.0 the default was 41, the usual maximum for Sanger/Illumina 1.8+ files, which rejected PacBio HiFi and nanopore files outright.
--fastq_qmin integer: Specify the minimal quality score accepted when reading fastq sequences. Stop with an error message if a quality score lower than the specified value is read. The offset (see --fastq_ascii) plus the minimal score must be at least 33, the first printable ASCII character: scores down to 0 with offset 33, down to -31 with offset 64. The value may therefore be negative, which is what the negative scores of older formats require, but note that those formats are not supported (see vsearch-fastq(5)). The default is 0, which is usual for recent Sanger/Illumina 1.8+ files.
--fastq_qmaxout integer: Specify the maximum quality score used when writing fastq files. The default is the highest score the output offset can represent (93 with --fastq_asciiout 33, 62 with offset 64), so scores read from the input are written back unchanged. Before version 2.32.0 the default was 41, the usual maximum for Sanger/Illumina 1.8+ files, which silently reduced any higher score. Older formats may use a maximum quality score of 40. Two commands are exceptions and keep the old default of 41, because they generate the score they clamp instead of passing one through: --fasta2fastq, which has no input quality and uses this option as the value to write, and --fastq_mergepairs, which caps the computed posterior quality of a merged base.

For --fastq_mergepairs the offset in question is --fastq_ascii, not --fastq_asciiout: that command writes fastq but does not accept --fastq_asciiout, so a merged quality symbol carries the same offset the input was read with. The sum rule is stated against --fastq_ascii there, and against --fastq_asciiout everywhere else.

--fastq_qminout integer: Specify the minimum quality score used when writing fastq files. The offset (see --fastq_asciiout) plus the minimum score must be at least 33, so the value may be negative when the output offset is 64. The default is 0, which is usual for recent Sanger/Illumina 1.8+ files. Older formats may use scores between -5 and 2.

A negative minimum matters where the score passes through vsearch and may itself be negative: --fastq_convert on an offset-64 file read with a lowered --fastq_qmin writes the score back unchanged with --fastq_qminout -5, where the default 0 would raise it. It has no effect where vsearch computes the score it clamps, because such a score is derived from an error probability and is therefore never negative: --fastq_mergepairs, and --fastq_convert --fastq_solexa, which converts to the Phred scale before any output clamp applies.

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.
--eeout: Add the expected error count to each sequence header in output fasta or fastq files, as the annotation ;ee=float. Synonym of --fastq_eeout. Use --xee to remove this annotation from headers. See vsearch-expected_error(7).
--eetabbedout filename: Write expected error statistics for each successfully merged pair to filename (pairs that could not be merged are not reported), in a tab-separated format with four columns: the expected errors in the forward read, the expected errors in the reverse read, the observed differences in the forward read within the overlap region, and the observed differences in the reverse read within the overlap region. See vsearch-expected_error(7).
--fasta_width positive integer: Set the maximal width of sequences when writing fasta files. Longer sequences are folded and written on several lines. Default width is 80 nucleotides. Set to zero (0) to suppress folding.
--fastaout_notmerged_fwd filename: Write forward reads that could not be merged to filename, in fasta format.
--fastaout_notmerged_rev filename: Write reverse reads that could not be merged to filename, in fasta format.
--fastq_eeout: Add the expected error count to each sequence header in output fasta or fastq files, as the annotation ;ee=float. Synonym of --eeout. Use --xee to remove this annotation from headers. See vsearch-expected_error(7).
--fastq_maxee positive real: Discard sequences with an expected error greater than real. The expected error is the sum of error probabilities for all positions in the sequence, and is strictly positive (zero or negative arguments are rejected, as they would discard every sequence). Applied after trimming. See vsearch-expected_error(7).
--fastq_maxlen positive integer: Discard sequences longer than positive integer bases. Applied after trimming.
--fastq_maxns positive integer: Discard sequences containing more than positive integer ambiguous bases (N).
--fastq_minlen positive integer: Discard sequences shorter than positive integer bases. Applied after trimming. Default is 1.
--fastq_truncqual positive integer: Truncate reads starting at the first base whose quality score is at or below positive integer. Accepted values range from 0 to 93.
--fastqout_notmerged_fwd filename: Write forward reads that could not be merged to filename, in fastq format (see vsearch-fastq(5)).
--fastqout_notmerged_rev filename: Write reverse reads that could not be merged to filename, in fastq format (see vsearch-fastq(5)).
--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.
--label_suffix string: Add the suffix string to sequence headers when writing fasta or fastq files. For example, with --label_suffix ";status=healthy", sequence header ‘>seq1’ becomes ‘>seq1;status=healthy’.
--lengthout: Add a sequence length annotation (;length=integer) to each sequence header when writing fasta or fastq files.
--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).
--quiet: Suppress messages to the standard output stdout(3) and standard error stderr(3), except for warnings and error messages.
--relabel string: Replace sequence headers with the prefix string and a ticker (1, 2, 3, etc.). For example, with --relabel "cluster:", the first sequence header becomes ‘>cluster:1’, the second sequence header becomes ‘>cluster:2’, and so on. To retain annotations, use their corresponding options (--lengthout, --eeout, and --sizeout). Use --relabel_keep to also retain old sequence identifiers.
--relabel_keep: Retain old sequence identifiers by including them at the end of the new headers, after a space.
--relabel_md5: Replace each sequence header with the MD5 digest derived from the sequence itself. The sequence is converted to upper case, and each ‘U’ is replaced with a ‘T’ before computation of the digest. The MD5 digest is a 128-bit value (16 bytes), represented using a string of 32 ASCII characters. Each pair of characters encodes an hexadecimal value, ranging from x00 to xff. See md5(3) for more details, and --relabel_sha1 for an alternative hashing algorithm. To retain annotations, use their corresponding options (--lengthout, --eeout, and --sizeout). Use --relabel_keep to also retain old sequence identifiers.
--relabel_self: Replace each sequence header with the sequence itself. To retain annotations, use their corresponding options (--lengthout, --eeout, and --sizeout). Use --relabel_keep to also retain old sequence identifiers.
--relabel_sha1: Replace each sequence header with the SHA1 digest derived from the sequence itself. The sequence is converted to upper case, and each ‘U’ is replaced with a ‘T’ before computation of the digest. The SHA1 digest is a 160-bit value (20 bytes), represented using a string of 40 ASCII characters. Each pair of characters encodes an hexadecimal value, ranging from x00 to xff. See sha1(3) for more details, and --relabel_md5 for an alternative hashing algorithm. To retain annotations, use their corresponding options (--lengthout, --eeout, and --sizeout). Use --relabel_keep to also retain old sequence identifiers.
--sample string: Add the given sample identifier string to sequence headers when writing fasta or fastq files. For instance, if string is ‘ABC’, the text ;sample=ABC will be added to the headers. string is silently truncated at the first ‘;’ or whitespace character (space, tab, newline, carriage return, vertical tab or form feed), so such characters should not be used in string. Other characters (alphabetical, numerical and punctuations) are accepted. When nothing is left after truncation — an empty string, or one starting with ‘;’ or a blank character — vsearch issues a warning and writes a bare ;sample= annotation.
--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.
--sizeout: Add abundance annotations to sequence headers when writing fasta or fastq files. Add the pattern ;size=integer. Existing ;size= annotations are reported unchanged; entries without one receive ;size=1. For this command --sizein is not needed: abundance annotations are always parsed from the input headers.
--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.
--xee: Strip expected error (ee) annotations from sequence headers when writing fasta or fastq files. Search for the pattern [>@;]ee=float[;]. Expected error annotations are added by the synonymous options --fastq_eeout and --eeout described in vsearch-fastx_filter(1).
--xlength: Strip sequence length annotations from sequence headers when writing fasta or fastq files. Search for the pattern [>@;]length=integer[;]. Sequence length annotations are added by the --lengthout option.
--xsize: Strip abundance annotations from sequence headers when writing fasta or fastq files. Search for the pattern [>@;]size=integer[;]. Abundance annotations are added by the --sizeout option.

EXAMPLES

Merge paired-end reads and write merged sequences to a fastq file:

vsearch \
    --fastq_mergepairs fwd.fastq \
    --reverse rev.fastq \
    --fastqout merged.fastq

Merge with a stricter overlap and mismatch threshold, and save unmerged reads for inspection:

vsearch \
    --fastq_mergepairs fwd.fastq \
    --reverse rev.fastq \
    --fastq_minovlen 20 \
    --fastq_maxdiffs 5 \
    --fastqout merged.fastq \
    --fastqout_notmerged_fwd unmerged_fwd.fastq \
    --fastqout_notmerged_rev unmerged_rev.fastq

Allow staggered pairs and filter on expected error after merging:

vsearch \
    --fastq_mergepairs fwd.fastq \
    --reverse rev.fastq \
    --fastq_allowmergestagger \
    --fastq_maxee 1.0 \
    --fastqout merged.fastq

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/.

COPYRIGHT

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:

Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
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.