Building a pipeline¶

Everything here runs today. Nothing in it executes a step: the whole point is to be told your pipeline is wrong before you wait for it.

Write a step¶

A step is an annotated function with @lg.step on it:

import ligase as lg
from ligase.types import Alignment, Reads, Reference, Variants


@lg.step(tool="minimap2")
def align(reads: Reads, ref: Reference) -> Alignment:
    return lg.sh("minimap2 -ax sr {ref} {reads.r1} {reads.r2} | samtools sort -o {out}")


@lg.step(tool="bcftools")
def call(aln: Alignment, ref: Reference, min_freq: float = 0.01) -> Variants:
    return lg.sh("bcftools mpileup -f {ref} {aln} | bcftools call -mv -Oz -o {out}")

The annotations are not documentation. They are the only description of the step that Ligase ever reads, because it never runs the body while building a graph:

a parameter annotated with a Ligase data type becomes a connection in the graph
anything else (min_freq: float = 0.01) is a plain parameter, which joins the cache key, so changing it reruns exactly the steps it could have changed
the return annotation is what says the step produces an Alignment

tool= names the external program the step shells out to — or several of them, tool=("minimap2", "samtools"), for a body that pipes one into the other. Each one's version belongs in the cache key, so upgrading minimap2 reruns the alignments it produced. See running a pipeline.

Leave an annotation off and you are told immediately, as the module is imported:

>>> @lg.step
... def align(reads: Reads, ref) -> Alignment: ...
StepDefinitionError: step `align` does not annotate `ref`. An unannotated
parameter might be data, and data the graph does not know about is an
unchecked connection

Connect them¶

Calling a step records the call and returns a handle. Nothing runs:

hxb2 = lg.reference("data/refs/HXB2.fasta", organism="HIV-1", build="HXB2")
reads = lg.reads("data/runs/SRR32217714_1.fastq.gz", "data/runs/SRR32217714_2.fastq.gz")

pipe = lg.Pipeline(call(align(reads, hxb2), hxb2))
pipe.check()        # raises on the first problem; nothing is computed

lg.Pipeline walks back from the output through the handles, so the graph is built from the outputs you ask for. A file passed to two steps is one node, because it is one file:

>>> pipe = lg.Pipeline(call(align(reads, hxb2), hxb2))
>>> pipe.steps
('align', 'call')
>>> len(pipe)       # two inputs, two steps -- not three inputs
4

Check it¶

check() raises the first problem and attaches the rest. errors() returns them without raising, which is usually what you want in a notebook:

>>> for error in pipe.errors():
...     print(error)
...     print()
DimensionError in step `call`:
  aln.reference = HXB2 (from step `align`)
  ref.reference = hg38 (from FASTA contigs: data/refs/hg38.fasta)
  Contig `1`: same name and length but different sequence (MD5 a965a71a against 44e7acaf).
  These must match. Did you pass the wrong reference?

One call step, so the heading is just the step. Once you fan out, one written step is many steps and the heading says which — see many samples below.

aln.reference came out of a step that has not run. The dimensions of a step's output are worked out from its inputs, so a mistake made at the top of a pipeline is reported wherever it first matters, with the step it travelled through named.

Every error carries the whole list, so catching one shows you all of them:

>>> try:
...     pipe.check()
... except lg.CheckError as caught:
...     print(len(caught.errors), "problems")
2 problems

Catch	To get
`lg.CheckError`	everything
`lg.DimensionError`	data that disagrees
`lg.ContigNamingError`	`chr1` against `1` — add a renaming step
`lg.ContigOrderError`	contigs in a different dictionary order
`lg.UndeclaredDimensionError`	something to declare where you load data
`lg.StepTypeError`	a step wired to the wrong kind of data

Many samples¶

lg.samples() globs, groups mates into samples and loads each one. .map() applies a chain to every sample:

samples = lg.samples("data/runs/*.fastq.gz", paired=True)
pipe = lg.Pipeline(samples.map(lambda s: call(align(s, hxb2), hxb2)))

>>> samples.ids
('SRR32217714', 'SRR32217715', 'SRR32217716')
>>> len(pipe)       # 3 samples + 1 reference + 3 aligns + 3 calls
10

When every sample fails the same way, the report says so once and names them all:

$ ligase check
checked 6 steps in pipeline.py

