Applied - Reproducible Environment

In one line: pinning your Python packages is the easy half — in bioinformatics the reference data versions determine your results just as much, and almost nobody records them.

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

SEED = 42
rng = np.random.default_rng(SEED)

def bootstrap(data, n, rng):        # ← pass the rng, never re-seed inside
    ...

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