Steps you do not have to write

ligase.steps is curated rather than convenient: a step is in it only once it declares its tool’s version capture, its citation, the dimensions of its output and how that output canonicalises (D42). Four of them, and the clause has already kept two of the obvious candidates out — see what a curated step is.

steps.minimap2.align

align(
    reads: 'Reads',
    ref: 'Reference',
    sample_id: 'str',
    threads: 'int' = 4,
) -> 'Alignment'

Align short reads to a reference, coordinate-sorted and read-grouped.

This is minimap2 -ax sr followed by the two samtools invocations that make its output usable, and the reason it is one step rather than three is that the three are not independently useful: an alignment without a read group has no sample, and an unsorted one cannot be piled up or indexed. Both are declared in tool=, so both versions join the cache key.

Why the read group is written here

minimap2 writes no @RG, and a BAM with no read group has no sample name. bcftools then falls back to naming the VCF’s sample column after the input file’s path — which in a Ligase run is a file inside the content-addressed store, so the calls would carry a cache path where a sample identifier belongs. That is wrong three times over: wrong in a methods paragraph, wrong for the sample-set dimension, and wrong for reproducibility, because the canonical form of a VCF includes its sample names and a cache path is not the same in two places.

So the sample is not optional and is not inferred. It is the one thing this step asks of you that the files cannot answer.

Paired-end reads

This step aligns a pair. Handed single-end reads it fails before any compute, with the engine saying so: {reads.r2}: reads is a single-end sample, so it has no second mate. That the requirement is not in the signature is a gap rather than a choice — Reads is one kind covering both shapes, and nothing lets a step require a value of the pairing dimension the way coordinates= lets it declare one. The architecture page  has the rest of that story; it is the same missing mechanism that keeps a paired-end trimmer out of this library.

Why threads is a parameter like any other

Changing it reruns this step, because the engine cannot know in advance that a thread count will not change the answer. What it can know, once the step has rerun, is whether the result changed: an alignment is hashed by what it means rather than by its bytes, and the @PG line where minimap2 records -t 8 is not part of what it means. So the alignment comes out with the same canonical hash and everything below this step is skipped.

Parameters

Returns

A recipe, typed as the coordinate-sorted alignment the step produces.

steps.bwa.mem

mem(
    reads: 'Reads',
    ref: 'Reference',
    sample_id: 'str',
    threads: 'int' = 4,
) -> 'Alignment'

Align short reads with BWA-MEM, coordinate-sorted and read-grouped.

The same output as ligase.steps.minimap2.align and deliberately not the same step. Two aligners are not interchangeable in provenance: they are different programs with different versions and different papers, so swapping one for the other changes the cache key and changes the methods paragraph. Neither of those is arranged here — both fall out of tool= naming bwa where the other names minimap2.

The read group is written by samtools addreplacerg rather than by bwa’s own -R, for no deeper reason than that it is then written the same way in both steps, and one mechanism has one place to be wrong. Why it has to be written at all is in ligase.steps.minimap2.align’s docstring.

The index, and what it costs

bwa mem cannot read a bare FASTA: it needs the five files bwa index writes. In a Ligase run the reference is a file in the content-addressed store, where a step has no business writing anything — an object there is named after its own bytes, and a step that left five siblings beside one would be putting unaddressed files into an addressed store.

So this step builds the index into a temporary directory of its own and throws it away afterwards. That is correct and it is slow: the index is rebuilt on every run of the step, which for a 2 kb reference costs nothing and for a human genome costs an hour. minimap2 has no such cost because it indexes from the FASTA in memory.

The honest fix is not a flag on this step. It is the same engine feature two other things want — a step that can produce more than one file — after which an index is something a step produces once and others consume. Until then: this step is the right one for a small reference, and for a large one, align with minimap2 or write a step of your own that points at an index you keep.

Paired-end reads

As with the other aligner, this one aligns a pair; see its docstring for why the requirement is not in the signature.

Parameters

Returns

A recipe, typed as the coordinate-sorted alignment the step produces.

steps.samtools.sort

sort(aln: 'Alignment', threads: 'int' = 4) -> 'Alignment'

Put an alignment in coordinate order.

Both aligners in this library sort what they produce, so this step is for an alignment that came from somewhere else: a BAM a collaborator sent, one an older pipeline wrote, one a step of your own produced. Nearly everything downstream wants coordinate order — bcftools mpileup does, and so does any form of indexing.

Sorting changes the bytes and not the meaning

