Iterated consensus sequence

NOTE This package is only a day old, so you should expect occasional breaking changes!

Iteratively build a consensus sequence: call a consensus from a BAM (or an initial reference), build a mapper index from it, remap the original reads against it, call a new consensus, and repeat until the consensus stops changing.

Mapping and consensus-calling are both fully user-configurable via a TOML pipeline config -- you bring your own bowtie2/bwa/minimap2/... and ivar/samtools/... commands, iterated-consensus just drives the loop, tracks convergence, and writes the results.

The code here was written entirely by Claude Sonnet 5 (model ID claude-sonnet-5), from Anthropic's Claude 5 family. This took about 3.5 hours from giving Claude an initial description, through planning and writing code, tests, and documentation.

Install

uv add iterated-consensus     # or: pip install iterated-consensus

Quick start

iterated-consensus config-template                  # list bundled presets
iterated-consensus config-template bowtie2-ivar > pipelines.toml
# edit pipelines.toml: fill in [input], adjust commands/threads as needed
iterated-consensus run --config pipelines.toml --out-dir results/ --dry-run  # preview
iterated-consensus run --config pipelines.toml --out-dir results/ --progress

--progress prints a one-line summary after each iteration (reads mapped, consensus length, identity to the previous consensus, time taken). Every run also writes results/index.html -- open it in a browser for a summary and full per-iteration detail, no --progress needed.

Errors (bad config, a failing command, a mismatched reference, ...) normally print a short error: ... message. Add --traceback to instead let them crash with the full Python traceback, for debugging.

Config format

A config has one or more [[mapper]] tables, one [consensus] table, and optionally [input], [output], and [run].

[[mapper]]
name = "bowtie2"
index_cmd = ["bowtie2-build", "{reference}", "{index_prefix}"]
map_cmd = "bowtie2 -x {index_prefix} -1 {reads_1:,} -2 {reads_2:,} -p {threads} | samtools sort -o {bam}"

# A second [[mapper]] table can be added too -- if more than one mapper is
# configured, every mapper runs each iteration and their BAMs are merged
# before the consensus step sees them.

[consensus]
steps = [
    "samtools mpileup -aa -A -d 0 -Q 0 -f {reference} {bam} | ivar consensus -p {consensus_prefix} -t 0.5",
]
output = "{consensus_prefix}.fa"   # where to find the result -- see note below

[input]
reads_1 = ["a_R1.fastq.gz", "b_R1.fastq.gz"]  # 0 or more paired sets
reads_2 = ["a_R2.fastq.gz", "b_R2.fastq.gz"]
reads_single = []                             # 0 or more single-end files
reference_fasta = "starting_reference.fasta"  # a local file...
# reference_id = "chr2"        # ...and/or a name -- see "Reference resolution" below
# --- or, instead of the FASTQ block above, start from a BAM: ---
# bam = "input.bam"
# reference_id = "chr2"        # only needed if the BAM has >1 reference
# reference_fasta = "chr2.fasta"  # optional -- see "Reference resolution" below
# bam_reads = "ref"            # ref | ref+unal | all -- see below

[output]
# consensus_fasta = "final_consensus.fasta"   # copy the last iteration's
#   consensus here once the run finishes, converged or not -- see below
# consensus_id = "my-sample-name"              # optional new FASTA header

[run]
threads = 4                    # or "auto" -- see "Threads and custom [run]
                                # variables" below
max_iterations = 20
convergence_identity = 100.0   # stop once consensus identity to the previous
                                # iteration reaches this percent...
convergence_streak = 1         # ...for this many iterations in a row

# Anything else here becomes a {name} placeholder in every command, e.g.:
# min_depth = 10               # -> {min_depth} in [consensus] steps
# sample_name = "patient-42"   # -> {sample_name} anywhere

[input] can instead (or partly) be supplied on the command line -- see iterated-consensus run --help. CLI values override the config's [input] field-by-field, so a config can be fully self-contained or left generic and pointed at different data per invocation.

