Ligase

Type-safe, reproducible computational biology pipelines. A Rust core with a Python surface.

The idea

Biological data types carry dimensions — organism, reference build, coordinate system, sample set — and the engine refuses to connect steps whose dimensions do not match, the way a physics engine refuses to add metres to seconds.

The worst comp bio bugs are silent. hg19 against hg38. 0-based against 1-based. chr1 against 1. Sample IDs that shifted by one somewhere in a join. None of these crash; they produce plausible numbers that are wrong. Workflow managers like Snakemake and Nextflow move files between steps very well, but they do not understand what is in them. Ligase does, and it checks before it computes.

pipe.check()      # type-checks dimensions; no compute
pipe.run()        # parallel, cached, per-sample fan-out
pipe.methods()    # generated methods paragraph + reproducibility bundle

When a pipeline is wrong, you should find out in milliseconds rather than after an overnight run:

DimensionError in step `call`:
  aln.reference = HXB2 (from BAM @SQ header)
  ref.reference = hg38
  These must match. Did you pass the wrong reference?

Status

Pre-alpha, at v0.1's last milestone. M0 (scaffold), M1 (types and dimensions), M2 (the graph), M3 (execution and cache), M4 (semantic hashing and early cutoff), M5 (provenance) and M6 (validation) are complete.

The reference pipeline now interprets as well as computes: its resistance calls reproduce published Stanford HIVdb results exactly at the 15% and 5% thresholds, with every remaining discrepancy explained.

You can write a pipeline of annotated steps, have the whole graph type-checked without running any of it, be told every mistake at once, and then run it — in parallel, cached, one branch per sample:

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)))
pipe.check()        # type-checks every step. Nothing runs.

$ ligase check
checked 8 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 3 other samples: SRR32217715, SRR32217716, SRR32217726.

4 errors in 4 steps.

The step bodies never ran, aln.reference came out of a step that does not exist yet, and a problem that hit the whole cohort is stated once with every affected sample named.

Point it at the right reference and it runs:

$ ligase run
running 12 steps in pipeline.py

ran       `align` (sample SRR8525904)        4.5s
ran       `call` (sample SRR8525904)         513ms
ran       `resistance` (sample SRR8525904)   3ms
...

12 steps: 12 ran.

Run it again and nothing is recomputed, because a step's key is built from its code, its parameters, its inputs' contents and its tools' versions — and none of those moved:

$ ligase run
cached    `align` for SRR8525904, SRR8525908, SRR8525920, SRR8525942
cached    `call` for SRR8525904, SRR8525908, SRR8525920, SRR8525942
cached    `resistance` for SRR8525904, SRR8525908, SRR8525920, SRR8525942

12 steps: 12 from cache.

Change min_freq and only call and resistance come back. Add a sample and only that sample is processed.

Run the aligner again on eight cores instead of one and it reruns — a thread count is a parameter, and nothing can know in advance that it will not change an answer — but the alignments come out the same, so nothing below it is recomputed:

$ ligase run
ran       `align` (sample SRR8525904)   3.3s
ran       `align` (sample SRR8525908)   6.9s
ran       `align` (sample SRR8525920)   3.7s
ran       `align` (sample SRR8525942)   6.0s
skipped   `call` for SRR8525904, SRR8525908, SRR8525920, SRR8525942
skipped   `resistance` for SRR8525904, SRR8525908, SRR8525920, SRR8525942

12 steps: 4 ran, 8 skipped.

"The same" is not byte equality — minimap2 writes -t 8 and the output path into the BAM's @PG line — but what the file means: its records sorted and the fields recording when and how it was made dropped. And when you want the opposite question answered, ligase audit <step> runs a step twice and says whether it repeats.

And when you need to say how the results were produced, ligase methods does — naming the version every tool actually ran at, every parameter as it was hashed, the reference by accession and checksum, and the papers to cite:

$ ligase methods
`align` ran on 4 samples (SRR8525904, SRR8525908, SRR8525920, SRR8525942)
with minimap2 2.31-r1302 and samtools 1.24, as `minimap2 -t {threads} -ax sr
{ref} {reads.r1} {reads.r2} | samtools sort -@ {threads} -o {out}`, with
threads = 4, against the HXB2 reference (K03455.1, 9719 bp, MD5
c66838498aab02f786a6ec2f2d62cdc1).

It can say that about a run where nothing was recomputed, which is the point: the command is recorded when a step runs, because since M4 it is not in the output's hash and a cache hit can no longer prove it. --bundle packages all of it so someone else can check the results, and --verify is that check.

What remains is M6 (validation): whether the resistance calls agree with the published ones. See milestones for the order and the acceptance test for each step.

Where to go

If you want to	Read
Install it and run the checks	Installation
Understand how the project is put together	Project layout
Know why the workspace is shaped this way	Workspace architecture
See the target API	Golden example
Understand why data carries dimensions	Types and dimensions
Declare dimensions and read the errors	Declaring dimensions
Write a pipeline and check it	Building a pipeline
Run one, and understand what was cached	Running a pipeline
Understand why steps are never run to build a graph	Pipelines and the graph
Know how the checker works inside	Type system
Know how the graph works inside	Graph
Understand why a step did or did not rerun	Running, and not running
Know why two different files can be the same result	What a file means
Check that a step repeats itself	Semantic hashing and ligase audit
Know how a run is scheduled	Execution
Know how the cache and its keys work	Cache
Know how a canonical form is built	Canonicalisation
Say how a result was produced, for a paper	How was this produced?
Understand a noticed line in a run report	When a tool names a sample after a file
Fix an output named after a file rather than a sample	Reading what the engine noticed
Write a methods paragraph and a bundle	Writing up a pipeline
Know how provenance is stored and checked	Provenance
Turn variant calls into drug-resistance calls	Resistance interpretation
Write a resistance step, and pick a threshold	Interpreting resistance
Know how the HIVdb rules are read and scored	Resistance interpretation
See what the calls were checked against	Validation
Know which dimensions exist and how they are inferred	Dimensions
Follow the build order	Milestones
Understand a design choice	Decisions

Names

Thing	Name
Rust core crate	ligase
PyPI distribution	ligase-bio
Python import	ligase

Licence

Dual-licensed under either MIT or Apache-2.0, at your option.