Which is a stronger statement than it sounds, and it is the engine’s rather than this step’s. An alignment’s canonical form sorts its records before hashing (D4), because the order records sit in a file is not what the file is about. So a sorted alignment and its unsorted original have the same canonical hash — and inserting this step into a pipeline therefore reruns the sort and skips everything below it, by early cutoff, with no special case anywhere.

The practical reading: a sort costs you one step’s work once, and never costs you the pipeline below it.

Why samtools index is not in this library

Because it cannot say what it produces, and D42 admits nothing that cannot. samtools index writes a .bai beside its input. Declared Alignment → Alignment it would be false twice: the output file is an index rather than an alignment, and the engine would not catch it — a file that does not parse as its kind falls back to having its bytes hashed, so an index’s bytes would stand in for an alignment’s meaning and every step below would accept it as an alignment.

Nor is the answer a new kind. An index brings no dimension, and under the scope rule a kind that brings no dimension does not earn a place here (D35). What an index actually wants is to travel beside the alignment it indexes, and a step’s output is one file, so there is nowhere for it to travel. That is an engine gap, recorded as one on the architecture page , rather than a step shipped with a false type on it.

Parameters

Returns

A recipe, typed as the coordinate-sorted alignment the step produces.

steps.bcftools.mpileup_call

mpileup_call(
    aln: 'Alignment',
    ref: 'Reference',
    min_freq: 'float' = 0.01,
) -> 'Variants'

Report every allele seen in at least min_freq of the reads.

There is no bcftools call in this step, and the name says mpileup_call because that is what people call this operation rather than because a caller runs. The pileup’s own allele depths are the evidence, and min_freq is applied straight to them.

Why the caller is gone

bcftools call -mv is a consensus caller: it asks what the sample’s genotype is, and a threshold below one half rarely removes anything, because an allele the majority does not carry is not a call it will make. For anything where the minority population is the question — drug resistance, subclonal variants, a mixed infection — that is the wrong question, and the reference demo is where the cost was measured.

On one real sample the reads at a resistance position were 238 reference and 192 alternate: a 45% mutation, clinically the whole answer, and bcftools call -mv emitted no record for it at all, because the majority allele agreed with the reference. The pipeline reported the sample susceptible to every drug. It was right given that input, and it matched the published majority-consensus report exactly — and the published 5% report for the same sample read high-level resistance to four drugs. A pipeline whose answer flips between “susceptible” and “resistant to most of a class” on a threshold it never exposed is not doing what it says it does.

So the threshold is exposed and means what its name says. It is a typed, hashed parameter, so changing it reruns exactly the steps it could have changed.

The index it does not leave behind

bcftools mpileup -f needs a .fai for the reference and writes one beside the path it was given if there is none. In a Ligase run that path is an object in the content-addressed store, where every name is the hash of its own bytes — so writing an index there would be putting an unaddressed file into an addressed store, and on a store that is read-only or remote it would not be writing anything at all, it would be failing. So the step indexes a link in a directory of its own and throws both away.

What it does not decide

min_freq is about what the reads show. What counts as present for the purpose of acting on it is a separate question with a separate parameter, wherever the interpretation happens — ligase.hivdb.interpret’s mixture, for the demo. Conflating the two is how a threshold ends up unexposed.

And a position this step reports nothing about is ambiguous between agreeing with the reference and nothing having aligned there. The two have opposite meanings and a VCF cannot tell them apart, so whatever reads these calls needs the alignment too; ligase.hivdb.interpret takes it for exactly that reason.

Parameters

Returns

A recipe, typed as the variant calls the step produces.

steps.LIBRARY

steps.LIBRARY: tuple = (LibraryStep(step=<Step align -> Alignment>, module='ligase.steps.minimap2', canonical='canonical', page='/guide/steps/minimap2/'), LibraryStep(step=<Step mem -> Alignment>, module='ligase.steps.bwa', canonical='canonical', page='/guide/steps/bwa/'), LibraryStep(step=<Step sort -> Alignment>, module='ligase.steps.samtools', canonical='canonical', page='/guide/steps/samtools/'), LibraryStep(step=<Step mpileup_call -> Variants>, module='ligase.steps.bcftools', canonical='canonical', page='/guide/steps/bcftools/'))

steps.LibraryStep

class LibraryStep(
    step: 'Step[Any,
    Any]',
    module: 'str',
    canonical: 'str',
    page: 'str',
) -> None

One entry in the library, and what it claims about itself.

The claims are kept as data rather than as prose so that one test can iterate them. An admission criterion nothing checks is a comment, and the four in this module’s docstring are meant to hold a contributor who has not read it.

Attributes

steps.LibraryStep.name

name: str  # property

How the step is written, module and all: minimap2.align.