Skip to content

Types and dimensions

What this is

A Ligase dataset is not just a file path. It is a kind of biological data — reads, a reference, an alignment, variant calls — carrying the dimensions of what is inside it: which organism, which reference build, which coordinate convention, which samples in which order.

import ligase as lg

hxb2 = lg.reference("data/refs/HXB2.fasta", organism="HIV-1", build="HXB2")
aln = lg.alignment("data/alignments/SRR32217714.bam")

lg.check_step("call", {"aln": aln, "ref": hxb2})  # nothing computed

check_step compares the dimensions and raises if they disagree. No compute runs, no files are written, and the whole thing takes milliseconds.

Why it exists

The worst computational biology bugs are silent. Aligning to hg19 and annotating against hg38. Feeding 0-based half-open intervals to a tool that expects 1-based closed. chr1 against 1. A sample column that shifted by one somewhere in a join. None of these crash. All of them produce plausible numbers that are wrong, and the wrongness can survive peer review.

Workflow managers move files between steps extremely well, but a file path is an opaque token to them: they cannot tell you that the BAM in slot one and the FASTA in slot two describe different genomes. Ligase reads the headers, learns what the dimensions are, and refuses the connection.

The physics analogy is exact. You would not add metres to seconds. A reference build is a unit, and so is a coordinate convention.

The dimensions

The full catalogue, including which formats can state which dimension, is in Dimensions. In brief:

Dimension What it captures
organism HIV-1, SARS-CoV-2, human
reference The build, by contig identity — not by name
contig_naming chr1 against 1
coordinates 0-based half-open against 1-based closed
sample_set The sample identifiers, in order
pairing Single-end against paired-end
quality_encoding Phred+33 against Phred+64

Each dimension is inferred from the file where the format records it, and declared by the user once where it does not. A FASTA cannot tell you its organism; a BAM @SQ header can tell you its reference exactly.

Three ideas that make this work

Reference identity is contigs, never filenames

A reference's identity is its contig names, their lengths, and their MD5 checksums where available. Not the filename, and not a label you gave it.

This matters because filenames lie constantly. genome.fa on two machines is two different genomes. A directory called hg38 can hold hg19. Worse, a label can be right and still useless: hg19 and hg38 both have a chr1, and distinguishing them needs the sequence itself.

So Ligase computes each contig's MD5 from the FASTA, reads the @SQ M5 an aligner wrote into a BAM, and compares them. If an aligner claims your BAM is on a reference whose chr1 has a different checksum from the chr1 in the FASTA you just passed, that is caught before anything runs.

A label such as HXB2 is kept, but only to make error messages readable. When two references describe themselves identically — two files both labelled hg38 — the message falls back to showing their fingerprints, because a message that prints the same value twice cannot be acted on.

Provenance travels inside the value

Every dimension value carries where it came from. Not reconstructed when an error is built; stored alongside the value from the moment it is read.

That is the difference between an error you can act on and one that sends you hunting:

DimensionError in step `call`:
  aln.reference = HXB2 (from BAM @SQ header: data/sample.bam)
  ref.reference = hg38 (from FASTA contigs: refs/hg38.fasta)
  These must match. Did you pass the wrong reference?

You know which two files to look at and which part of each one said what.

Unknown is a real answer

Where a format cannot state a dimension, Ligase says so and says why, rather than guessing:

>>> ref = lg.reference("data/refs/HXB2.fasta")
>>> ref.dimension("organism") is None
True
>>> ref.unknown_reason("organism")
'a FASTA does not record its organism; declare it'

A guess that happens to be wrong passes check() silently, which is the exact failure mode this project exists to remove. So a FASTQ whose quality scores fall in the range Phred+33 and Phred+64 share reports its encoding as unknown; guessing would shift every quality value by 31 and nothing downstream would notice.

"Unknown" is distinct from "does not apply". A reference has no quality encoding at all — that is a category error, not a gap — and check_step never asks you to fill it in.

When a dimension must be known

Ligase draws a line between two things:

  • Checked — whenever two inputs both state a value, those values must agree. Always, for every dimension. This costs nothing and demands nothing.
  • Required — additionally, the value must not be unknown. Only where a disagreement would silently corrupt results.

Checking an alignment against a reference requires the reference dimension: an alignment that cannot say what it was aligned to is not something to shrug at. But aligning reads to a reference requires nothing, because introducing the reference is what alignment does.

This line is deliberate, and it errs toward demanding less. A checker that cries wolf gets bypassed, and a bypassed checker catches nothing. Correct input that gets rejected is not a safe failure — it is how users learn to stop trusting the tool.

Coordinates, and the type that needed them

coordinates sat in the catalogue doing almost nothing for a milestone, because SAM, BAM and VCF are all 1-based closed and a lossless re-encoding between containers changes nothing (D9).

Intervals is what it was waiting for. BED is 0-based half-open and GFF is 1-based closed, and feeding one to a tool that expected the other produces results that are off by exactly one base and look entirely plausible:

>>> regions = lg.intervals("data/regions.bed")
>>> regions.dimension("coordinates").value
'0-based half-open'

An interval file's contigs come from its records, so they carry no length and no checksum — a BED records neither, and inventing them would be the lie this whole system guards against. They are still compared against a reference by name, which is enough to catch a chr-prefixed BED against a bare reference, or an interval on a contig the reference does not have.

Their order, on the other hand, is not compared at all. The contigs of a sequence dictionary have a meaningful order — GATK requires two dictionaries to agree on it — but the contigs of a BED's records appear in whatever order the intervals were sorted in. Reporting that as a contig order mismatch would be a false alarm, and false alarms are how a correctness tool gets switched off.

What is not here yet

Dimensions are checked where steps meet. They are not converted: a typed liftover that the engine could search for to fix a mismatch automatically is v2, and the design note is D5.

Next