samtools-consensus(1)

📦
version 1.24 · 9 versions
, in epel 10, nixpkgs unstable

✨ NAME ✨

okay so bestie, welcome to samtools-consensus(1) — the manual for samtools consensus, aka the tool that produces a consensus FASTA/FASTQ/PILEUP. it’s giving reference-building main character energy, no cap 💅

🌈 SYNOPSIS 🌈

samtools consensus [-saAMq] [-r region] [-f format] [-l line-len] [-d min-depth] [-C cutoff] [-c call-fract] [-H het-fract] in.bam

🦋 DESCRIPTION 🦋

so the whole point: generate a consensus from a SAM, BAM, or CRAM file based on what’s actually inside the alignment records. the consensus gets written as either FASTA, FASTQ, or a pileup-oriented format — and you pick which one with the -f FORMAT option. slay ✨

the default output for FASTA and FASTQ formats includes one base per non-gap consensus. so insertions with respect to the aligned reference will be included, and deletions get removed. this behaviour can be controlled with the --show-ins and --show-del options. this could be used to compute a new reference from sequences assemblies to realign against — iconic little trick 🫶

the pileup-style format strictly adheres to one row per consensus location, differing from the one row per reference base used in the related samtools mpileup command. this means the base quality values for inserted columns are reported. the base quality value of gaps — either within an insertion or otherwise — are determined as the average of the surrounding non-gap bases. the columns shown are: reference name, position, nth base at that position (zero if not an insertion), consensus call, consensus confidence, sequences and quality values. note that even when a reference is supplied, the consensus base is always reported (if non-zero depth) in the 5th column. okay periodt, column five stays booked and busy 💖

two consensus calling algorithms are offered. the default computes a heterozygous consensus in a Bayesian manner, derived from the "Gap5" consensus algorithm. quality values are also tweaked to take into account other nearby low quality values. this can also be disabled with the --no-adj-qual option. lowkey thoughtful of it 🎀

this method also utilises the mapping qualities, unless the --no-use-MQ option is used. mapping qualities are also auto-scaled to take into account the local reference variation by processing the MD:Z tag, unless --no-adj-MQ is used. mapping qualities can be capped between a minimum (--low-MQ) and maximum (--high-MQ), although the defaults are liberal and trust the data to be true. finally, an overall scale on the resulting mapping quality can be supplied (--scale-MQ, defaulting to 1.0). this has the effect of favouring more calls with a higher false positive rate (values greater than 1.0) or being more cautious with higher false negative rates and lower false positive (values less than 1.0). choose your fighter, bestie ⭐

the second method is a simple frequency counting algorithm, summing either +1 for each base type or +qual if the --use-qual option is specified. this is enabled with the --mode simple option. cute and straightforward 🌸

the summed share of a specific base type is then compared against the total possible, and if this is above the --call-fract fraction parameter then the most likely base type is called, or "N" otherwise (or absent if it is a gap). the --ambig option permits generation of ambiguity codes instead of "N", provided the minimum fraction of the second most common base type to the most common is above the --het-fract fraction. iykyk 🧚

💅 OPTIONS 💅

