Quickstart

pip install ligase-bio
mkdir ligase-quickstart && cd ligase-quickstart

1. Two files that do not agree

A tiny reference, and an alignment against it:

mkdir refs aln
cat > refs/genome.fasta <<'EOF'
>1 contig one
ACGTACGTACGTACGTACGT
>2 contig two
TTTTCCCCGGGGAAAATTTT
EOF

printf '@HD\tVN:1.6\tSO:coordinate
@SQ\tSN:chr1\tLN:20\tM5:a965a71aa3690f605935c54d320905ab
@SQ\tSN:chr2\tLN:20\tM5:89c2d6e069f386bafb1a0c811ba821f4
@RG\tID:rg0\tSM:sample01\tPL:ILLUMINA
read1\t0\tchr1\t1\t60\t10M\t*\t0\t0\tACGTACGTAC\tIIIIIIIIII\tRG:Z:rg0
' > aln/sample01.sam

Look carefully and you can see the problem: the FASTA calls its contigs 1 and 2, the SAM calls them chr1 and chr2. Same sequence, same lengths, different names. In a real project those two files come from different places months apart, and nobody looks carefully.

2. Inspect them first

Before writing anything at all, ask what Ligase makes of the two files:

ligase inspect refs/genome.fasta aln/sample01.sam

That is the whole of Ligase in one command, and it needed no pipeline. Read it in three parts.

Each file gets a block headed by what it was read as. Inside it, the dimensions Ligase worked out — 2 contigs, 40 bp, chr-prefixed (UCSC style), [sample01] — each with where it read that from, which is the part that lets you go and look rather than guess.

Then the dimensions that have no value, with the reason. Those are an invitation: organism is unknown (a FASTA does not record its organism; declare it) is how you find out what declaring is for. And the last line of each block names the dimensions that do not apply at all — a FASTA is not missing a quality encoding, it has no such property to have.

Finally the comparison: the same ContigNamingError a pipeline would have given you, and a line saying which dimensions the verdict covers. It exits non-zero, so ligase inspect refs/*.fasta aln/*.bam && snakemake refuses to start a run over files that disagree.

3. A pipeline

"""A two-step pipeline, with a mistake in it."""
 
from pathlib import Path
 
import ligase as lg
 
 
@lg.step
def call(aln: lg.Alignment, ref: lg.Reference) -> lg.Variants:
    """Pretend to call variants. Real pipelines shell out; this one does not."""
 
    def write(out: Path) -> None:
        out.write_text(
            "##fileformat=VCFv4.3\n"
            "##contig=<ID=1,length=20>\n"
            "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tsample01\n"
            "1\t5\t.\tA\tG\t60\tPASS\t.\tGT\t1\n"
        )
 
    return lg.write(write)
 
 
ref = lg.reference("refs/genome.fasta", organism="test")
aln = lg.alignment("aln/sample01.sam", organism="test")
 
pipeline = lg.Pipeline(call(aln, ref))

Three things are doing work here.

The type annotations are not decoration. aln: lg.Alignment and ref: lg.Reference are what let the engine know what each argument is and which dimensions apply to it, and -> lg.Variants is how it knows what comes out. The graph is built from these signatures.

lg.write hands the engine a function to call with the path to write to, so this step needs no external tool. A real step would use lg.sh and name the tool it runs.

lg.Pipeline(...) builds the graph. It does not run anything — calling call(aln, ref) records a call rather than performing one.

4. Check it

ligase check

The same mismatch as in step 2, and that is the point: inspect compared two files you named, while check compared the inputs the pipeline itself wires together — every connection, every step, every sample, under the argument names your step declares. On a pipeline over four hundred samples that is a different question with the same answer, and the one run refuses to proceed past.

It took milliseconds, and the step body never ran. Without the check, this pipeline would have succeeded: tools match contigs by name, so the variant caller would have found no overlap between the alignment and the reference, written a valid empty VCF, and exited zero.

5. Fix it

Rename the contigs in the SAM to match the reference:

sed -i.bak 's/SN:chr/SN:/; s/\tchr1\t/\t1\t/' aln/sample01.sam
ligase check

6. Run it

ligase run

Now run it again, without changing anything:

Nothing recomputed. A step’s key is built from its code, its parameters, its inputs’ contents and its tools’ versions — and none of those moved. Edit the step body, or a byte of the input, and it runs again.

7. Say how it was produced

ligase methods

On a real pipeline that paragraph names every tool at the version it actually ran at, every parameter as it was hashed, and every reference by accession and checksum. --bundle packages it so somebody else can check the results.

What to read next