Applied - Reproducible Environment

Core Principle: Computational reproducibility in bioinformatics requires strict version control across Python dependencies, containerized runtime environments, and external reference genome database releases.

Uses: Versions and Compatibility · UTA and Data Providers · SeqRepo · Biopython Deprecations · NumPy Random Generator

The three layers

1. CODE        your scripts: git
2. SOFTWARE    Python + packages: a lockfile
3. DATA        genome build, annotation release, UTA/SeqRepo snapshot,
               gnomAD version, MANE version   ← the layer people forget

Layer 3 is where bioinformatics differs from general software. An analysis that pins hgvs==1.5.7 but not the UTA snapshot is not reproducible, because normalization output depends on the sequence data. Same for GENCODE releases, gnomAD versions, and genome patches.

Layer 2: pinning software

# pyproject.toml
[project]
requires-python = ">=3.12,<3.14"
dependencies = [
    "numpy>=2.4,<3",
    "pandas>=3.0,<4",
    "pyarrow>=15",
    "matplotlib>=3.10,<4",
    "seaborn>=0.13,<0.14",
    "biopython>=1.86,<2",
    "hgvs>=1.5,<1.6",
    "scipy>=1.14",
    "statsmodels>=0.14",
    "pysam>=0.22",
]

Ranges for development, an exact lockfile for the actual run:

uv lock                       # uv.lock: exact versions and hashes
uv sync --frozen
uv pip freeze > requirements.lock.txt

Or with conda/mamba, which is usually necessary in bioinformatics because you also need non-Python tools:

# environment.yml
name: bioenv
channels: [conda-forge, bioconda]
dependencies:
  - python=3.12
  - numpy=2.5
  - pandas=3.0
  - biopython=1.87
  - hgvs=1.5.7
  - pysam=0.22
  - samtools=1.21
  - bcftools=1.21
  - mafft=7.525
  - pip

mamba env create -f environment.yml
mamba env export --no-builds > environment.lock.yml

Bioconda is where the command-line tools live, and mixing pip and conda in one environment reliably breaks. Pick one as primary.

Why the pins matter here specifically: see Versions and Compatibility:

• pandas 3.0 made chained assignment a silent no-op
• Biopython 1.86 removed Bio.Application
• NumPy 2 changed scalar promotion (NEP 50), altering float32 results
• hgvs asks you to pin to a minor range because normalization can change

Layer 3: pinning data

This is the part that distinguishes a reproducible bioinformatics analysis.

REFERENCE_DATA = {
    "genome_build": "GRCh38.p14",
    "genome_fasta": "GCA_000001405.29_GRCh38.p14_genomic.fna.gz",
    "genome_md5": "a08035b6a6...",
    "annotation": "GENCODE v45",
    "annotation_url": "https://ftp.ebi.ac.uk/pub/databases/gencode/...v45.annotation.gtf.gz",
    "mane": "MANE v1.4",
    "uta_schema": "uta_20210129",
    "seqrepo_snapshot": "2024-02-20",
    "gnomad": "v4.1",
    "clinvar": "2026-08-01",
}

"We used hg38" is not sufficient. Patch releases add and remove alt contigs, and a variant on a contig present in p13 but not p12 simply vanishes.

Capture it at runtime:

import os, json, subprocess, hashlib, sys
from datetime import datetime, timezone
from pathlib import Path
import numpy, pandas, matplotlib, seaborn, Bio, hgvs