[consensus].output is only used right after the steps run, to find and read whatever file your tool actually wrote (different tools name it differently -- ivar consensus -p PREFIX writes PREFIX.fa, which is why the example above is output = "{consensus_prefix}.fa", not {consensus_prefix} alone). What gets read there is then copied to this iteration's own consensus.fasta (see "Output" below) -- and it's that fixed-name copy, not the path output pointed to, that later becomes {reference} for the next iteration. This is why a mapper's index_cmd in --dry-run output references consensus.fasta even if your output pattern produces a different filename or extension: output only has to match what your consensus tool actually writes, nothing downstream reads that path directly.

Command steps: list or shell string

Every command (index_cmd, map_cmd, each [consensus] step) can be written as a list of argv tokens (run directly, no shell -- safest, use this whenever you're just running one program) or as a single string (run via a shell -- needed for pipes, as in the samtools mpileup | ivar consensus example above).

After each mapper's map_cmd runs (and after merging, if more than one mapper is configured), iterated-consensus checks for a .bam.bai index next to the resulting BAM and creates one with samtools index if it's missing -- most consensus tools need one, so you don't have to remember to add an indexing step to map_cmd yourself. Indexing needs the BAM coordinate-sorted first, so the BAM header's own SO tag is checked before doing anything else: if it's already marked SO:coordinate, indexing runs directly (fast, and avoids a pointless sort pass over a large, already-sorted BAM); if not (map_cmd forgot a samtools sort, or a mapper wrote it unsorted), it's sorted in place automatically first, so map_cmd doesn't strictly need its own sort step either. (A BAM whose header lies about being sorted -- rare, but not impossible -- is still caught: if indexing the "already-sorted" file fails anyway, it's sorted for real and indexing is retried.) This only covers the BAM a mapping step actually produces; it doesn't apply to iter_000 of a BAM-start run, which has no mapping step (see "Reference resolution" above) -- add your own samtools index {bam} step to [consensus] if that BAM needs indexing too (the bundled bwa-samtools preset does exactly this).

Placeholders

• {reference} -- the current reference FASTA: the previous iteration's consensus, or the starting reference for iteration 0. For a BAM-start run, iteration 0's reference is only available if it could be resolved -- see "Reference resolution" below; if not, and [consensus] uses {reference} anyway, that's a config error caught before anything runs, not a silent guess.
• {index_prefix} -- path prefix for this mapper's index this iteration.
• {bam} -- path this mapper should write its BAM to.
• {consensus_prefix} -- path prefix for the consensus step's output (see the note on [consensus].output above -- this is not the same file that later becomes {reference}).
• {threads} -- from [run] threads. See "Threads and custom [run] variables" below for threads = "auto" and defining your own placeholders alongside it (e.g. for splitting a thread budget across a piped command).
• {reads_1}, {reads_2}, {reads_single} -- the read-file lists, matching [input]'s reads_1/reads_2/reads_single one-for-one. Only present if that category is non-empty for this run -- referencing e.g. {reads_single} in a run with no unpaired reads is a config error, so a mapper template should only reference the categories it actually expects.

Threads and custom [run] variables

Every [run] key is available as a {name} placeholder in every index_cmd, map_cmd, and [consensus] step -- both the ones with dedicated meaning ({threads}, {max_iterations}, {convergence_identity}, {convergence_streak}, {threads_reserve}) and any custom ones you add. This is handy for more than just thread counts -- e.g. a logging step that records what a run was configured with:

[run]
threads = 8
min_depth = 10
sample_name = "patient-42"

[consensus]
steps = [
    "samtools mpileup -d 0 {bam} | ivar consensus -t 0.5 -m {min_depth} -p {consensus_prefix}",
    'echo "Built {sample_name} consensus at {threads} threads, target {convergence_identity}% over {max_iterations} iterations" >> log.txt',
]