DimensionError in step `call` (sample SRR32217714):
  aln.reference = HXB2 (from step `align`)
  ref.reference = hg38 (from FASTA contigs: data/refs/hg38.fasta)
  Contig `K03455.1` is not in hg38.
  These must match. Did you pass the wrong reference?
  The same problem affects 2 other samples: SRR32217715, SRR32217716.

3 errors in 3 steps.

pipe.errors() is not collapsed — one error per failing step, each with its own sample in error.scope, which is the heading in a form you can compute with:

>>> error.scope
('sample', 'SRR32217714')
>>> sorted({e.scope[1] for e in pipe.errors() if e.scope})
['SRR32217714', 'SRR32217715', 'SRR32217716']

pipe.report().text is the collapsed form the command line prints, if you want to show it somewhere yourself.

scope is None whenever the step's own name already identifies it, which is every step in a pipeline without fan-out. Where a name is shared, the discriminator is whatever actually separates those steps — ('sample', ...) under fan-out, ('file', ...) for same-named steps over different data from one sample — and None if nothing does, because a suffix that repeats would look like an identifier without being one.

Mates are paired by filename — <sample>_1/<sample>_2 or <sample>_R1/<sample>_R2, with anything after the marker dropped, so Illumina's SAMPLE_R1_001.fastq.gz works. A file whose mate cannot be identified is an error rather than a guess:

>>> lg.samples("data/runs/*.fastq.gz", paired=True)
LoadError: cannot group reads into samples: sample `SRR32217714` has only one
of its two mates (data/runs/SRR32217714_1.fastq.gz); pass paired=False for
single-end reads

That is deliberately unhelpful about carrying on. Pairing two samples' files by accident produces results for a mixture of two people and nothing crashes.

Intervals¶

A BED or GFF loads like anything else, and brings the coordinates dimension with it:

>>> regions = lg.intervals("data/regions.bed")
>>> regions.dimension("coordinates").value
'0-based half-open'
>>> [c.name for c in regions.contigs]
['K03455.1']

The contigs come from the records, so they have no length and no checksum — a BED records neither. They are still compared against a reference by name, which catches a chr-prefixed BED against a bare reference, and an interval on a contig the reference does not have.

Feed those intervals to a step that also takes 1-based data and you get the error the dimension exists for:

DimensionError in step `restrict`:
  v.coordinates = 1-based closed (from step `call`)
  regions.coordinates = 0-based half-open (from BED coordinate convention: data/regions.bed)
  These must match. Did you mean to convert between coordinate systems?

If a step of yours writes intervals, say which convention it writes. Alignments and variant calls are 1-based closed by specification, but a tool that writes BED from a VCF changes the convention, and inheriting the input's would be a silent lie:

@lg.step(tool="bedtools", coordinates="0-based half-open")
def windows(v: Variants) -> Intervals:
    return lg.sh("bedtools makewindows ...")

From the command line¶

Put the pipeline in a file and let the CLI check it:

# pipeline.py
import ligase as lg
from ligase.types import Alignment, Reads, Reference, Variants

@lg.step(tool="minimap2")
def align(reads: Reads, ref: Reference) -> Alignment:
    return lg.sh("minimap2 -ax sr {ref} {reads.r1} {reads.r2} | samtools sort -o {out}")

@lg.step(tool="bcftools")
def call(aln: Alignment, ref: Reference, min_freq: float = 0.01) -> Variants:
    return lg.sh("bcftools mpileup -f {ref} {aln} | bcftools call -mv -Oz -o {out}")

hxb2 = lg.reference("data/refs/HXB2.fasta", organism="HIV-1", build="HXB2")
samples = lg.samples("data/runs/*.fastq.gz", paired=True)

pipe = lg.Pipeline(samples.map(lambda s: call(align(s, hxb2), hxb2)))

ligase check

$ ligase check
checked 2 steps in pipeline.py

no problems found.

With no argument it reads pipeline.py in the working directory; pass a path for anything else. The file is imported — that is how a Python pipeline declares itself — so the lg.Pipeline has to be assigned at the top level where the CLI can find it. The step bodies are still not run.

Exit codes separate the two failures worth separating:

Code	Meaning
0	every pipeline in the file checks clean
1	a pipeline has errors
2	the file could not be checked at all

That makes it a CI gate:

- run: ligase check

Next¶

Pipelines and the graph — why it works this way
Declaring dimensions — filling the gaps a file cannot state
Graph architecture — propagation, requirements and attribution