def provenance(hdp=None, extra=None):
    p = {
        "run_at": datetime.now(timezone.utc).isoformat(),
        "python": sys.version.split()[0],
        "platform": sys.platform,
        "packages": {
            "numpy": numpy.__version__, "pandas": pandas.__version__,
            "matplotlib": matplotlib.__version__, "seaborn": seaborn.__version__,
            "biopython": Bio.__version__, "hgvs": hgvs.__version__,
        },
        "git_commit": subprocess.run(
            ["git", "rev-parse", "HEAD"], capture_output=True, text=True
        ).stdout.strip() or "not-a-repo",
        "git_dirty": bool(subprocess.run(
            ["git", "status", "--porcelain"], capture_output=True, text=True
        ).stdout.strip()),
        "reference_data": REFERENCE_DATA,
        "seqrepo_dir": os.environ.get("HGVS_SEQREPO_DIR"),
        "uta_db_url": os.environ.get("UTA_DB_URL", "remote-default"),
        "random_seed": SEED,
    }
    if hdp is not None:
        p["uta_data_version"] = hdp.data_version()
    if extra:
        p.update(extra)
    return p

Path("provenance.json").write_text(json.dumps(provenance(hdp), indent=2))

git_dirty is the underrated field. A commit hash from a dirty working tree is a lie, and this is how you find out.

Checksums on inputs

def sha256(path, chunk=1 << 20):
    h = hashlib.sha256()
    with open(path, "rb") as fh:
        while (b := fh.read(chunk)):
            h.update(b)
    return h.hexdigest()

INPUT_CHECKSUMS = {p.name: sha256(p) for p in Path("data/raw").glob("*")}

def verify(path, expected):
    actual = sha256(path)
    if actual != expected:
        raise ValueError(f"{path} changed: {actual} != {expected}")

Reference files get silently updated, re-downloaded, or symlinked to a different version. A checksum turns that into an error at load time rather than an unexplained difference in your results six months later.

Random seeds

Use np.random.default_rng, seed once at the top, pass the generator down, and record the seed in the provenance. For parallel workers use SeedSequence.spawn, not seed + i. See NumPy Random Generator.

Also seed anything else that samples: random.seed, sklearn random_state, PYTHONHASHSEED.

Containers

FROM mambaorg/micromamba:1.5

COPY environment.lock.yml /tmp/env.yml
RUN micromamba install -y -n base -f /tmp/env.yml && micromamba clean --all --yes

ENV HGVS_SEQREPO_DIR=/data/seqrepo/2024-02-20
ENV UTA_DB_URL=postgresql://anonymous@uta:5432/uta/uta_20210129

# docker-compose.yml
services:
  uta:
    image: biocommons/uta:uta_20210129        # ← the tag IS the data version
    ports: ["5432:5432"]
  analysis:
    build: .
    volumes: ["./data:/data", "./results:/results"]
    depends_on: [uta]

Containerising UTA is the neat part: the image tag pins the data version as firmly as a package version.

Reference data is too large for an image; mount it and record its version in the environment. Pin base images by digest, not by tag: tags move.

Project layout

project/
├── environment.lock.yml
├── pyproject.toml / uv.lock
├── Snakefile              # or nextflow.nf
├── config.yaml            # ← all paths and versions, no hardcoding
├── src/
│   ├── plotstyle.py       # matplotlib/seaborn defaults: see [[rcParams and Style Sheets]]
│   └── pipeline.py
├── data/
│   ├── raw/               # read-only, checksummed, never edited
│   ├── interim/           # .parquet caches: see [[Reading and Writing Data]]
│   └── external/          # reference files, version in the path
├── results/
│   ├── figures/
│   ├── tables/
│   └── provenance.json
└── tests/

data/external/gencode_v45/...: put the version in the path. A directory called annotation/ tells you nothing in a year.

data/raw/ is read-only. Every transformation writes somewhere else. This is the single most valuable convention in the list.

Workflow managers

# Snakefile
rule annotate:
    input:  vcf="data/raw/{sample}.vcf.gz",
            gnomad="data/external/gnomad_v4.1/af.parquet"
    output: "results/tables/{sample}_annotated.tsv"
    conda:  "envs/hgvs.yml"
    log:    "logs/{sample}_annotate.log"
    shell:  "python src/annotate.py {input.vcf} {input.gnomad} {output} 2> {log}"