A custom variable can't reuse a name iterated-consensus already sets itself (the dedicated [run] fields above, plus the per-iteration placeholders reference, index_prefix, bam, consensus_prefix, reads_1, reads_2, reads_single) -- that's rejected at config load time rather than silently shadowed.

threads = "auto" resolves to the number of CPUs actually available to the process (respecting container/cgroup/taskset limits on Linux, where that's exposed; the installed core count elsewhere) at config-load time, once per run -- not re-detected per iteration. Pair it with threads_reserve (an integer, only valid alongside threads = "auto") to leave some cores free for other work on the machine:

[run]
threads = "auto"
threads_reserve = 2   # use (detected CPUs - 2), never less than 1

Splitting a thread budget across a pipe. {threads} is one number, but a piped command like bwa mem | samtools sort or samtools mpileup | ivar consensus runs two programs concurrently, each of which could use its own thread count -- and since they're running at the same time, those counts add up against your actual core count, they don't each get to use the whole budget. Using {threads} unmodified for more than one stage of the same pipe oversubscribes the machine. iterated-consensus doesn't try to auto-split {threads} for you -- pipeline stages have wildly different threading characteristics (some don't support it at all, some scale linearly, some plateau early), so a generic split would just be a guess. Instead, treat {threads} (or a threads = "auto" budget) as the total, and partition it yourself into named [run] variables that add up to no more than that:

[run]
threads = "auto"
threads_reserve = 2   # e.g. resolves to 6 on an 8-core machine
map_threads = 5       # give most of the budget to the mapper...
sort_threads = 1       # ...and a thread or two to samtools sort running
                        # alongside it in the same pipe

[[mapper]]
name = "bwa"
map_cmd = "bwa mem -t {map_threads} {index_prefix} {cat:reads_1} {cat:reads_2} | samtools sort -@ {sort_threads} -o {bam}"

{threads} is still the right placeholder for any command that's just one program (most index_cmds, or a [consensus] step with no pipe) -- reach for named variables like map_threads/sort_threads only where a pipe means two or more programs are genuinely running at once.

Read-list expansion syntax

A read-list placeholder (reads_1, reads_2, reads_single) can be written plain or with modifiers, to match whatever multi-file syntax your mapper wants:

Whether the first colon-separated part is a prefix or the list name itself is inferred from whether it names a known read list. {cat:name} is a reserved special case in the prefix position -- concatenation only happens once per run (not per iteration) and only for mappers whose template actually uses {cat:...}.

bam_reads: which reads to use when starting from a BAM

Every iteration remaps the same original read pool (extracted once, reused throughout) -- it never shrinks to just whatever mapped last time. When that pool comes from an input BAM rather than FASTQ files, bam_reads controls its scope:

• ref (default, strictest) -- only reads aligned to the chosen reference/contig.
• ref+unal -- that, plus reads that didn't map anywhere (candidates for mapping once the consensus improves).
• all -- every read in the BAM, regardless of what it mapped to.