🌟 General options that apply to both algorithms 🌟

  • -r REG, --region REG — limit the query to region REG. this requires an index. it is incompatible with the --regions-file option. no sharing, bestie 🙅
  • --regions-file BED-FILE — creates consensus sequences for each region listed in the BED-FILE, although not necessarily in the same chromosomal order. as with single region queries, without -a the sequences may be truncated or absent if not covered. note this does not deduplicate regions or merge overlapping regions. this requires an index. this is incompatible with the -r option. messy but valid 🫶
  • -f FMT, --format FMT — produce format FMT, with "fastq", "fasta" and "pileup" as permitted options.
  • -l N, --line-len N — sets the maximum line length of line-wrapped fasta and fastq formats to N.
  • -o FILE, --output FILE — output consensus to FILE instead of stdout.
  • -m STR, --mode STR — select the consensus algorithm. valid modes are "simple" frequency counting and the "bayesian" (Gap5) methods, with Bayesian being the default. (note case does not matter, so "Bayesian" is accepted too.) there are a variety of Bayesian methods. straight "bayesian" is the best set suitable for the other parameters selected. the choice of internal parameters may change depending on the "--P-indel" score. this method distinguishes between substitution and indel error rates. the old Samtools consensus in version 1.16 did not distinguish types of errors, but for compatibility the "bayesian_116" mode may be selected to replicate this. choosing a mode is literally a personality test 💀
  • -a — outputs all bases, from start to end of reference, even when the aligned data does not extend to the ends. this is most useful for construction of a full length reference sequence. full-length queen behavior 👑
  • -a -a, -aa — output absolutely all positions, including references with no data aligned against them. absolutely unbothered by empty coverage ✨
  • --rf, --incl-flags STR|INT — only include reads with at least one FLAG bit set. defaults to zero, which filters no reads.
  • --ff, --excl-flags STR|INT — exclude reads with any FLAG bit set. defaults to "UNMAP,SECONDARY,QCFAIL,DUP".
  • --min-MQ INT — filters out reads with a mapping quality below INT. this defaults to zero.
  • --min-BQ INT — filters out bases with a base quality below INT. this defaults to zero.
  • --show-del yes/no — whether to show deletions as "*" (yes) or to omit them from the output (no). defaults to no.
  • --show-ins yes/no — whether to show insertions in the consensus. defaults to yes.
  • --mark-ins — insertions, when shown, are normally recorded in the consensus with plain 7-bit ASCII (ACGT, or acgt if heterozygous). however, this makes it impossible to identify the mapping between consensus coordinates and the original reference coordinates. if fasta output is selected then the option adds an underscore before every inserted base, plus a corresponding character in the quality for fastq format. when used in conjunction with -a --show-del yes, this permits an easy derivation of the consensus to reference coordinate mapping. coordinate girlies rise 🎀
  • -A, --ambig — enables IUPAC ambiguity codes in the consensus output. without this the output will be limited to A, C, G, T, N and *.
  • -d D, --min-depth D — the minimum depth required to make a call. defaults to 1. failing this depth check will produce consensus "N", or absent if it is an insertion. note this check is performed after filtering by flags and mapping/base quality. order matters, bestie 🔥
  • -T ref.fa, --reference ref.fa — for base positions with zero coverage, use the supplied reference instead of "N". note this does not replace minimum depth or minimum quality filters, as the base is known but considered low quality, so the ambiguity is retained. she knows the tea but keeps the question mark 🫶
  • --ref-qual INT — when --reference is given, this specifies the quality value to use for reference-derived bases. this defaults to zero.

🌸 The following options apply only to the simple consensus mode 🌸

  • -q, --use-qual — for the simple consensus algorithm, this enables use of base quality values. instead of summing 1 per base called, it sums the base quality instead. these sums are also used in the --call-fract and --het-fract parameters too. quality values are always used for the "Gap5" consensus method and this option has no effect. note currently quality values only affect SNPs and not inserted sequences, which still get scores with a fixed +1 per base type occurrence. quality is having a moment, but only for SNPs 💅
  • -H H, --het-fract H — for consensus columns containing multiple base types, if the second most frequent type is at least H fraction of the most common type, then a heterozygous base type will be reported in the consensus. otherwise the most common base is used, provided it meets the --call-fract parameter (otherwise "N"). the fractions computed may be modified by the use of quality values if the -q option is enabled. note although IUPAC has ambiguity codes for A,C,G,T vs any other A,C,G,T, it does not have codes for A,C,G,T vs gap (such as in a heterozygous deletion). given the lack of any official code, we use lower-case letter to symbolise a half-present base type. lowercase = half-present, iconic and slightly unhinged 🌈
  • -c C, --call-fract C — only used for the simple consensus algorithm. require at least C fraction of bases agreeing with the most likely consensus call to emit that base type. this defaults to 0.75. failing this check will output "N".
  • -@ NTHREADS — specify the number of additional threads to use for computing the consensus. note if no index is present, threads will only be used for parallel decompression, meaning asking for more than 2 threads is unlikely to speed up processing. with an index, the consensus is computed for multiple regions simultaneously, offering near linear speed ups. multithreading girlies, assemble ⚡
  • -Z BASE_COUNT — when using multiple threads, this specifies the number of bases per threading job. the default is 500,000 bp for fasta/fastq output and 100,000 for pileup output. larger blocks may yield improved threading performance at a cost of more memory.