Snakemake and Nextflow give you dependency tracking, resumability, per-rule environments, and provenance reports. Worth the learning cost past about five steps.

Caching intermediates

def cached(path, compute):
    p = Path(path)
    if p.exists():
        return pd.read_parquet(p)
    df = compute()
    p.parent.mkdir(parents=True, exist_ok=True)
    df.to_parquet(p)
    return df

variants = cached("data/interim/variants.parquet",
                  lambda: read_vcf("data/raw/cohort.vcf.gz"))

Parquet preserves dtypes; CSV re-infers them on every read, which means a column can change type between runs. See Reading and Writing Data.

Invalidate the cache when the input's mtime or checksum changes.

Testing

def test_projection_known():
    """Catches configuration drift in UTA/SeqRepo/assembly."""
    var = hp.parse_hgvs_variant("NM_000546.6:c.215C>G")
    assert str(am.c_to_g(var)) == "NC_000017.11:g.7676154G>C"

def test_coordinate_conversion():
    assert gff_to_slice(100, 200) == slice(99, 200)

def test_pipeline_smoke(tmp_path):
    out = run_pipeline("tests/data/mini.vcf", tmp_path)
    assert len(out) == 42

A handful of known-answer tests catch the errors that matter most: a changed reference version, a coordinate-convention slip, an unnoticed library behaviour change.

The reproducibility checklist

• Code in git, committed, clean working tree
• Exact package versions in a lockfile
• Genome build with patch level
• Annotation release recorded
• UTA schema and SeqRepo snapshot recorded
• gnomAD / ClinVar / MANE versions recorded
• Input checksums verified
• Random seeds set and recorded
• provenance.json written next to the results
• Figures regenerable from scripts, not hand-edited
• Known-answer tests pass
• The whole thing runs end-to-end from a clean checkout

That last item is the real test, and the one most projects fail.

Common mistakes

• Pinning packages but not reference data.
• "We used hg38" with no patch level.
• Not recording the UTA/SeqRepo snapshot.
• Editing files in data/raw/.
• Hand-editing figures in Illustrator with no record of the edit.
• git rev-parse HEAD from a dirty tree.
• Re-seeding RNGs inside functions.
• Latest-tag base images.
• Mixing pip and conda.
• No end-to-end run from a clean checkout.

See also

Versions and Compatibility · UTA and Data Providers · SeqRepo · NumPy Random Generator · Reading and Writing Data · Saving Figures · rcParams and Style Sheets · Applied - Variant Normalization Pipeline

Exercise

Audit a lockfile against the declared requirements, and produce the digest that lets a later run prove it used the same lock.

Write audit(text), returning a list of output lines. The input holds the requirement section, then a line containing exactly --, then the lock section. A requirement is a distribution name optionally followed by comma-separated constraints, each an operator among ==, !=, >=, <=, > and < and a dotted numeric version, as in numpy>=2.4,<3. A lock line is always name==version. Ignore blank lines and lines beginning with #. Normalise every name by lower-casing it and rewriting underscores as hyphens, so Bio_Python and bio-python are one distribution. Compare versions component by component as integers, padding the shorter one with zeros, so 1.5 and 1.5.0 are equal and 2.10 is above 2.9.

Classify each requirement as satisfied when the locked version meets every one of its constraints, violated when the distribution is locked but fails at least one, and missing when it is not locked at all. Print requirements=<n> locked=<n>, then ok=<n> violated=<n> missing=<n> extra=<n>, where extra counts locked distributions with no requirement. Then one line per violated requirement, sorted by name, as violated<TAB>name<TAB>constraints<TAB>locked version, with the constraints rejoined by commas or * when there were none. Then one line per missing requirement, sorted by name, as missing<TAB>name. Finally print lock_sha256=<digest>, the first twelve hexadecimal characters of the SHA-256 of the normalised lock entries written as name==version, sorted ascending and joined with single newlines.