Your input BAM itself is never rewritten, regardless of any of this. Some of the above needs it indexed, which needs it coordinate-sorted; if it's already sorted, an index is created directly beside it if missing (that's harmless -- purely additive, same as any tool would do), but if it actually needs sorting, that happens to a separate copy under --out-dir, not to your file. (This is specifically about the BAM you pass in via bam =; BAMs iterated-consensus generates itself, like each iteration's mapping output, are sorted in place freely -- those are its own working files.)

Reference resolution

Two [input] keys between them cover every way of specifying a starting reference, for both FASTQ-start and BAM-start:

• reference_id -- just a name, never a file. It can be a record ID within reference_fasta, a contig name within bam, and/or an NCBI accession -- the same string can serve more than one of these roles at once (see the examples below).
• reference_fasta -- a pre-existing local FASTA file. If it has exactly one sequence, that sequence is used automatically; if it has more than one, reference_id must be given to pick which.

Whether you need one, the other, both, or neither depends on the situation:

• Neither: BAM-start, the BAM has only one reference, and its name is itself an NCBI accession (e.g. NC_045512.2) -- it's fetched automatically.
• reference_id only: BAM-start, the BAM has several references, you pick one with reference_id, and that name is itself an accession.
• reference_fasta only: FASTQ-start (or BAM-start), and the file has just one sequence.
• Both: FASTQ-start (or BAM-start) with a multi-sequence reference_fasta -- reference_id picks which record. For BAM-start specifically, this is also how you'd give the actual reference the BAM was aligned against, rather than relying on auto-fetch.
• Neither, and unresolvable: iteration 0 just has no {reference}. That's fine for a BAM-start run unless [consensus] actually uses {reference} -- in which case it's reported as a config error before anything runs, since that combination can never succeed (iteration 0 always runs first, before any consensus this tool computed exists to fall back on). A FASTQ-start run always needs some reference to build iteration 0's mapping index against, so this case is always an error there.

For a BAM-start run, whichever way a reference is obtained, it's validated against the BAM before use: the FASTA record's id must match the resolved contig name exactly, and its sequence length must match the BAM header's length for that contig exactly. A mismatch aborts the run immediately -- proceeding would mean calling a pileup-based consensus against a reference the BAM's own coordinates don't actually match, silently producing garbage.

Convergence: convergence_identity and convergence_streak

Every iteration from 1 onward computes the identity between its new consensus and the previous one (iteration 0 has nothing to compare against, so it's skipped). That identity feeds a streak counter: it increments whenever identity >= convergence_identity, and resets to 0 otherwise. The run stops, reported as converged, once the streak reaches convergence_streak -- i.e. once convergence_streak consecutive iterations have each been at or above convergence_identity. If that never happens, the run stops anyway once max_iterations is reached, just reported as not converged.

With the defaults (convergence_identity = 100.0, convergence_streak = 1), a run stops as soon as one iteration produces a consensus identical to the one before it. Raising convergence_streak (e.g. to 2 or 3) guards against declaring convergence on a fluke -- a sequence could hit exactly 100% once by chance (e.g. a low-coverage region that happens to resolve to the same majority base) and then drift again next iteration; requiring several consecutive matches is a stronger signal that it's genuinely settled. Lowering convergence_identity below 100 is also legitimate -- useful for a sample that may never perfectly stabilize (e.g. a genuinely heterogeneous/mixed population), where "close enough" is a more realistic stopping condition than exact equality.

Final output: [output]

Everything a run produces lives under --out-dir regardless, findable as iter_NNN/consensus.fasta for whichever iteration ran last (see "Output" below) -- [output] is an optional convenience on top of that, for when you want the final result copied somewhere specific rather than having to know which iter_NNN was the last one.

• consensus_fasta -- where to copy it to. Can be anywhere, not necessarily under --out-dir. Its parent directory is created if missing.
• consensus_id -- the FASTA header to give the copy. Optional; if omitted, it keeps whatever id the consensus tool itself assigned. Requires consensus_fasta to also be given -- renaming with nowhere to write doesn't mean anything on its own.

This always uses whichever iteration ran last -- converged or not, since even a run that hit max_iterations without converging usually still has a usable "current best" consensus worth having on hand. It's written (or rewritten) at the end of every run() call, including a resumed run that turns out to already be converged -- so re-running the same command is always safe. Unlike reference_initial.fasta (see "Reference resolution" above), this is always a real copy, never a symlink: it's meant to be a small, standalone, portable deliverable, not a pointer back into --out-dir's own working files.

Output

Iterations are numbered from 0. iter_000 is always the bootstrap step: it produces the first consensus from whatever mapping was already available at the start, rather than one this tool built by iterating. Its exact shape depends on how the run started:

• FASTQ-start: iter_000 builds an index from the reference you gave, maps the reads against it, and calls a consensus -- same shape as every later iteration, just against a reference you supplied rather than one this tool computed.
• BAM-start: iter_000 calls a consensus directly from the input BAM, with no mapping step -- the alignment already exists.

From iter_001 onward, every iteration has the same shape regardless of how the run started: build an index from the previous iteration's consensus.fasta, remap the (always-the-same, extracted-once) reads against it, and call a new consensus. iter_001 is therefore always the first iteration with an identity_to_previous value, since iter_000 has nothing before it to compare against.

results/
  reads/                     extracted/concatenated read files (built once)
  reference_initial.fasta    starting reference (FASTQ-start always; BAM-start
                             if one was resolved) -- see note below
  iter_000/
    <mapper>_index.*         index files (FASTQ-start only -- see above)
    <mapper>.bam              that mapper's mapping output (FASTQ-start only)
    merged.bam                 (only if >1 mapper) merged BAM the consensus step sees
    consensus.fasta             this iteration's consensus
    stats.json                  reads mapped, length, identity to previous, base composition
    logs/                        captured stdout+stderr of every command run
  iter_001/
    <mapper>_index.*         index files, one set per configured mapper
    <mapper>.bam              that mapper's mapping output
    merged.bam                 (only if >1 mapper) merged BAM the consensus step sees
    consensus.fasta             this iteration's consensus
    stats.json                  reads mapped, length, identity to previous, base composition
    logs/                        captured stdout+stderr of every command run
  iter_002/
    ...
  metrics.tsv                 one row per iteration
  summary.json                 iterations run, converged?, total time
  index.html                    human-readable report rendered from summary.json

reference_initial.fasta is a relative symlink straight to your original reference_fasta (or the NCBI-fetched cache file) whenever that's safe -- i.e. it already contains exactly the one sequence needed, nothing else -- rather than a copy, so a large reference genome doesn't get needlessly duplicated. If reference_id had to pick one record out of a multi-sequence reference_fasta, symlinking isn't possible (the file has other sequences in it too), so that case still writes a real single-record copy. Either way, {reference} behaves identically -- every tool that reads it follows the symlink transparently. The BAM itself is never symlinked here: getting from a full input BAM to what iter_000 actually needs (indexed, coordinate-sorted, and -- for bam_reads other than all -- filtered down to reads for the chosen reference) is a real transformation, not a copy, so iteration_0_source.bam under reads/ is always a genuine file.

Resuming a run

If a run stops because it hit max_iterations without converging, raise max_iterations in the config and rerun the exact same command (same --config, same --out-dir): it picks up from the next iteration rather than starting over. This is detected automatically from summary.json in --out-dir -- there's no separate flag. Already-extracted reads and already-completed iterations aren't redone.

If a run already converged, rerunning it against the same --out-dir is a no-op: it reports the existing result without redoing anything.

Resume trusts that --out-dir corresponds to the same logical run -- pointing it at a config with a different mapper, consensus pipeline, or input isn't validated or rejected, it'll just continue on top of whatever's there. Also, resuming is driven entirely by summary.json; a run that crashed before writing it (e.g. killed mid-iteration) can't be resumed and should be started fresh in a new --out-dir.

Known limitations

• --dry-run shows iterations 0 and 1 in full (both always run), but can't show iteration 2 onward -- those commands depend on files that don't exist yet, and whether the run even reaches them depends on convergence.
• A mapper's command template must match the input categories actually present for a given run (e.g. don't reference {reads_single} in a config meant to also run on paired-only data). There's no conditional templating -- write separate configs for meaningfully different input shapes.

Development

uv sync
uv run pytest

Optionally, uv run pre-commit install sets up a pre-commit hook that runs the test suite (and keeps uv.lock in sync with pyproject.toml automatically, regenerating and staging it if a commit -- e.g. a version bump -- leaves it stale) before each commit. This is per-clone setup: it writes into .git/hooks/, which isn't itself tracked by git, so it doesn't happen automatically just from cloning the repo.