🔥 The following options apply only to Bayesian consensus mode enabled (default on) 🔥

  • -C C, --cutoff C — only used for the Gap5 consensus mode, which produces a Phred style score for the final consensus quality. if this is below C then the consensus is called as "N". the cutoff is giving strict gatekeeper energy 🚪
  • --use-MQ, --no-use-MQ — enable or disable the use of mapping qualities. defaults to on.
  • --adj-MQ, --no-adj-MQ — if mapping qualities are used, this controls whether they are scaled by the local number of mismatches to the reference. the reference is unknown by this tool, so this data is obtained from the MD:Z auxiliary tag (or ignored if not present). defaults to on. MD:Z stays doing the most, fr fr 💜
  • --NM-halo INT — how far to either side of the base call we peek when counting local mismatches ✨ like, it sets the radius of the vibe check around each base 💅
  • --low-MQ MIN, --high-MQ MAX — the floor and the ceiling for mapping quality. ⚠️ Heads up bestie: these are not filters, they just cap the values from above and below. Defaults are 0 and 60 🎀
  • --scale-MQ FLOAT — a general multiplicative mapping-quality scaling factor 🦋 it globally raises or lowers the quality values that the consensus algorithm actually uses. Defaults to 1.0, which is the "leave my values alone" setting 😌
  • --P-het FLOAT — controls how likely any given position is to be a heterozygous site 👀 this feeds the priors for the Bayesian calculations, and it makes barely any difference on deep data. Defaults to 1e-3. Smaller numbers make the algorithm more likely to call a pure base type. NGL the algorithm always computes the probability of a base being homozygous vs heterozygous, no matter whether the output gets reported as ambiguous (it'll come out as "N" if it's deemed heterozygous and you didn't turn on --ambig mode) 💀
  • --P-indel FLOAT — controls how likely small indels are 🔥 again, priors for the Bayesian maths, and again it barely matters on deep data. Defaults to 2e-4 ✨
  • --het-scale FLOAT — a multiplicative correction applied per base quality before it's added to the heterozygous hypotheses. Turning it down = fewer heterozygous calls. This often leads to a significant reduction in false positive het calls, at the cost of some increase in false negatives (aka mislabelling real heterozygous sites as homozygous 😭). It's usually a good idea to reduce this on instruments where a significant proportion of bases may be aligned in the wrong column due to insertions and deletions, which leads to alignment errors and reference bias. Think of it as a het sensitivity tuning knob 🎛️ Defaults to 1.0 (nop).
  • -p, --homopoly-fix — some technologies that call runs of the same base type together always stick the lowest quality calls at one end 🙈 which causes drama when you reverse complement and compare alignments with indels. This option averages the qualities at both ends so orientation biases don't ruin your day. Recommended for old 454 or PacBio HiFi data sets 👑
  • --homopoly-score FLOAT — the -p option also drops confidence values inside homopolymers, because there's an extra likelihood of sequence-specific errors. The quality values get multiplied by FLOAT. Defaults to 0.5, but it's not used if -p wasn't specified. Lowkey convenient: adjusting this score automatically enables -p ✨
  • -t, --qual-calibration FILE — loads a quality calibration table from FILE 📂 the format is a series of lines with these fields, each line starting with the literal text QUAL:
QUAL value substitution undercall overcall

Lines starting with a # are ignored 🙅‍♀️ each line maps a recorded quality value to the Phred equivalent score for substitution, undercall and overcall errors. Quality values are expected to be sorted in increasing numerical order, but they're allowed to skip values. This lets the consensus algorithm know the most likely cause of an error, and whether the instrument is more likely to have indel errors (more common in some long-read technologies) or substitution errors (more common in clocked short-read instruments) 🧚

Some pre-defined calibration tables are built in 💖 these are specified with a fake filename starting with a colon. See the -X option for more details.

Note: due to the additional heuristics the consensus algorithm applies, these recalibration tables are not a true reflection of the instrument error rates — they're a work in progress 🚧
  • -X, --config STR — picks predefined sets of configuration parameters 🎀 acceptable values for STR are below, along with the exact list of parameters they're equivalent to:
    • hiseq — --qual-calibration :hiseq
    • hifi — --qual-calibration :hifi --homopoly-fix 0.3 --low-MQ 5 --scale-MQ 1.5 --het-scale 0.37
    • r10.4_sup — --qual-calibration :r10.4_sup --homopoly-fix 0.3 --low-MQ 5 --scale-MQ 1.5 --het-scale 0.37
    • r10.4_dup — --qual-calibration :r10.4_dup --homopoly-fix 0.3 --low-MQ 5 --scale-MQ 1.5 --het-scale 0.37
    • ultima — --qual-calibration :ultima --homopoly-fix 0.3 --low-MQ 10 --scale-MQ 2 --het-scale 0.37

🌈 EXAMPLES 🌈

  • Make a modified FASTA reference that has a 1:1 coordinate correspondence with the original reference used in alignment 👑 like, every coordinate lines up perfectly, no cap:
samtools consensus -a --show-ins no --show-del yes in.bam -o ref.fa
  • Make a FASTQ file for the contigs with aligned data, insertions included 🦋:
samtools consensus -f fastq in.bam -o cons.fq

💜 AUTHOR 💜

Written by James Bonfield from the Sanger Institute 🫶 absolute icon, we stan.

⭐ SEE ALSO ⭐

samtools(1), samtools-mpileup(1),

Samtools website: <http://www.htslib.org/> 🌸

okay that's the whole page bestie, go forth and call some consensus bases 💅✨ this one's samtools-1.24, dated 9 July 2026, and it's the samtools-consensus(1) manual page — periodt 🎀

📦 where this page comes from