diff --git a/.github/ISSUE_TEMPLATE/bug_report.yml b/.github/ISSUE_TEMPLATE/bug_report.yml
new file mode 100644
index 0000000..436c662
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/bug_report.yml
@@ -0,0 +1,151 @@
+name: Bug report
+description: A run fails, or BAMpiro produces a result that looks wrong
+title: "[bug] "
+labels: ["bug"]
+body:
+ - type: markdown
+ attributes:
+ value: |
+ Thanks for the report. Most first-run failures are already covered in
+ [Troubleshooting](https://github.com/PathoGenOmics-Lab/BAMpiro/blob/main/docs/troubleshooting.md);
+ a common one is a request that exceeds the machine ("Process requirement exceeds
+ available memory"), which `--max_cpus` and `--max_memory` cap.
+
+ When a task fails, Nextflow prints its **work directory**. `cd` there before filing:
+ `.command.err`, `.command.sh` and `.command.log` are what make the bug reproducible for us.
+
+ - type: textarea
+ id: what-happened
+ attributes:
+ label: What happened
+ description: What you ran into, and at which point of the run.
+ validations:
+ required: true
+
+ - type: textarea
+ id: expected
+ attributes:
+ label: What you expected instead
+ description: For a wrong result rather than a crash, say what the right answer looks like and how you know.
+ validations:
+ required: true
+
+ - type: textarea
+ id: command
+ attributes:
+ label: The exact command
+ description: The full `nextflow run` line, including every `--` parameter and `-` option.
+ render: shell
+ placeholder: |
+ nextflow run main.nf --tsv samples.tsv --outdir results_bampiro -profile local,docker
+ validations:
+ required: true
+
+ - type: dropdown
+ id: profile
+ attributes:
+ label: Profile(s) used with -profile
+ multiple: true
+ options:
+ - standard (SLURM, the default)
+ - local
+ - docker
+ - generic (non-tuberculosis organism)
+ - test / test_full
+ - a custom profile or -c config of my own
+ validations:
+ required: true
+
+ - type: dropdown
+ id: engine
+ attributes:
+ label: Container engine
+ options:
+ - Docker
+ - Singularity
+ - Apptainer
+ - None (-stub-run or tools on PATH)
+ validations:
+ required: true
+
+ - type: dropdown
+ id: where
+ attributes:
+ label: Where it runs
+ options:
+ - Laptop or single workstation
+ - SLURM cluster
+ - Another scheduler
+ - Cloud VM
+ validations:
+ required: true
+
+ - type: input
+ id: nextflow-version
+ attributes:
+ label: Nextflow version
+ description: The output of `nextflow -version` (the version line is enough).
+ placeholder: "24.04.2"
+ validations:
+ required: true
+
+ - type: input
+ id: bampiro-version
+ attributes:
+ label: BAMpiro version
+ description: >-
+ The release tag or `git rev-parse --short HEAD` if you cloned. It is also `manifest.version`
+ in `nextflow.config` and appears in the QC report header.
+ placeholder: "1.0.1"
+ validations:
+ required: true
+
+ - type: input
+ id: failing-process
+ attributes:
+ label: Failing process
+ description: The process name Nextflow reports, e.g. `MAP_READS`, `FREEBAYES_CALL`, `QC_REPORT`.
+ placeholder: "CONSENSUS_FASTA (Consensus: T00001)"
+
+ - type: textarea
+ id: command-err
+ attributes:
+ label: .command.err from the task's work directory
+ description: The end of the file is usually the useful part. Trim it if it is huge.
+ render: shell
+
+ - type: textarea
+ id: command-sh
+ attributes:
+ label: .command.sh from the task's work directory
+ description: The resolved script Nextflow actually ran, with your parameters substituted in.
+ render: shell
+
+ - type: textarea
+ id: nextflow-log
+ attributes:
+ label: Nextflow's error block, or .command.log
+ description: >-
+ The `ERROR ~ Error executing process ...` block Nextflow prints, or the relevant part of
+ `.nextflow.log` / the task's `.command.log`.
+ render: shell
+
+ - type: textarea
+ id: samplesheet
+ attributes:
+ label: The samplesheet row(s) involved
+ description: >-
+ Header plus the affected rows. Replace real paths and identifiers if the cohort is not
+ public. Please do not paste patient data.
+ render: text
+
+ - type: checkboxes
+ id: checks
+ attributes:
+ label: Before submitting
+ options:
+ - label: I read [Troubleshooting](https://github.com/PathoGenOmics-Lab/BAMpiro/blob/main/docs/troubleshooting.md) and my problem is not one of the entries there
+ required: true
+ - label: I searched the existing issues
+ required: true
+ - label: I re-ran with `-resume` and it fails the same way
diff --git a/.github/ISSUE_TEMPLATE/config.yml b/.github/ISSUE_TEMPLATE/config.yml
new file mode 100644
index 0000000..1cf84e7
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/config.yml
@@ -0,0 +1,11 @@
+blank_issues_enabled: false
+contact_links:
+ - name: Troubleshooting guide
+ url: https://github.com/PathoGenOmics-Lab/BAMpiro/blob/main/docs/troubleshooting.md
+ about: Common first-run errors and their fixes. Check this before filing a bug.
+ - name: Documentation
+ url: https://github.com/PathoGenOmics-Lab/BAMpiro/tree/main/docs
+ about: Installation, quick start, the samplesheet format, the full parameter reference, the QC report and the tutorial series.
+ - name: Questions about running BAMpiro
+ url: https://github.com/PathoGenOmics-Lab/BAMpiro/issues
+ about: Discussions are not enabled on this repository, so please open an issue and label it as a question.
diff --git a/.github/ISSUE_TEMPLATE/feature_request.yml b/.github/ISSUE_TEMPLATE/feature_request.yml
new file mode 100644
index 0000000..48288b0
--- /dev/null
+++ b/.github/ISSUE_TEMPLATE/feature_request.yml
@@ -0,0 +1,74 @@
+name: Feature request
+description: Suggest a new option, process, output or report panel
+title: "[feature] "
+labels: ["enhancement"]
+body:
+ - type: markdown
+ attributes:
+ value: |
+ BAMpiro is a small lab tool, so the useful proposals are the ones tied to an analysis you
+ are actually trying to run. Describe that first; the implementation is the easy part to
+ discuss afterwards.
+
+ If the thing you need is already possible through a parameter, the
+ [configuration reference](https://github.com/PathoGenOmics-Lab/BAMpiro/blob/main/docs/configuration.md)
+ lists every one of them.
+
+ - type: textarea
+ id: problem
+ attributes:
+ label: What are you trying to do
+ description: The analysis or workflow that BAMpiro does not currently support, and why the current behaviour gets in the way.
+ validations:
+ required: true
+
+ - type: textarea
+ id: proposal
+ attributes:
+ label: What you would like it to do
+ description: A new parameter, a new process, an extra output file, another report panel. Be as concrete as you can.
+ validations:
+ required: true
+
+ - type: dropdown
+ id: area
+ attributes:
+ label: Which part of the pipeline
+ multiple: true
+ options:
+ - Samplesheet, parameters or profiles
+ - Reference preparation, repeat and mappability masking
+ - Read QC, trimming, Kraken2 contamination screening
+ - Mapping and deduplication
+ - Variant calling and filtering
+ - Annotation, canonical numbering or the coordinate liftover
+ - Consensus sequences
+ - Pathotypr lineage and drug-resistance typing
+ - The interactive QC report
+ - Documentation
+ - Tests or CI
+ - Something else
+ validations:
+ required: true
+
+ - type: input
+ id: organism
+ attributes:
+ label: Organism and reference
+ description: >-
+ BAMpiro's defaults are tuned for MTBC but it is organism-agnostic, and that often decides
+ whether a feature belongs behind a flag.
+ placeholder: "M. tuberculosis / H37Rv, or Legionella pneumophila with my own reference + GFF"
+
+ - type: textarea
+ id: alternatives
+ attributes:
+ label: What you tried instead
+ description: Existing parameters, a custom `-c` config, or post-processing outside the pipeline, and why that was not enough.
+
+ - type: checkboxes
+ id: contribute
+ attributes:
+ label: Would you like to implement it
+ options:
+ - label: I would be willing to open a pull request for this (see [CONTRIBUTING.md](https://github.com/PathoGenOmics-Lab/BAMpiro/blob/main/CONTRIBUTING.md))
diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md
new file mode 100644
index 0000000..c70ca34
--- /dev/null
+++ b/.github/PULL_REQUEST_TEMPLATE.md
@@ -0,0 +1,22 @@
+## What changed and why
+
+
+
+## How it was verified
+
+
+
+```bash
+tests/run_tests.sh
+```
+
+## Checklist
+
+- [ ] `tests/run_tests.sh` passes locally (or the legs I could run do, and I said which)
+- [ ] `ruff check .` is clean, and I did not restyle code the change does not touch
+- [ ] New or changed behaviour has a test, and tests pin what the code does today
+- [ ] A new process in `modules/` has a `stub:` block, and `-profile test_full` reaches it
+- [ ] `tests/data/` was not edited by hand (fixtures come from `make_test_data.py`)
+- [ ] User-visible changes have a `CHANGELOG.md` entry and, where relevant, a docs update
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 0000000..4372e37
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,106 @@
+name: CI
+
+# The gate for every change: linting, unit tests for the Python under bin/, tests for the report's
+# hand-written front-end statistics, and a full stub run of the DAG. (The docs build lives in
+# docs.yml, which runs on pull requests too.)
+#
+# Nothing here needs a container, a reference genome or real reads: the fixture cohort under
+# tests/data/ is 170 kB and `-stub-run` executes each process's stub: block instead of the tool.
+
+on:
+ pull_request:
+ push:
+ branches: [main]
+ workflow_dispatch:
+
+# A new push supersedes the previous run of the same branch. The repository is private on a free
+# plan, so Actions minutes are metered; do not burn them on results nobody will read.
+concurrency:
+ group: ci-${{ github.ref }}
+ cancel-in-progress: true
+
+permissions:
+ contents: read
+
+jobs:
+ lint:
+ name: Lint (ruff)
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+ - uses: actions/setup-python@v5
+ with:
+ python-version: "3.12"
+ - run: pip install ruff
+ - run: ruff check .
+
+ unit:
+ name: Unit tests (Python ${{ matrix.python }})
+ runs-on: ubuntu-latest
+ strategy:
+ fail-fast: false
+ matrix:
+ # The oldest version the scripts must still run on, and the one in the container.
+ python: ["3.10", "3.13"]
+ steps:
+ - uses: actions/checkout@v4
+ - uses: actions/setup-python@v5
+ with:
+ python-version: ${{ matrix.python }}
+ - run: pip install pytest numpy
+ - run: python -m pytest tests/unit -q
+
+ fixtures:
+ name: Fixtures are reproducible
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+ - uses: actions/setup-python@v5
+ with:
+ python-version: "3.12"
+ # The committed cohort must be exactly what the generator produces, or the seed has drifted
+ # and every test that reads tests/data/ is measuring something nobody can regenerate.
+ - name: Regenerate and diff
+ run: |
+ python tests/data/make_test_data.py
+ git diff --exit-code --stat tests/data/
+
+ report-frontend:
+ name: Report front-end (Node)
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+ - uses: actions/setup-node@v4
+ with:
+ node-version: "20"
+ # No dependencies: the tests use node:test and load the report bundle through node:vm.
+ - run: node --test tests/js/*.test.mjs
+
+ pipeline:
+ name: Stub run (Nextflow ${{ matrix.nxf_ver }})
+ runs-on: ubuntu-latest
+ strategy:
+ fail-fast: false
+ matrix:
+ # The minimum the manifest declares (v1 parser) and the current release (strict v2
+ # parser). Both must pass: main.nf is written to satisfy the stricter of the two.
+ nxf_ver: ["24.04.2", "latest-stable"]
+ env:
+ NXF_ANSI_LOG: "false"
+ steps:
+ - uses: actions/checkout@v4
+ - uses: actions/setup-java@v4
+ with:
+ distribution: temurin
+ java-version: "17"
+ - uses: actions/setup-python@v5
+ with:
+ python-version: "3.12"
+ - run: pip install pytest
+ - uses: nf-core/setup-nextflow@v2
+ with:
+ version: ${{ matrix.nxf_ver }}
+ - run: python -m pytest tests/pipeline -q
+
+# The docs are validated by docs.yml, which already installs the cairo/pango libraries the
+# social plugin needs and now runs on pull requests too. Nothing to duplicate here.
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 76e0343..25efb2b 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -15,7 +15,14 @@ name: docs
on:
push:
- branches: [main, indel-mask]
+ branches: [main]
+ paths:
+ - docs/**
+ - mkdocs.yml
+ - .github/workflows/docs.yml
+ # Validate the site on pull requests too, so a broken link fails review rather than main.
+ # Only `build` runs there: `deploy` is gated on a push below.
+ pull_request:
paths:
- docs/**
- mkdocs.yml
diff --git a/.gitignore b/.gitignore
index 02728c6..91baf73 100644
--- a/.gitignore
+++ b/.gitignore
@@ -15,6 +15,11 @@ dag*.dot
results_bampiro/
mappability_cache/
+# Test-suite scratch (tests/run_tests.sh clears these)
+results_test/
+.test_cache/
+.pytest_cache/
+
# Python
__pycache__/
*.pyc
diff --git a/.zenodo.json b/.zenodo.json
new file mode 100644
index 0000000..961c7f0
--- /dev/null
+++ b/.zenodo.json
@@ -0,0 +1,34 @@
+{
+ "title": "BAMpiro: bacterial short-read mapping, variant calling and lineage/drug-resistance typing",
+ "description": "
A Nextflow (DSL2) pipeline for short-read bacterial and Mycobacterium tuberculosis complex analysis: read QC and contamination screening, length-aware mappability masking, mapping, variant calling, consensus generation, SnpEff annotation, alignment-free lineage and WHO drug-resistance typing with Pathotypr, and a single self-contained interactive QC report.
",
+ "upload_type": "software",
+ "license": "GPL-3.0-or-later",
+ "access_right": "open",
+ "creators": [
+ {
+ "name": "Ruiz-Rodriguez, Paula",
+ "affiliation": "I2SysBio (University of Valencia - CSIC)"
+ },
+ {
+ "name": "Coscolla, Mireia",
+ "affiliation": "I2SysBio (University of Valencia - CSIC)"
+ }
+ ],
+ "keywords": [
+ "bioinformatics",
+ "Nextflow",
+ "workflow",
+ "variant calling",
+ "Mycobacterium tuberculosis",
+ "drug resistance",
+ "lineage typing",
+ "quality control"
+ ],
+ "related_identifiers": [
+ {
+ "identifier": "https://github.com/PathoGenOmics-Lab/BAMpiro",
+ "relation": "isSupplementTo",
+ "scheme": "url"
+ }
+ ]
+}
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 045b4a5..b5ba47e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,125 @@ All notable changes to this project are documented in this file. The format is
based on [Keep a Changelog](https://keepachangelog.com/), and the project follows
[Semantic Versioning](https://semver.org/).
+## [Unreleased] - targeting 1.1.0
+
+Headline: a **test suite and CI**, and a pipeline that runs correctly on a machine
+that is not the authors' cluster.
+
+`manifest.version` already reads `1.1.0`: that is the number this work is heading
+for, not a release that happened. Nothing is tagged and no release is published.
+See [RELEASING.md](RELEASING.md) for the order the actual release has to follow.
+
+The **container stays at the 1.0.1 image**, pinned by its digest. Nothing here
+changes a bundled tool, and the scripts added under `bin/` are staged by Nextflow
+from the pipeline directory rather than baked into the image, so there is nothing
+to rebuild. A release only needs a new image when the Dockerfile changes.
+
+### Added
+
+- **Test suite** (`tests/`) - around 1,100 tests in three legs, all runnable with
+ `tests/run_tests.sh` and none of them needing a container, a reference genome or
+ a network connection:
+ - `tests/unit/` - the Python under `bin/`: the consensus decision tree, the QC
+ verdict engine, every parser, the k-mer liftover and the read filter.
+ - `tests/js/` - the report's hand-written ES5 statistics (Spearman,
+ Kruskal-Wallis, chi-square, Benjamini-Hochberg), checked against **SciPy and
+ statsmodels** reference values rather than a snapshot of our own output, plus
+ the integrity of the asset bundle.
+ - `tests/pipeline/` - samplesheet validation and a full `-stub-run` of the DAG.
+- **A `stub:` block for every process**, so `-stub-run` walks the whole workflow in
+ seconds with no data. A test fails if a new process arrives without one.
+- **Fixture cohort** (`tests/data/`, 170 kB) - a 5 kb reference with three genes and
+ three samples covering paired-end, single-end and a sample split over two runs,
+ all derived from one seed. CI regenerates it and fails on any diff.
+- **`test` and `test_full` profiles** - the default feature set, and every optional
+ branch on so a single run instantiates every process.
+- **Continuous integration** (`.github/workflows/ci.yml`) - lint, unit tests on two
+ Python versions, the front-end tests, and a stub run against both the minimum
+ supported Nextflow and the current release, on every pull request.
+- **`slurm` and `garnatxa` execution profiles**. `-profile garnatxa` replaces the
+ hand-written `-c cluster.config` a user previously needed, and works when the
+ pipeline is pulled straight from GitHub.
+- **`--help`**, listing every parameter with its default.
+- **`--max_cpus` / `--max_memory` / `--max_time`** - cap every process to what the
+ machine can give, so an over-sized request (Kraken2 asks for 80 GB) fits a laptop
+ without editing the modules.
+- **Community files** - CONTRIBUTING, a Code of Conduct, a security policy, issue
+ forms, a pull-request template, `RELEASING.md` and `.zenodo.json`.
+
+### Changed
+
+- **`-profile standard` no longer submits to SLURM.** It is now the portable default
+ and runs on the current host, so a fresh clone works anywhere. Use `-profile slurm`
+ for a generic cluster or `-profile garnatxa` for the I2SysBio one. **This changes
+ what an existing `-profile standard` command does.**
+- The base config no longer carries one site's paths: the SLURM QoS, the
+ `module load singularity`, the Singularity `--bind` paths and the Kraken database
+ location all moved into `conf/garnatxa.config`.
+- **`--tsv` is required** and says what to pass; the default was a leftover
+ samplesheet name. `--kraken2_db` has no default either, and without one the
+ contamination screen is skipped.
+- **The container is pinned by digest** rather than the mutable `1.0.1` tag, so two
+ runs a month apart cannot silently use different tool versions.
+- **An unrecognised `--parameter` now stops the run** with a suggestion, instead of
+ being accepted and ignored. The declared set is read from `nextflow.config`, so
+ the two cannot drift.
+- The startup banner reports the resolved profile, executor, container and Kraken
+ state, because where a run lands is the easiest thing to get wrong.
+
+### Fixed
+
+- **The QC report could not be produced on the default settings.** Its three optional
+ inputs all fall back to a placeholder, and all three features are off by default,
+ so every run handed `QC_REPORT` the same `NO_FILE` path three times and Nextflow
+ refused to stage them (`input file name collision`). The placeholders are now
+ distinct files under `assets/`.
+- **`main.nf` did not compile under the strict Nextflow parser**, the default from
+ 25.10, and only ran through the deprecated v1 fallback. The samplesheet parsing is
+ now a function, every statement lives inside the workflow, and each dynamic
+ `publishDir` is a closure. Verified on 24.04.2 and 26.04.6.
+- Removed `process.publishDirMode`, which is not a Nextflow directive and only
+ produced a warning on every run.
+- A blank line in a bedgraph raised `ValueError` in `build_min_unique_len.py`
+ instead of being skipped.
+- **A GFF attribute was matched as a substring**, so `locus_tag=` also matched
+ inside `old_locus_tag=` and `gene=` inside `pseudogene=`, both routine in RefSeq
+ and Prokka output. Genes were labelled with an obsolete tag, and the curated
+ H37Rv gene map could be silently overridden.
+- `mask_profile` summed raw interval lengths for its headline total, so an
+ unmerged BED could report above 100% masked.
+- Four parsers in `stats_to_legacy.py` wrapped a whole read loop in one
+ `try/except`, so one malformed line silently discarded every line after it; a
+ partial genome size then propagated into depth, breadth, evenness and every
+ density bin.
+- Genuinely-zero metrics were written as `NA`, making a failed sample
+ indistinguishable from an unmeasured one.
+- Genotype classification assumed diploid spellings, so a haploid no-call was
+ counted as heterozygous. This matters for `-profile generic`, which is haploid.
+- `safe_tabix` had drifted between its copies: the one in `annotation.nf` still
+ used a `zgrep` check that cannot tell a truncated VCF from an empty one, and
+ wrote a fake index without saying so.
+- A missing `FORMAT/DP` reached the consensus as `ADP=.`, which reads as zero
+ depth and turns a called SNP into a gap.
+
+### Refactored
+
+Behaviour-preserving throughout, and each step verified rather than assumed:
+
+- **The science came out of the process scripts.** Around 500 lines of bash and
+ awk lived inside `modules/*.nf`, where no test could reach it: a stub run
+ replaces the script block and nothing can import it. The genotype
+ re-validation, the backbone pileup parser, the hom/het bcftools expressions and
+ `safe_tabix` now live in `bin/`, with 62 tests over them, including a run
+ against a real bcftools. `variants.nf` drops from 485 lines to 344.
+- **`qc_report.py` became a package.** 1697 lines and 53 functions are now a
+ 319-line CLI over `bin/qcreport/` (parsers, metrics, panels, render). The old
+ and new versions produce byte-identical output on the demo cohort.
+- **The workflow is one sub-workflow per stage.** `main.nf` goes from 556 lines to
+ 277, with the stage sequence reduced to eleven calls. Equivalence checked by
+ diffing execution traces against the previous version: identical task counts and
+ identical published outputs on both test profiles.
+
## [1.0.1] - 2026-07-19
Headline: a consolidated **interactive QC report**, alignment-free **lineage &
diff --git a/CITATION.cff b/CITATION.cff
index 7dc3b2f..9c8a4d6 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -16,7 +16,7 @@ authors:
repository-code: "https://github.com/PathoGenOmics-Lab/BAMpiro"
url: "https://github.com/PathoGenOmics-Lab/BAMpiro"
license: GPL-3.0
-version: 1.0.1
+version: 1.1.0
# Add `doi:` here once the release is archived on Zenodo (mints a persistent software DOI).
keywords:
- bioinformatics
diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md
new file mode 100644
index 0000000..239bba4
--- /dev/null
+++ b/CODE_OF_CONDUCT.md
@@ -0,0 +1,83 @@
+# Contributor Covenant Code of Conduct
+
+## Our Pledge
+
+We as members, contributors, and leaders pledge to make participation in our community a harassment-free experience for everyone, regardless of age, body size, visible or invisible disability, ethnicity, sex characteristics, gender identity and expression, level of experience, education, socio-economic status, nationality, personal appearance, race, caste, color, religion, or sexual identity and orientation.
+
+We pledge to act and interact in ways that contribute to an open, welcoming, diverse, inclusive, and healthy community.
+
+## Our Standards
+
+Examples of behavior that contributes to a positive environment for our community include:
+
+* Demonstrating empathy and kindness toward other people
+* Being respectful of differing opinions, viewpoints, and experiences
+* Giving and gracefully accepting constructive feedback
+* Accepting responsibility and apologizing to those affected by our mistakes, and learning from the experience
+* Focusing on what is best not just for us as individuals, but for the overall community
+
+Examples of unacceptable behavior include:
+
+* The use of sexualized language or imagery, and sexual attention or advances of any kind
+* Trolling, insulting or derogatory comments, and personal or political attacks
+* Public or private harassment
+* Publishing others' private information, such as a physical or email address, without their explicit permission
+* Other conduct which could reasonably be considered inappropriate in a professional setting
+
+## Enforcement Responsibilities
+
+Community leaders are responsible for clarifying and enforcing our standards of acceptable behavior and will take appropriate and fair corrective action in response to any behavior that they deem inappropriate, threatening, offensive, or harmful.
+
+Community leaders have the right and responsibility to remove, edit, or reject comments, commits, code, wiki edits, issues, and other contributions that are not aligned to this Code of Conduct, and will communicate reasons for moderation decisions when appropriate.
+
+## Scope
+
+This Code of Conduct applies within all community spaces, and also applies when an individual is officially representing the community in public spaces. Examples of representing our community include using an official e-mail address, posting via an official social media account, or acting as an appointed representative at an online or offline event.
+
+## Enforcement
+
+Instances of abusive, harassing, or otherwise unacceptable behavior may be reported to the community leaders responsible for enforcement at paula.ruiz.rodriguez@csic.es. All complaints will be reviewed and investigated promptly and fairly.
+
+All community leaders are obligated to respect the privacy and security of the reporter of any incident.
+
+## Enforcement Guidelines
+
+Community leaders will follow these Community Impact Guidelines in determining the consequences for any action they deem in violation of this Code of Conduct:
+
+### 1. Correction
+
+**Community Impact**: Use of inappropriate language or other behavior deemed unprofessional or unwelcome in the community.
+
+**Consequence**: A private, written warning from community leaders, providing clarity around the nature of the violation and an explanation of why the behavior was inappropriate. A public apology may be requested.
+
+### 2. Warning
+
+**Community Impact**: A violation through a single incident or series of actions.
+
+**Consequence**: A warning with consequences for continued behavior. No interaction with the people involved, including unsolicited interaction with those enforcing the Code of Conduct, for a specified period of time. This includes avoiding interactions in community spaces as well as external channels like social media. Violating these terms may lead to a temporary or permanent ban.
+
+### 3. Temporary Ban
+
+**Community Impact**: A serious violation of community standards, including sustained inappropriate behavior.
+
+**Consequence**: A temporary ban from any sort of interaction or public communication with the community for a specified period of time. No public or private interaction with the people involved, including unsolicited interaction with those enforcing the Code of Conduct, is allowed during this period. Violating these terms may lead to a permanent ban.
+
+### 4. Permanent Ban
+
+**Community Impact**: Demonstrating a pattern of violation of community standards, including sustained inappropriate behavior, harassment of an individual, or aggression toward or disparagement of classes of individuals.
+
+**Consequence**: A permanent ban from any sort of public interaction within the community.
+
+## Attribution
+
+This Code of Conduct is adapted from the [Contributor Covenant][homepage], version 2.1, available at [https://www.contributor-covenant.org/version/2/1/code_of_conduct.html][v2.1].
+
+Community Impact Guidelines were inspired by [Mozilla's code of conduct enforcement ladder][Mozilla CoC].
+
+For answers to common questions about this code of conduct, see the FAQ at [https://www.contributor-covenant.org/faq][FAQ]. Translations are available at [https://www.contributor-covenant.org/translations][translations].
+
+[homepage]: https://www.contributor-covenant.org
+[v2.1]: https://www.contributor-covenant.org/version/2/1/code_of_conduct.html
+[Mozilla CoC]: https://github.com/mozilla/diversity
+[FAQ]: https://www.contributor-covenant.org/faq
+[translations]: https://www.contributor-covenant.org/translations
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
new file mode 100644
index 0000000..d4f9e14
--- /dev/null
+++ b/CONTRIBUTING.md
@@ -0,0 +1,250 @@
+# Contributing to BAMpiro
+
+BAMpiro is a Nextflow (DSL2) pipeline maintained by a small group at I2SysBio (University of
+Valencia - CSIC) as part of the [PathoGenOmics-Lab](https://github.com/PathoGenOmics-Lab). Bug
+reports, questions and patches are all welcome, and none of them need to be large: a corrected
+default, a clearer error message or a missing line in the docs is a useful contribution.
+
+This document covers how to get a working setup, how to run the tests, what CI will check, and the
+few local conventions that are not obvious from reading the code.
+
+Everything here happens under the [Code of Conduct](CODE_OF_CONDUCT.md). Contributions are accepted
+under the project's [GPL-3.0 license](LICENSE).
+
+## Getting set up
+
+You need two things to run the pipeline itself:
+
+- **Nextflow** `>= 24.04.2`, which needs **Java 17+**
+- **Docker** or **Singularity / Apptainer** for the container
+
+The cleanest way to get both without touching your system Java is a conda environment, which is what
+[Your first run](docs/tutorials/first-run.md) recommends:
+
+```bash
+conda create -n bampiro -c conda-forge -c bioconda nextflow
+conda activate bampiro
+nextflow -version # confirm it is >= 24.04.2
+```
+
+Then clone the repository and check that it runs:
+
+```bash
+git clone https://github.com/PathoGenOmics-Lab/BAMpiro
+cd BAMpiro
+nextflow run main.nf -profile test -stub-run
+```
+
+That last command executes the whole DAG against the bundled 170 kB fixture cohort in a couple of
+seconds. It pulls no container, reads no reference genome and touches the network only if Nextflow
+needs to fetch a plugin.
+
+For the test suite and the linter you also want:
+
+- **Python** `>= 3.10` with `pytest` and `numpy` (`pip install pytest numpy ruff`)
+- **Node** `>= 18` for the report front-end tests
+
+Working on the docs additionally needs MkDocs Material: `pip install -r docs/requirements.txt`, then
+`make docs` (a strict build) or `make docs-serve` (live reload on http://127.0.0.1:8000).
+
+## Running the tests
+
+One script runs everything, the same way CI does:
+
+```bash
+tests/run_tests.sh # everything
+tests/run_tests.sh lint # ruff only
+tests/run_tests.sh unit # Python unit tests for bin/ only
+tests/run_tests.sh js # report front-end only
+tests/run_tests.sh pipeline # Nextflow stub runs only
+```
+
+There are around 1,200 tests in three legs:
+
+| Leg | What it covers | Needs |
+| :--- | :--- | :--- |
+| `tests/unit/` | Everything under `bin/`: the consensus decision tree, the QC verdict engine, the parsers, the k-mer liftover, the read filter, and the awk and filter expressions the process scripts call | `pytest`, `numpy` |
+| `tests/js/` | The report's hand-written ES5 statistics, checked against SciPy and statsmodels reference values, plus the integrity of the asset bundle | `node >= 18` |
+| `tests/pipeline/` | Samplesheet validation and a full `-stub-run` of the DAG on both test profiles | `nextflow`, Java 17+ |
+
+A leg whose tool is not installed skips itself instead of failing, so a partial local setup still
+gives useful signal. Nothing needs a container, a reference genome, real reads or a network
+connection. [`tests/README.md`](tests/README.md) documents the fixture cohort, the two test profiles
+and the Node harness in more detail.
+
+Two things to know before you add a test:
+
+- Anything written to disk goes in `tmp_path`. Do not add files to `tests/data/` by hand. If a
+ fixture is worth keeping, add it to `tests/data/make_test_data.py`, which derives the whole cohort
+ from one seed. CI regenerates the cohort and fails on any diff, so a hand-added file breaks the
+ build.
+- Unit tests import a script from `bin/` with `load_script("name")` from `tests/conftest.py`; `bin/`
+ is not a package. Importing `qc_report` reads `bin/report_assets/`, so it only works in place.
+
+## What CI runs
+
+[`.github/workflows/ci.yml`](.github/workflows/ci.yml) gates every pull request with five jobs:
+
+| Job | What it does |
+| :--- | :--- |
+| Lint (ruff) | `ruff check .` |
+| Unit tests | `pytest tests/unit` on Python 3.10 and 3.13 |
+| Fixtures are reproducible | Regenerates `tests/data/` and fails on any diff |
+| Report front-end (Node) | `node --test tests/js/*.test.mjs` on Node 20 |
+| Stub run | `pytest tests/pipeline` against Nextflow `24.04.2` and `latest-stable` |
+
+The two Python versions are the oldest the scripts must still run on and the one in the container.
+The two Nextflow versions are the minimum the manifest declares and the current release, the latter
+as an early warning rather than a hard requirement.
+
+[`.github/workflows/docs.yml`](.github/workflows/docs.yml) additionally runs `mkdocs build --strict`
+on pull requests that touch `docs/`, `mkdocs.yml` or the workflow itself, so a broken link or a page
+missing from the nav fails review rather than `main`.
+
+Actions minutes are metered on this repository, so a new push cancels the previous run of the same
+branch. Run `tests/run_tests.sh` locally before you push rather than using CI as your test runner.
+
+## Linting: narrow on purpose
+
+`ruff` is configured in [`pyproject.toml`](pyproject.toml) with a deliberately narrow rule set:
+
+```toml
+select = ["F", "E9", "B006", "B008", "B017", "B904"]
+line-length = 120
+```
+
+That catches real defects (undefined names, unused imports and variables, broken f-strings, syntax
+errors, mutable default arguments) and nothing else. The Python under `bin/` chose a compact, dense
+layout, and the rule set is narrow so that choice is not relitigated file by file.
+
+**Do not restyle the codebase.** No repository-wide formatter run, no import reordering, no
+reflowing of code you are not otherwise changing. A diff that mixes a one-line fix with 400 lines of
+reformatting cannot be reviewed. If you think a rule should be added, propose it in an issue on its
+own, separately from the change that prompted it.
+
+`bin/extract_kraken_reads.py` is excluded from linting entirely. See
+[Two things worth knowing](#two-things-worth-knowing) below.
+
+## Commit messages
+
+The convention in this repository is a lowercase area prefix, a colon, and a short lowercase summary
+of what the change does, followed by a blank line and two to four sentences of prose explaining it.
+No bullet lists, no headed sections, no trailers. Read `git log --oneline -30` for the shape of it.
+
+The prefix is usually a conventional-commit type (`fix:`, `test:`, `docs:`, `ci:`, `refactor:`), but
+the part of the pipeline you touched works just as well and is used often (`report:`, `liftover:`,
+`mask:`, `pathotypr:`, `config:`, `container:`). A parenthesised scope is fine where it helps:
+`docs(tutorials):`, `ci(docs):`.
+
+The body says why the change is the way it is, not what the diff already shows:
+
+```
+ci: validate the docs on pull requests instead of duplicating the build
+
+docs.yml already runs mkdocs build --strict with the system libraries the social
+plugin needs, so it just gains a pull_request trigger; only the build job runs
+there, since deploy stays gated on the Pages variable. Drops the stale
+indel-mask branch trigger now that the branch is merged.
+```
+
+## Where code goes
+
+The pipeline is layered, and the layer decides whether your change is testable:
+
+| Layer | What belongs there |
+| :--- | :--- |
+| `main.nf` | Reading the samplesheet and calling each stage. Nothing else. |
+| `subworkflows/*.nf` | One per pipeline stage, with an explicit `take:` / `emit:`. All the channel plumbing, and the `if (params.x)` that switches a stage off. |
+| `modules/*.nf` | The processes: their inputs, outputs, resources and the command line they run. |
+| `bin/` | Everything the command line actually does, once it is more than a call to a tool. |
+
+That last row is the one worth internalising. **A process script cannot be tested.** `-stub-run`
+replaces it wholesale, and nothing can import it, so a decision written inside a `script:` block has
+no coverage and cannot get any. If your change involves a threshold, a filter expression, a parse or
+any other decision, put it in a file under `bin/` and have the process call it. That is why
+`vcf_filter_rules.py`, `format_snps_for_backbone.awk` and `backbone_allpos.awk` exist.
+
+A sub-workflow's `emit:` name must not match a `def` local in the same scope: `emit: x = y` is a
+real assignment, so a local `x` shadows the binding property Nextflow then looks for, and the
+failure is confusing. Give the emit a different name.
+
+## Adding a process
+
+Every process in `modules/` needs a `stub:` block. A stub `touch`es each declared output and nothing
+else:
+
+```groovy
+stub:
+"""
+touch ${sampleId}.${refId}.consensus.fasta
+touch ${sampleId}.${refId}.consensus.log
+"""
+```
+
+Without one, `-stub-run` aborts and the whole pipeline leg of CI goes red.
+`tests/pipeline/test_stub_run.py::test_there_is_a_stub_block_for_every_process` also checks this
+statically, so you get a clear failure instead of a confusing Nextflow error.
+
+If the process sits on an optional branch, make sure `-profile test_full` reaches it. That profile
+turns on every optional feature (Kraken, pathotypr typing, canonical annotation, the liftover, the
+virgin consensus) precisely so one stub run instantiates every declared process, and
+`test_the_full_profile_reaches_every_process` fails when something drops out of that coverage.
+
+## Tests pin observed behaviour
+
+The suite records what the code does **today**, not what it ought to do. When a test looks like it is
+asserting something wrong, that is usually deliberate and there is a comment saying so.
+
+If you find a genuine defect, do not leave a failing test behind to mark it. Split it in two:
+
+1. A test that pins the current, wrong behaviour, with a comment explaining why it is wrong.
+2. A separate change that fixes the code and flips the test.
+
+Both can be in the same pull request, as separate commits. What matters is that `main` is never red:
+a suite that is expected to fail is a suite everyone learns to ignore.
+
+## Changelog
+
+[`CHANGELOG.md`](CHANGELOG.md) follows [Keep a Changelog](https://keepachangelog.com/) and the
+project follows [Semantic Versioning](https://semver.org/). Anything a user would notice (a new
+parameter, a changed default, a fixed wrong result, a new report panel) gets an entry under
+`Added`, `Changed` or `Fixed` in the section for the version being prepared. Add an `## [Unreleased]`
+section at the top if there is not one yet.
+
+Internal work with no user-visible effect (refactors, tests, CI, most docs changes) does not need an
+entry.
+
+## Two things worth knowing
+
+**`main.nf` has to satisfy the strict Nextflow parser.** From 25.10 that parser is the default, and
+it is stricter than the language you may be used to: no statements at the top level of a script
+(everything lives in a `workflow` or a function), no `while` loops, no C-style `for`, no assignment
+used as an expression, and a dynamic process directive must be a closure -
+`publishDir path: { "..." }`, not `publishDir "..."`. If a run dies with
+`Statements cannot be mixed with script declarations`, that is what you have hit. CI runs both the
+minimum supported 24.04.2 and `latest-stable`, so either parser rejecting your change fails the
+build.
+
+**`bin/extract_kraken_reads.py` is vendored and must stay byte-identical to upstream.** It comes
+from [KrakenTools](https://github.com/jenniferlu717/KrakenTools) (the URL is `params.krakentools_url`
+in `nextflow.config`, and `modules/qc.nf` falls back to downloading it when the local copy is
+missing). It is not ours to restyle, refactor or "fix", and `pyproject.toml` excludes it from ruff
+for that reason. Updating it means replacing it with a newer upstream copy, in its own commit.
+
+## Opening a pull request
+
+Branch off `main`, keep the change focused, and fill in
+[the pull request template](.github/PULL_REQUEST_TEMPLATE.md). Run `tests/run_tests.sh` before you
+push. If you are unsure whether an idea fits, open an issue first and ask; that is cheaper than
+writing code that gets turned down.
+
+## Reporting bugs and asking questions
+
+Use the [issue templates](https://github.com/PathoGenOmics-Lab/BAMpiro/issues/new/choose). The bug
+template asks for the things that actually make a BAMpiro failure reproducible: your Nextflow
+version, the exact command, the profile, the container engine, the failing process name and the
+contents of its work directory (`.command.err`, `.command.sh`, `.command.log`). Nextflow prints the
+work directory of a failed task, so `cd` there first.
+
+Before you file, check [Troubleshooting](docs/troubleshooting.md). Most first-run failures are one of
+the entries in that table.
diff --git a/README.md b/README.md
index 37b8105..1f38089 100644
--- a/README.md
+++ b/README.md
@@ -4,9 +4,10 @@
+[](https://github.com/PathoGenOmics-Lab/BAMpiro/actions/workflows/ci.yml)
[](LICENSE)
[](https://www.nextflow.io/)
-[](CHANGELOG.md)
+[](CHANGELOG.md)
[](https://hub.docker.com/r/paururo/bampiro)
[](https://github.com/PathoGenOmics-Lab)
[](https://pathogenomics-lab.github.io/BAMpiro/)
@@ -44,7 +45,7 @@ tuberculosis* but is **organism-agnostic** - point it at any reference genome +
- Alignment-free MTBC lineage + WHO drug-resistance typing ([Pathotypr](docs/pathotypr.md))
- A self-contained, interactive [HTML QC report](docs/qc-report.md) with 21 panels
- Dual amino-acid numbering (used reference + H37Rv / Mycobrowser)
-- One pinned container with every tool and marker panel built in
+- One digest-pinned container with every tool and marker panel built in
## Features
@@ -56,19 +57,25 @@ tuberculosis* but is **organism-agnostic** - point it at any reference genome +
| 🩺 Lineage & drug resistance | Alignment-free MTBC lineage + WHO DR typing ([Pathotypr](docs/pathotypr.md)), reference-agnostic |
| 📊 Interactive QC report | Self-contained HTML dashboard, [21 panels](docs/qc-report.md) + per-sample `qc_flags.tsv` |
| 🔤 Dual amino-acid numbering | Protein changes in both the used reference and H37Rv/Mycobrowser numbering |
-| 📦 Containerized & reproducible | A single pinned image with every tool + bundled marker panels ([details](docs/installation.md)) |
+| 📦 Containerized & reproducible | A single image, pinned by digest, with every tool + bundled marker panels ([details](docs/installation.md)) |
| ⚙️ Fully configurable | Every step exposed as a Nextflow parameter ([reference](docs/configuration.md)) |
## Installation
Requires **Nextflow ≥ 24.04.2** and **Docker** or **Singularity**. The pipeline
-pulls a pinned `paururo/bampiro` image with every tool built in - nothing else to
-install.
+pulls the `paururo/bampiro` image, pinned by digest, with every tool built in - nothing
+else to install.
```bash
-nextflow run main.nf --tsv samples.tsv --outdir results_bampiro -profile standard
+nextflow run main.nf --tsv samples.tsv --outdir results_bampiro -profile local,docker
```
+`--tsv` is required, and `-profile` says where the work runs: `local` on this machine,
+`slurm` on any SLURM cluster, `garnatxa` on the I2SysBio one. With no `-profile` everything
+runs on the current host, which on a cluster login node means the login node - so choose
+one deliberately there. `nextflow run main.nf --help` lists every parameter with its
+default.
+
Full requirements and the bundled software versions: [Installation](docs/installation.md).
## Quick Start
@@ -79,7 +86,7 @@ amino-acid numbering) with:
```bash
nextflow run main.nf \
- --tsv samples.tsv --outdir results_bampiro -profile standard \
+ --tsv samples.tsv --outdir results_bampiro -profile local,docker \
--run_pathotypr true --annotate_canonical true
```
@@ -88,7 +95,7 @@ See the [Quick Start guide](docs/quickstart.md) for the full format and multi-ru
merging, and [Outputs](docs/outputs.md) for the result layout.
**New to BAMpiro?** Follow the hands-on
-[**Jupyter tutorial**](docs/tutorial/bampiro_tutorial.ipynb) — install → samplesheet →
+[**Jupyter tutorial**](docs/tutorial/bampiro_tutorial.ipynb) - install → samplesheet →
run → explore the outputs with `pandas` / `matplotlib` on a bundled 17-sample example
cohort (runnable without running the pipeline first).
@@ -96,13 +103,13 @@ cohort (runnable without running the pipeline first).
-
+
▶ Open the live interactive demo report
- — a full 17-sample demo cohort, every panel populated, right in your browser
+ - a full 17-sample demo cohort, every panel populated, right in your browser
Every run writes a single self-contained `
_qc_report.html` (no internet,
@@ -118,13 +125,13 @@ Full panel list and interactive features: [Interactive QC Report](https://pathog
## Documentation
-📖 **[Read the full documentation online → pathogenomics-lab.github.io/BAMpiro](https://pathogenomics-lab.github.io/BAMpiro/)**
-— a searchable site with a pipeline diagram, a hover glossary, a step-by-step
+📖 **[Read the full documentation online → pathogenomics-lab.github.io/BAMpiro](https://pathogenomics-lab.github.io/BAMpiro/)**:
+a searchable site with a pipeline diagram, a hover glossary, a step-by-step
[**tutorial series**](https://pathogenomics-lab.github.io/BAMpiro/tutorials/), a hands-on
Jupyter notebook, and the [live QC-report demo](https://pathogenomics-lab.github.io/BAMpiro/examples/qc_report_demo.html).
> [!TIP]
-> **New here?** Follow the [**Tutorials**](https://pathogenomics-lab.github.io/BAMpiro/tutorials/) —
+> **New here?** Follow the [**Tutorials**](https://pathogenomics-lab.github.io/BAMpiro/tutorials/):
> a guided path from your first run to a finished phylogeny.
| Document | Description |
@@ -146,6 +153,17 @@ locally with `make docs-serve`
(needs [MkDocs Material](https://squidfunk.github.io/mkdocs-material/):
`pip install -r docs/requirements.txt`).
+## Tests
+
+```bash
+tests/run_tests.sh
+```
+
+Unit tests for the Python under `bin/`, tests for the report front-end's hand-written statistics
+(checked against SciPy), and a full `-stub-run` of the DAG over a 170 kB fixture cohort. No
+container, no reference genome, no network. [`tests/README.md`](tests/README.md) has the details;
+CI runs the same checks on every pull request.
+
## Why "BAMpiro"?
The name is a play on words combining bioinformatics and folklore:
@@ -172,7 +190,7 @@ If you use BAMpiro in your research, please cite:
title = {BAMpiro: bacterial short-read mapping, variant calling and lineage/drug-resistance typing},
author = {Ruiz-Rodriguez, Paula and Coscoll{\'a}, Mireia},
url = {https://github.com/PathoGenOmics-Lab/BAMpiro},
- version = {1.0.1},
+ version = {1.1.0},
license = {GPL-3.0}
}
```
diff --git a/RELEASING.md b/RELEASING.md
new file mode 100644
index 0000000..8bb2d92
--- /dev/null
+++ b/RELEASING.md
@@ -0,0 +1,70 @@
+# Releasing
+
+A release is a git tag. Everything else follows from it, so the order matters.
+
+## 1. Prepare
+
+- [ ] `tests/run_tests.sh` is green, and CI is green on `main`.
+- [ ] `CHANGELOG.md` has a section for the new version with today's date, and the
+ `## [Unreleased]` section is empty.
+- [ ] The version is bumped in **three** places, which must agree:
+ `nextflow.config` (`manifest.version`), `CITATION.cff` (`version`), and the
+ version badge in `README.md`.
+
+## 2. Publish the container first
+
+The pipeline pins its image **by digest**, so the image has to exist before the
+tag that references it.
+
+```bash
+gh workflow run docker-publish.yml -f tag=
+```
+
+It builds once and pushes to both GHCR and Docker Hub. Copy the `sha256:` digest
+the action prints into `nextflow.config`:
+
+```groovy
+container = "docker://paururo/bampiro@sha256:"
+```
+
+Commit that, then verify a run resolves it:
+
+```bash
+nextflow run main.nf -profile test -stub-run
+```
+
+!!! note
+ Whichever registry the cluster pulls from has to be **public**, or Singularity
+ needs `APPTAINER_DOCKER_USERNAME` / `APPTAINER_DOCKER_PASSWORD` in the
+ environment. Docker Hub is the public one today.
+
+## 3. Tag and release
+
+```bash
+git tag
+git push --tags
+gh release create --title "" --notes-file
+```
+
+The tag triggers `docker-publish.yml` again, which is harmless: it rebuilds the
+same content and republishes the same tags.
+
+## 4. Archive for citation
+
+Currently **blocked**: Zenodo can only archive a public repository, and this one
+is private. When it goes public:
+
+1. Sign in to [Zenodo](https://zenodo.org) with GitHub and enable the switch for
+ `PathoGenOmics-Lab/BAMpiro`.
+2. Cut a release. Zenodo archives it and mints a DOI, using the metadata in
+ `.zenodo.json`.
+3. Add the **concept DOI** (the one that always resolves to the newest version,
+ not the per-version DOI) to `CITATION.cff` as `doi:` and to the README badge
+ row.
+
+## 5. Docs
+
+`docs.yml` publishes GitHub Pages on push to `main`, gated on the repository
+variable `ENABLE_PAGES`. Pages on a private repository needs a paid plan, so
+while this repository is private the site is off and the README's docs badges do
+not resolve. Nothing to do at release time beyond checking the build passed.
diff --git a/SECURITY.md b/SECURITY.md
new file mode 100644
index 0000000..ea16566
--- /dev/null
+++ b/SECURITY.md
@@ -0,0 +1,71 @@
+# Security Policy
+
+BAMpiro is a research pipeline, not a network service. It runs on a laptop or an HPC cluster, under
+your own account, on data you point it at. It listens on no port, has no users, no authentication and
+no persistent state beyond the files it writes to `--outdir` and the Nextflow work directory. That
+rules out most of what "vulnerability" usually means, so this policy is short and specific to what is
+actually plausible here.
+
+## Supported versions
+
+Only the latest release is supported. Fixes go into the next release; there are no backports to
+earlier tags.
+
+| Version | Supported |
+| :--- | :--- |
+| Latest release (see [CHANGELOG.md](CHANGELOG.md) and the [releases page](https://github.com/PathoGenOmics-Lab/BAMpiro/releases)) | Yes |
+| Anything older | No |
+
+## Reporting a vulnerability
+
+Please report privately rather than opening a public issue:
+
+- **GitHub private vulnerability reporting**: the *Security* tab of this repository, then
+ *Report a vulnerability*. This keeps the report visible only to the maintainers.
+- **Email**: paula.ruiz.rodriguez@csic.es
+
+Include what you found, how to reproduce it, and what an attacker would gain. If it involves an input
+file, a minimal example helps more than a real dataset (and please do not send us patient data).
+
+This is a small academic team with no on-call rotation, so we cannot promise a response time. We will
+acknowledge your report as soon as we see it and tell you what we intend to do about it. If you would
+like credit in the changelog entry for the fix, say so.
+
+## What is realistically in scope
+
+**The container's pinned dependencies.** Runs execute inside `paururo/bampiro`, which bundles
+BWA-MEM2, Samtools, BCFtools, FreeBayes, SnpEff, FastP, Kraken2, MUMmer, Python and a Java runtime,
+all pinned (versions are listed in [docs/installation.md](docs/installation.md)). A vulnerability in
+one of those reaches BAMpiro users through the image. Tell us if a bundled tool needs an urgent
+version bump.
+
+Two related notes. `nextflow.config` references the image by its mutable version **tag**, so the same
+tag can be rebuilt; re-pin to `docker://paururo/bampiro@sha256:` if you need byte-for-byte
+reproducibility, as documented in the installation page. And `modules/qc.nf` falls back to
+downloading `extract_kraken_reads.py` from `params.krakentools_url` over the network when the
+vendored copy under `bin/` is missing, which is a supply-chain surface worth knowing about; the
+vendored copy is there so the normal path never touches the network.
+
+**Untrusted input files.** The pipeline parses whatever you hand it: the TSV samplesheet, FASTQ
+reads, a reference FASTA, a GFF3, and optionally a Kraken2 database. Paths from the samplesheet reach
+shell commands inside process scripts. Running BAMpiro on files from a source you do not trust is
+therefore closer to running an untrusted script than to opening a data file. Treat a crafted
+samplesheet or reference that causes something other than a clean parse error as a bug worth
+reporting.
+
+**The QC report is your data.** Each run writes a single self-contained
+`_qc_report.html` that embeds sample identifiers, per-sample metrics and any metadata
+columns you supplied. It contains no remote references and phones nothing home, but it is not
+anonymised. Publishing or sharing one shares whatever was in your samplesheet.
+
+## What is out of scope
+
+- Denial of service through resource exhaustion. Kraken2 alone asks for around 80 GB of RAM by
+ design; an under-sized host being killed by the OOM reaper is a sizing problem, covered in
+ [Troubleshooting](docs/troubleshooting.md).
+- Anything requiring an attacker who already has your shell account, your cluster credentials or
+ write access to the work directory.
+- The security of Nextflow, Docker, Singularity or your scheduler. Report those upstream.
+
+Ordinary bugs, wrong results and crashes are not security issues. Open a
+[normal issue](https://github.com/PathoGenOmics-Lab/BAMpiro/issues/new/choose) for those.
diff --git a/assets/NO_FILE b/assets/NO_FILE
new file mode 100644
index 0000000..e69de29
diff --git a/assets/NO_FILE.README.md b/assets/NO_FILE.README.md
new file mode 100644
index 0000000..72bfc82
--- /dev/null
+++ b/assets/NO_FILE.README.md
@@ -0,0 +1,26 @@
+# `NO_FILE*` placeholders
+
+Nextflow has no notion of an optional process input: a `path` input must always
+receive something. When one of BAMpiro's optional features is off, `main.nf`
+stages one of these empty files instead, and the process recognises it by the
+`NO_FILE*` prefix and skips the corresponding argument.
+
+| File | Stands in for | Off when |
+| :--- | :--- | :--- |
+| `NO_FILE` | pathotypr lineage summary (`GENERATE_LEGACY_STATS`) | pathotypr did not type that sample |
+| `NO_FILE_DR` | drug-resistance calls (`QC_REPORT`) | `--run_pathotypr false` |
+| `NO_FILE_H37RV` | canonically annotated VCFs (`QC_REPORT`) | `--annotate_canonical false` |
+| `NO_FILE_LIFTOVER` | variant coordinate map (`QC_REPORT`) | `--variant_liftover false` |
+
+Two properties matter, and both are load-bearing:
+
+- **The names are distinct.** `QC_REPORT` takes three of these at once and all
+ three features are off by default. Nextflow refuses to stage two inputs of the
+ same task under the same filename, so sharing one `NO_FILE` fails the run with
+ `input file name collision`.
+- **The files exist and are empty.** A path that does not exist stages as a
+ dangling symlink, which happens to work with the default staging mode but fails
+ outright under `stageInMode 'copy'`. Empty keeps the `[ -s ... ]` guards false.
+
+They are referenced as `${projectDir}/assets/...`, so a file called `NO_FILE`
+sitting in someone's launch directory is never picked up by mistake.
diff --git a/assets/NO_FILE_DR b/assets/NO_FILE_DR
new file mode 100644
index 0000000..e69de29
diff --git a/assets/NO_FILE_H37RV b/assets/NO_FILE_H37RV
new file mode 100644
index 0000000..e69de29
diff --git a/assets/NO_FILE_LIFTOVER b/assets/NO_FILE_LIFTOVER
new file mode 100644
index 0000000..e69de29
diff --git a/bin/WGS_fasta_allpos.py b/bin/WGS_fasta_allpos.py
index 5859a55..62a418e 100644
--- a/bin/WGS_fasta_allpos.py
+++ b/bin/WGS_fasta_allpos.py
@@ -136,8 +136,11 @@ def _int(key: str) -> Optional[int]:
if ad and ad not in {".", ""}:
parts = ad.split(",")
try:
- ro = int(parts[0])
- ao = sum(int(x) for x in parts[1:]) if len(parts) > 1 else 0
+ # Both or neither: assigning ro before the ao sum can raise would leave a
+ # half-initialized (5, None) on a malformed AD such as "5,x".
+ ad_ro = int(parts[0])
+ ad_ao = sum(int(x) for x in parts[1:]) if len(parts) > 1 else 0
+ ro, ao = ad_ro, ad_ao
except ValueError:
pass
return ro, ao
@@ -187,11 +190,14 @@ def _load(self, path: str) -> None:
line = line.strip()
if not line or line.startswith("#"):
continue
- if line.lower().startswith("chrom"):
- continue
parts = line.split("\t")
if len(parts) < 3:
continue
+ # A header row is one whose coordinate columns are not numbers. Detecting it by
+ # the first field instead would drop a real interval on any contig named
+ # chrom*, e.g. "chromosome_1", which would then never be masked.
+ if not (parts[1].strip().isdigit() and parts[2].strip().isdigit()):
+ continue
chrom = parts[0].strip()
try:
start = int(parts[1]); end = int(parts[2])
@@ -319,7 +325,11 @@ def _load(self, fasta_path: str) -> None:
if line.startswith(">"):
if name is not None:
self.seqs[name] = "".join(parts).upper()
- name = line[1:].split()[0]
+ fields = line[1:].split()
+ if not fields:
+ raise RuntimeError(
+ f"{fasta_path}: FASTA header with no sequence name (a bare '>')")
+ name = fields[0]
parts = []
else:
parts.append(line.strip())
diff --git a/bin/backbone_allpos.awk b/bin/backbone_allpos.awk
new file mode 100644
index 0000000..bdb9843
--- /dev/null
+++ b/bin/backbone_allpos.awk
@@ -0,0 +1,54 @@
+# Turn a samtools mpileup stream into the all-positions backbone VCF body.
+#
+# samtools mpileup -aa ... | awk -v MINCOV=30 -v GAPDP=7 -v BAQ=1 -v OFS='\t' \
+# -f backbone_allpos.awk
+#
+# MINCOV params.allpos_min_cov - depth at or above this is a confident reference call
+# GAPDP params.consensus_min_dp - only read when BAQ=1
+# BAQ 1 when the input carries a second, BAQ-adjusted depth in column 10 (the two pileups
+# pasted together), 0 when it is a plain single mpileup
+#
+# The backbone supplies one record per reference position, which is what makes an all-sites VCF and
+# therefore a phylogenetic supermatrix possible. Both branches of CALL_BACKBONE used to carry their
+# own copy of count_bases; tests/unit/test_backbone_allpos.py covers this one.
+
+# Count reference and alternate observations in an mpileup base string.
+#
+# The format is not one character per read: a read start is '^' followed by a mapping-quality
+# character that must not be counted, and an indel is a sign, a length, and that many bases which
+# belong to the PREVIOUS position. Getting this wrong silently inflates or deflates the depth
+# evidence the consensus relies on.
+function count_bases(s, i, n, c, num) {
+ RO_G = 0; AO_G = 0; i = 1; n = length(s)
+ while (i <= n) {
+ c = substr(s, i, 1)
+ if (c == "^") { i += 2 } # ^ + mapping-quality char (read start)
+ else if (c == "$") { i++ } # read end
+ else if (c == "+" || c == "-") { # indel: skip sign + count + that many bases
+ i++; num = ""
+ while (i <= n && substr(s, i, 1) ~ /[0-9]/) { num = num substr(s, i, 1); i++ }
+ i += (num + 0)
+ }
+ else if (c == "." || c == ",") { RO_G++; i++ } # match to reference
+ else if (c ~ /[ACGTacgt]/) { AO_G++; i++ } # mismatch = alt observation
+ else { i++ } # * (del), <>, N, etc. -> ignore
+ }
+}
+
+{
+ chrom = $1; pos = $2; ref = $3; dp = ($4 + 0)
+ count_bases($5) # RO/AO from the no-BAQ pileup (real base evidence)
+
+ wt = (dp >= MINCOV ? 1 : 0)
+ nc = (dp >= MINCOV ? 0 : 1)
+ gt = (dp >= MINCOV ? "0" : "./.")
+
+ # A site is flagged baq_dropout when it is well covered but BAQ collapses it to at or below
+ # consensus_min_dp (the indel-adjacent homopolymer artefact), so the consensus masks it with X
+ # rather than leaving a gap. Without the BAQ column there is nothing to compare and no flag.
+ flt = "."
+ if (BAQ == 1) { dp_baq = ($10 + 0); if (dp > GAPDP && dp_baq <= GAPDP) flt = "baq_dropout" }
+
+ info = "ADP=" dp ";WT=" wt ";HET=0;HOM=0;NC=" nc
+ print chrom, pos, ".", ref, ".", ".", flt, info, "GT:DP:AD", gt ":" dp ":" RO_G "," AO_G
+}
diff --git a/bin/build_min_unique_len.py b/bin/build_min_unique_len.py
index aa01fda..2ba38a8 100755
--- a/bin/build_min_unique_len.py
+++ b/bin/build_min_unique_len.py
@@ -13,6 +13,8 @@
reference_start, so the read filter looks up min_unique_len[read.reference_start] directly.
"""
import argparse
+import sys
+
import numpy as np
@@ -41,7 +43,8 @@ def main():
for K, path in items: # ascending K -> first (smallest) unique wins
kv = np.uint16(min(K, a.sentinel - 1)) # never store the sentinel itself as a real length
for line in open(path):
- if not line or line[0] == "#":
+ # A blank line reads as "\n", never "", so testing `not line` alone never fires.
+ if not line.strip() or line[0] == "#":
continue
c, s, e, v = line.rstrip("\n").split("\t")
arr = tracks.get(c)
@@ -60,6 +63,22 @@ def main():
tail = (tracks[n] == SENT) & (~covered[n])
tracks[n][tail] = 0 # span >= 0 always -> such reads are kept
+ # A never-scored position becomes 0, i.e. "uniquely placeable at any read length", so
+ # the read filter keeps everything there. At a contig END that is deliberate: genmap
+ # cannot score a k-mer that runs off the sequence. A hole in the MIDDLE is a different
+ # thing entirely, and it would silently disable masking over that stretch, so say so
+ # rather than let the permissive default pass for a measurement.
+ if tail.any():
+ idx = np.flatnonzero(tail)
+ head_run = np.flatnonzero(covered[n])
+ first, last = (head_run[0], head_run[-1]) if head_run.size else (L, -1)
+ interior = idx[(idx > first) & (idx < last)]
+ if interior.size:
+ sys.stderr.write(
+ "[build_min_unique_len] WARNING: %s has %d unscored position(s) away from "
+ "the contig ends (first at %d); masking is disabled there\n"
+ % (n, interior.size, int(interior[0]) + 1))
+
np.savez(a.out_prefix + ".min_unique_len.npz", **tracks)
# Optional repeat BED: a base is unmappable iff NO read of length <= mask_window covering it
diff --git a/bin/build_snp_matrix.py b/bin/build_snp_matrix.py
index 427695c..905d235 100644
--- a/bin/build_snp_matrix.py
+++ b/bin/build_snp_matrix.py
@@ -60,6 +60,10 @@ def _af_dp(fmt, val):
alleles = [a for a in d.get('GT', './.').replace('|', '/').split('/') if a not in ('.', '')]
if alleles:
af = sum(1 for a in alleles if a != '0') / len(alleles)
+ # A single-element AD cannot give an allele fraction, but it still states a depth: without this
+ # the matrix emits an empty DP cell for a site whose depth is right there in the record.
+ if dp is None and ad:
+ dp = sum(ad)
return af, dp
diff --git a/bin/filter_reads_mappability.py b/bin/filter_reads_mappability.py
index 3ce5616..34f51cd 100755
--- a/bin/filter_reads_mappability.py
+++ b/bin/filter_reads_mappability.py
@@ -58,6 +58,8 @@ def anchored(A, idx, sp):
w = sys.stdout.write
for line in sys.stdin:
+ if not line.strip(): # a blank line has no fields to index
+ continue
if line[0] == "@":
if line.startswith("@SQ"):
for tok in line.split("\t"):
diff --git a/bin/format_snps_for_backbone.awk b/bin/format_snps_for_backbone.awk
new file mode 100644
index 0000000..a7a3a33
--- /dev/null
+++ b/bin/format_snps_for_backbone.awk
@@ -0,0 +1,51 @@
+# Rewrite a FreeBayes SNP VCF into the reduced form the backbone merge and the consensus step
+# expect: INFO becomes ADP/WT/HET/HOM/NC and the sample column becomes GT:DP.
+#
+# awk -v HETMIN=0.10 -f format_snps_for_backbone.awk valid_snps.vcf
+#
+# HETMIN is params.het_min_frac. Both FreeBayes processes pipe through this, so the genotype
+# decision below lives in exactly one place; tests/unit/test_format_snps_for_backbone.py covers it.
+
+BEGIN { OFS = "\t" }
+
+/^##/ { print; next }
+/^#CHROM/ { print; next }
+
+{
+ dp = 0; ro = 0; ao = 0
+
+ # INFO is tokenised on ';' and each key anchored (^DP=/^RO=/^AO=) so DP is never grabbed from a
+ # look-alike field such as MQMDP= or ADP=. AO is per-alt; after `bcftools norm -m -` the record is
+ # biallelic, so the first value is the one that matters.
+ ni = split($8, I, ";")
+ for (k = 1; k <= ni; k++) {
+ if (I[k] ~ /^DP=/) { v = I[k]; sub(/^DP=/, "", v); dp = v + 0 }
+ else if (I[k] ~ /^RO=/) { v = I[k]; sub(/^RO=/, "", v); ro = v + 0 }
+ else if (I[k] ~ /^AO=/) { v = I[k]; sub(/^AO=/, "", v); split(v, ax, ","); ao = ax[1] + 0 }
+ }
+
+ # Depth falls back through FORMAT/DP, then RO+AO, then 1 so the consensus never divides by zero.
+ # `+ 0` matters: FORMAT/DP may be ".", and without the coercion the comparisons below are string
+ # comparisons that "." fails, so every fallback is skipped and ADP=. reaches the consensus, which
+ # reads it as zero depth and turns a called SNP into a gap.
+ if (dp == 0) { n = split($9, fmt, ":"); m = split($10, dat, ":"); for (i = 1; i <= n; i++) if (fmt[i] == "DP") dp = dat[i] + 0 }
+ if (dp == 0) dp = ro + ao
+ if (dp == 0) dp = 1
+
+ split($10, b, ":"); gtype = b[1]
+
+ # Re-validate het: `bcftools norm -m -` splits a multiallelic 1/2 into biallelic records whose
+ # reference index 0 is an ARTIFACT (RO~0). A het whose reference allele is essentially
+ # unsupported is not a real het -> homozygous-alt (matches the af-based hom/het rule and stops
+ # the reference base leaking into the consensus IUPAC).
+ if ((gtype == "0/1" || gtype == "1/0") && (ro + ao) > 0 && ro / (ro + ao) < HETMIN) gtype = "1/1"
+
+ het = 0; hom = 0
+ if (gtype == "0/1" || gtype == "1/0") het = 1
+ else if (gtype == "1/1") hom = 1
+
+ $8 = "ADP=" dp ";WT=0;HET=" het ";HOM=" hom ";NC=0"
+ $9 = "GT:DP"
+ $10 = gtype ":" dp
+ print
+}
diff --git a/bin/pathotypr_liftover.py b/bin/pathotypr_liftover.py
index 39bbe68..aac8087 100644
--- a/bin/pathotypr_liftover.py
+++ b/bin/pathotypr_liftover.py
@@ -63,7 +63,9 @@ def _positions(path):
continue
c = line.replace(",", "\t").split("\t")
try:
- if len(c) >= 3 and c[1].lstrip("-").isdigit() and c[2].lstrip("-").isdigit():
+ # isdigit() only, not lstrip("-"): a BED coordinate is never negative, and accepting
+ # one silently expands to zero and negative positions that no reference has.
+ if len(c) >= 3 and c[1].isdigit() and c[2].isdigit():
a, b = int(c[1]), int(c[2]) # BED 0-based half-open -> 1-based positions a+1..b
for p in range(a + 1, b + 1):
pos.add(p)
diff --git a/bin/qc_report.py b/bin/qc_report.py
index fcea3e8..51ac807 100644
--- a/bin/qc_report.py
+++ b/bin/qc_report.py
@@ -13,1425 +13,32 @@
- Click any row, dot or point to highlight a sample everywhere; click a flag chip to filter.
With --gate it exits non-zero if any sample FAILs (at the default thresholds), so a downstream target can block.
Pure standard library (no numpy / matplotlib / CDN). Study-specific signatures live in a separate downstream gate.
+
+The implementation lives in the qcreport package next to this file: qcreport.parsers
+reads the pipeline's artefacts, qcreport.metrics holds the metric registry and the
+verdict engine, qcreport.panels builds the optional analysis panels and qcreport.render
+assembles the HTML. This file is the CLI: arguments in, payload assembled, two files out.
"""
from __future__ import annotations
import argparse
-import base64
-import gzip
-import html
-import json
import os
-import random
-import re
import sys
-import zlib
from datetime import datetime, timezone
-DEF = dict(depth_min=10.0, breadth_min=90.0, missing_max=10.0, dup_max=40.0,
- iupac_max=2.0, mapping_min=80.0, titv_min=0.0, snp_z=3.0,
- het_max_frac=3.0, mixed_min_frac=2.0) # het_max_frac + mixed_min_frac are PERCENTS
-
-# Ancient (aDNA) sample threshold view: looser cov/breadth/missing/mapping/dup/iupac (a 5x mummy must not FAIL the
-# modern gate), plus damage_min_ct = the terminal 5' C>T deamination floor for the DAMAGE_LOW authentication screen.
-ANC_DEF = dict(depth_min=3.0, breadth_min=30.0, missing_max=70.0, dup_max=60.0,
- iupac_max=5.0, mapping_min=20.0, damage_min_ct=0.05)
-FAIL_FLAGS = {"LOW_DEPTH", "LOW_BREADTH", "HIGH_MISSING", "NO_DATA"}
-NBINS = 200 # bins along the reference for the genome callability-landscape heatmap
-REPO_URL = "https://github.com/PathoGenOmics-Lab/BAMpiro" # surfaced in the report header + footer
-
-# metric definitions + typical acceptable ranges (surfaced as in-report help)
-DEFS = {
- "dose": ["Per-sample dose from the samplesheet (a numeric annotation, not a QC metric). Feeds the scatter/correlation axes and the dose x treatment test.", ""],
- "raw_reads": ["Read pairs/reads before trimming.", "context-dependent"],
- "trimmed_reads": ["Reads after fastp quality/adapter trimming.", ""],
- "read_length": ["Mean read length after trimming.", ""],
- "duplication_pct": ["% reads marked as duplicates (PCR/optical).", "typically < 30-40%"],
- "mapped_pct": ["% reads mapping to the reference.", "> 90% for a clean library"],
- "properly_paired_pct": ["% read pairs mapping in the expected orientation/insert.", "> 90%"],
- "mean_depth": ["Mean sequencing depth over the genome (sensitive to a few high-copy repeats; see Median depth).", ">= 10-20x for confident calls"],
- "median_depth": ["Median sequencing depth over the genome (robust to high-copy repeats that inflate the mean).", ">= 10-20x"],
- "breadth_pct": ["% of the genome covered at >= 1x. Note: 1x is not enough to CALL a base; see Callable % for usable genome.", ">= 90% at 1x to enter a phylogeny"],
- "breadth5x_pct": ["% of the genome covered at >= 5x.", "high for a confident consensus"],
- "breadth10x_pct": ["% of the genome covered at >= 10x.", ">= 90% for confident SNP calls"],
- "coverage_cv": ["Coverage evenness = coefficient of variation of depth (sd / mean). Lower is more even; spikes flag repeats/contamination.", "< ~0.5-1 for an even library"],
- "avg_base_quality": ["Mean per-base quality of mapped reads.", ">= 30"],
- "unmapped_pct": ["% of reads that did not map to the reference (an off-target / contamination proxy).", "low for a clean, on-target library"],
- "singleton_pct": ["% of reads mapped without their mate (approx: reads mapped minus reads mapped-and-paired).", "low"],
- "mapq0_pct": ["% of reads with MAPQ 0 (ambiguous multi-mapping, e.g. repeats / IS elements).", "low"],
- "error_rate_pct": ["samtools mismatch (error) rate, as a percent; elevated with damage or contamination.", "< ~1%"],
- "insert_size": ["Mean fragment/insert size (bp) of properly paired reads. For aDNA a short distribution is expected.", "library-dependent"],
- "insert_size_sd": ["Standard deviation of insert size; a very wide spread signals chimeras/contamination.", ""],
- "q20_pct": ["% of trimmed bases with quality >= Q20 (fastp).", "> 90%"],
- "q30_pct": ["% of trimmed bases with quality >= Q30 (fastp).", "> 80-90%"],
- "gc_pct": ["GC content of the trimmed reads (fastp); a large deviation from the reference GC flags contamination.", "match the reference GC"],
- "est_genome_cov": ["Estimated fold coverage = trimmed bases / genome length (fastp-based, pre-mapping).", "tracks mean_depth"],
- "total_variants": ["Total variants (SNPs + indels) called vs the reference.", "context-dependent"],
- "snps": ["SNPs called vs the reference. Dominated by true phylogenetic distance, NOT quality; the SNP outlier flag is a cohort distributional signal, not a defect.", "expected to be higher for divergent lineages"],
- "snp_density": ["SNP calls per callable megabase = snps / (callable_pct/100 * genome_len/1e6). Normalises SNP burden for how much of the genome was callable; compare within a lineage. Suppressed for a toy genome (< 100 kb).", "compare within a lineage"],
- "het_variants": ["Heterozygous variant calls. In a clonal/haploid organism these are expected ~0; a nonzero rate is consistent with mixture, contamination, or (for aDNA) post-mortem deamination.", "~0"],
- "homo_indels": ["Homozygous indels called.", ""],
- "ti_tv": ["Transition/transversion ratio (a spectrum sanity check, not a rate). Typical genome-wide MTBC ~1.7-2.0; M. bovis (A4) runs higher (~2.6) from its transition-biased spectrum, and aDNA deamination inflates it. LOW Ti/Tv (< ~1.5) is the QC concern (mapping/base-call noise is transversion-heavy); high is biology or damage, not a defect.", "~1.5-2.1 for human lineages"],
- "missing_pct": ["% of the consensus that is '-' or N (no confident base).", "< 10%"],
- "iupac_pct": ["% of the consensus that is an IUPAC ambiguity (mixed call).", "< 2%"],
- "callable_pct": ["% of the consensus that is a confident A/C/G/T (the honest, depth-aware usable-genome measure).", "> 90%"],
- "longest_n_run": ["Longest single run of consecutive N/gap in the consensus (one big deletion vs many scattered holes).", "context-dependent"],
- "n_gaps": ["Number of separate N/gap runs in the consensus.", "context-dependent"],
- "het_frac": ["Fraction of VARIANT genotypes that are heterozygous (het / (het+hom)). An illustrative within-sample mixture signal, not the textbook heterozygosity rate; the denominator is variant-count, so read it with care. Suppressed below a minimum variant count, and NOT applied to aDNA (deamination inflates it).", "~0 for a clonal isolate"],
- "n_lineages": ["Number of lineages (collapsed to main lineage) above the mixture cut-off AND above a minimum marker count; > 1 is consistent with a mixed infection or cross-contamination.", "1"],
- "damage_ct": ["Terminal 5' C>T deamination frequency (mapDamage2). Elevated in authentic aDNA; consistent with post-mortem damage but not proof of authenticity or endogeneity on its own.", "aDNA typically >= 0.05-0.40"],
- "damage_ga": ["Terminal 3' G>A deamination frequency (mapDamage2).", "aDNA typically >= 0.05-0.40"],
- "fraglen": ["Mean fragment (read) length from mapDamage2. Short (~30-70 bp) is expected for aDNA.", "aDNA ~30-70 bp"],
- "ann_high": ["snpEff HIGH-impact variants (stop_gained, frameshift, start/stop_lost, splice). A cohort-relative excess is a contamination / misassembly / wrong-reference-annotation signal, not a per-sample defect on its own.", "cohort-relative robust-z; high flags a look-here sample"],
- "missense_silent": ["Missense / synonymous variant-count ratio: a within-sample pN/pS PROXY (a spectrum sanity check, NOT dN/dS and NOT a selection test; no site-degeneracy correction, entangled with divergence and with how snpEff bins effects). Read within a lineage. For a real pN/pS use the eskaks panel.", "compare within a lineage; robust-z high tail only"],
- "lof_pct": ["Fraction of protein-coding-effect variants that are loss-of-function (stop_gained + frameshift + start/stop_lost + splice). Elevated LoF can be real pseudogenisation OR a frameshift storm from an indel-calling / assembly / annotation artifact.", "cohort-relative robust-z"],
- "coding_pct": ["Percent of annotated variants falling in a transcript (HIGH+MODERATE+LOW; MODIFIER is the intergenic / intron / UTR bucket). Tracks how gene-dense the reference annotation is; not a quality gate.", "reference-dependent, descriptive"],
- "annotated_pct": ["Percent of this sample's real variants that received an ANN annotation from snpEff. A low value flags a GFF-to-database mismatch: the annotation may be silently garbage.", "cohort-relative low + absolute floor"],
- "snpeff_warn": ["Variants carrying a snpEff ERROR/WARNING/INFO subfield (e.g. WARNING_TRANSCRIPT_INCOMPLETE, ERROR_CHROMOSOME_NOT_FOUND). A value near the total variant count means the reference database was built wrong for this contig; treat every annotation with suspicion.", "low; near-total = broken DB"],
-}
-# Contamination/mixture cut-offs (het_max_frac, mixed_min_frac) and the aDNA damage floor are ILLUSTRATIVE sensitivity
-# defaults, not validated gates; they are WARN-only and live-adjustable. Tune to your depth and organism.
-
-# (key, label, kind, direction). kind: int|float|pct ; direction: hi_good|hi_bad|neu
-METRICS = [
- ("raw_reads", "Raw reads", "int", "neu"),
- ("trimmed_reads", "Trimmed", "int", "neu"),
- ("read_length", "Read len", "float", "neu"),
- ("q20_pct", "Q20 %", "pct", "hi_good"),
- ("q30_pct", "Q30 %", "pct", "hi_good"),
- ("gc_pct", "GC %", "pct", "neu"),
- ("duplication_pct", "Dup %", "pct", "hi_bad"),
- ("mapped_pct", "Mapped %", "pct", "hi_good"),
- ("properly_paired_pct", "Proper %", "pct", "hi_good"),
- ("unmapped_pct", "Unmapped %", "pct", "hi_bad"),
- ("singleton_pct", "Singleton %", "pct", "hi_bad"),
- ("mapq0_pct", "MAPQ0 %", "pct", "hi_bad"),
- ("error_rate_pct", "Error %", "pct", "hi_bad"),
- ("insert_size", "Insert bp", "float", "neu"),
- ("insert_size_sd", "Insert sd", "float", "neu"),
- ("mean_depth", "Depth", "float", "hi_good"),
- ("median_depth", "Median depth", "float", "hi_good"),
- ("breadth_pct", "Breadth %", "pct", "hi_good"),
- ("breadth5x_pct", "Breadth 5x %", "pct", "hi_good"),
- ("breadth10x_pct", "Breadth 10x %", "pct", "hi_good"),
- ("coverage_cv", "Cov CV", "float", "hi_bad"),
- ("avg_base_quality", "Base Q", "float", "hi_good"),
- ("est_genome_cov", "Est cov", "float", "hi_good"),
- ("total_variants", "Variants", "int", "neu"),
- ("snps", "SNPs", "int", "neu"),
- ("het_variants", "Het", "int", "hi_bad"),
- ("homo_indels", "Indels", "int", "neu"),
- ("ti_tv", "Ti/Tv", "float", "neu"),
- ("missing_pct", "Missing %", "pct", "hi_bad"),
- ("iupac_pct", "IUPAC %", "pct", "hi_bad"),
- ("callable_pct", "Callable %", "pct", "hi_good"),
- ("longest_n_run", "Longest N-run", "int", "hi_bad"),
- ("n_gaps", "N-gaps", "int", "neu"),
- ("ann_high", "HIGH impact", "int", "hi_bad"),
- ("ann_moderate", "MODERATE", "int", "neu"),
- ("ann_low", "LOW impact", "int", "neu"),
- ("ann_modifier", "MODIFIER", "int", "neu"),
- ("missense_silent", "Missense/silent", "float", "neu"),
- ("lof_pct", "LoF %", "pct", "hi_bad"),
- ("coding_pct", "Coding %", "pct", "neu"),
- ("annotated_pct", "Annotated %", "pct", "hi_good"),
- ("snpeff_warn", "snpEff warn", "int", "hi_bad"),
-]
-DIST = ["mean_depth", "median_depth", "breadth_pct", "breadth10x_pct", "coverage_cv", "mapped_pct",
- "unmapped_pct", "properly_paired_pct", "duplication_pct", "mapq0_pct", "error_rate_pct",
- "missing_pct", "callable_pct", "iupac_pct", "snps", "ti_tv",
- "ann_high", "missense_silent", "annotated_pct"]
-
-# columns that are identifiers / metadata / already-surfaced -> NOT auto-detected as "extra" metrics
-KNOWN_KEYS = {k for k, *_ in METRICS} | {
- "sample_id", "lineage", "lineage_counts", "drug_resistance", "dr_counts", "date", "ancient",
- "consensus_id", "length", "genome_len", "mapped_reads", "homo_variants",
- "snp_profile", "het_profile", "indel_profile", # per-position variant tracks (huge strings, not metrics)
- "eff_missense", "eff_synonymous", "eff_stop_gained", "eff_stop_lost", "eff_start_lost",
- "eff_frameshift", "eff_inframe_indel", "eff_splice", "eff_intergenic", "eff_regulatory",
- "ann_present", "ann_annotated_frac", "ann_total", "ann_annotated", "coding_count", "lof_count", "ann_db_error",
-}
-_PCT_HINTS = ("_pct", "_rate", "percent", "fraction", "frac")
-_INT_HINTS = ("reads", "count", "_n", "num_", "variants", "snps", "indels", "bases")
-
-
-def parse_summary(path):
- rows = {}
- with open(path) as f:
- header = f.readline().rstrip("\n").split("\t")
- for line in f:
- if not line.strip():
- continue
- parts = line.rstrip("\n").split("\t")
- if len(parts) < len(header):
- parts += [""] * (len(header) - len(parts))
- d = dict(zip(header, parts))
- sid = (d.get("sample_id") or parts[0]).strip()
- if not sid:
- continue
- rows[sid] = d
- return rows
-
-
-def consensus_stats(path):
- op = gzip.open if str(path).endswith(".gz") else open
- name, chunks = None, []
- with op(path, "rt", encoding="utf-8", errors="replace") as fh:
- for line in fh:
- if line.startswith(">"):
- if name is None: # keep the first record's id, but keep reading ALL records so a
- name = line[1:].split()[0] # multi-contig consensus is measured whole-genome, not contig 1 only
- else:
- chunks.append(line.strip())
- seq = "".join(chunks).upper()
- L = len(seq) or 1
- gap, n = seq.count("-"), seq.count("N")
- iup = sum(seq.count(c) for c in "RYSWKMBDHV")
- acgt = sum(seq.count(c) for c in "ACGT")
- # gap-length view: one big deletion vs many scattered holes are very different for a phylogeny
- runs = re.findall(r"[-N]+", seq)
- longest_run = max((len(r) for r in runs), default=0)
- # binned missingness profile ALONG the reference (consensus is reference-aligned) -> the genome landscape heatmap.
- # A callability proxy for coverage: per bin, the fraction of positions that are N / gap (no confident base).
- prof = []
- for b in range(NBINS):
- s0, e0 = (b * L) // NBINS, ((b + 1) * L) // NBINS
- sub = seq[s0:e0]
- ln = len(sub) or 1
- prof.append(round((sub.count("-") + sub.count("N")) / ln, 3))
- return {"consensus_id": name, "length": len(seq), "missing_pct": 100.0 * (gap + n) / L,
- "iupac_pct": 100.0 * iup / L, "callable_pct": 100.0 * acgt / L,
- "longest_n_run": longest_run, "n_gaps": len(runs), "miss_profile": prof}
-
-
-def to_float(x):
- try:
- if x in (None, "", "NA", "."):
- return None
- return float(x)
- except (ValueError, TypeError):
- return None
-
-
-def clean_str(x):
- """A metadata string, or None for the NA/empty sentinels (so 'NA' lineage does not read as a real lineage)."""
- if x is None:
- return None
- s = str(x).strip()
- return s if s and s not in ("NA", ".", "-") else None
-
-
-_PROFILE_WARNED = set()
-
-
-def parse_profile(s):
- """A comma-joined per-bin variant-count profile string -> [int]*NBINS, or None if absent/unparseable.
- Coerced (pad with 0 / truncate) to NBINS: the genome-landscape JS assumes every profile has exactly
- NBINS bins so that a brushed span maps 1:1 onto the hotspot table. Upstream (stats_to_legacy.BIN_N)
- already emits NBINS, so a mismatch means a foreign / re-binned summary; warn once per observed length."""
- if not s or str(s) in ("NA", ".", ""):
- return None
- try:
- v = [int(float(x)) for x in str(s).split(",")]
- except (ValueError, TypeError):
- return None
- if not v:
- return None
- if len(v) != NBINS:
- if len(v) not in _PROFILE_WARNED:
- sys.stderr.write(
- "[qc_report] warning: variant profile has %d bins, expected %d; "
- "padding/truncating so the genome landscape stays aligned.\n" % (len(v), NBINS))
- _PROFILE_WARNED.add(len(v))
- v = (v + [0] * NBINS)[:NBINS]
- return v
-
-
-def parse_lineage_colors(path):
- """Read a 'lineage#hexdesc' palette (e.g. mycolorsTB) -> {lineage: '#hex'}. Empty if absent."""
- if not path or not os.path.exists(path):
- return {}
- out = {}
- try:
- with open(path, encoding="utf-8") as fh:
- for line in fh:
- if line.startswith("#") or not line.strip():
- continue
- c = line.rstrip("\n").split("\t")
- if len(c) >= 2 and c[1].strip().startswith("#") and len(c[1].strip()) in (4, 7):
- out[c[0].strip()] = c[1].strip()
- except OSError:
- return {}
- return out
-
-
-def parse_bed(path):
- """Parse a BED (0-based half-open) of regions to mask -> [(start, end)]. Empty if absent."""
- if not path or not os.path.exists(path):
- return []
- iv = []
- try:
- op = gzip.open if str(path).endswith(".gz") else open
- with op(path, "rt", encoding="utf-8", errors="replace") as fh:
- for line in fh:
- if not line.strip() or line.startswith(("#", "track", "browser")):
- continue
- c = line.split("\t")
- if len(c) >= 3:
- try:
- iv.append((int(c[1]), int(c[2])))
- except ValueError:
- pass
- except OSError:
- return []
- return iv
-
-
-def mask_profile(intervals, genome_len, nbins=NBINS):
- """Per-bin masked fraction (0-1) along the reference + the total masked %."""
- if not intervals or genome_len <= 0:
- return None, 0.0
- binbp = genome_len / nbins
- cov = [0.0] * nbins
- total = 0
- for s, e in intervals:
- s = max(0, s); e = min(genome_len, e)
- if e <= s:
- continue
- total += e - s
- b0 = int(s // binbp); b1 = min(nbins - 1, int((e - 1) // binbp))
- for b in range(b0, b1 + 1):
- bs = b * binbp; be = (b + 1) * binbp
- cov[b] += max(0.0, min(e, be) - max(s, bs))
- prof = [round(min(1.0, cov[b] / binbp), 3) if binbp > 0 else 0.0 for b in range(nbins)]
- return prof, round(100.0 * total / genome_len, 2)
-
-
-def parse_gene_burden(path):
- """Cohort gene-burden aux TSV (from collect_summary --gene-burden-out) -> list of gene dicts sorted by
- total_impactful desc. Optional file; empty on any problem. Bounded to the top 200 genes."""
- if not path or not os.path.exists(path):
- return []
- out = []
- try:
- with open(path, encoding="utf-8", errors="replace") as fh:
- header = fh.readline().rstrip("\n").split("\t")
- idx = {h: i for i, h in enumerate(header)}
-
- def sget(c, k):
- i = idx.get(k)
- return c[i] if (i is not None and i < len(c) and c[i] != "") else None
-
- def nget(c, k):
- v = sget(c, k)
- try:
- return int(float(v))
- except (TypeError, ValueError):
- return None
- for line in fh:
- if not line.strip():
- continue
- c = line.rstrip("\n").split("\t")
- gene = sget(c, "gene")
- if not gene:
- continue
- hi, mod = nget(c, "high") or 0, nget(c, "moderate") or 0
- tot = nget(c, "total_impactful")
- out.append({"gene": gene, "high": hi, "moderate": mod,
- "low": nget(c, "low"), "modifier": nget(c, "modifier"),
- "dominant_effect": sget(c, "dominant_effect") or sget(c, "dominant") or "NA",
- "n_samples": nget(c, "n_samples"),
- "total_impactful": tot if tot is not None else hi + mod,
- "start": nget(c, "start"), "end": nget(c, "end")})
- except Exception:
- return []
- out.sort(key=lambda r: (-(r.get("total_impactful") or 0), r["gene"]))
- return out[:200]
-
-
-def parse_pnps(path):
- """Optional cohort dN/dS (eskaks) per-gene TSV -> [{gene, pn, ps, pnps, n, effect}]. Empty on any problem
- (panel hidden). Real eskaks-aggregate columns: gene, gene_name, n_samples, n_pairs, mean_dN, mean_dS, dNdS,
- effect_dominant. inf/nan/empty are read as None, never blindly float()-ed."""
- if not path or not os.path.exists(path):
- return []
-
- def num(v):
- if v is None:
- return None
- v = str(v).strip().lower()
- if v in ("", "na", "nan", "inf", "-inf", "+inf"):
- return None
- return to_float(v)
- rows = []
- try:
- with open(path, encoding="utf-8", errors="replace") as fh:
- header = fh.readline().rstrip("\n").split("\t")
- idx = {h.strip().lower(): i for i, h in enumerate(header)}
-
- def col(c, p):
- i = idx.get(c)
- return p[i] if (i is not None and i < len(p)) else None
- for line in fh:
- if not line.strip():
- continue
- p = line.rstrip("\n").split("\t")
- g = col("gene_name", p) or col("gene", p)
- if not g:
- continue
- rows.append({"gene": g, "pn": num(col("mean_dn", p)), "ps": num(col("mean_ds", p)),
- "pnps": num(col("dnds", p)), "n": num(col("n_pairs", p)),
- "effect": col("effect_dominant", p)})
- except OSError:
- return []
- return rows[:300]
-
-
-_DR_DRUG_ORDER = ["RIF", "INH", "EMB", "PZA", "STR", "STM", "FQ", "LFX", "MFX", "OFX", "KAN", "AMK",
- "CAP", "ETH", "PTO", "LZD", "BDQ", "CFZ", "BDQ_CFZ", "DLM", "PAS", "CS"]
-
-
-def parse_dr(path):
- """Per-sample drug-resistance calls -> {'samples','drugs','calls':[{s,drug,gene,mutation,grade,gn,marker,
- af,dp}]} or None. Tab-separated with a header; recognised columns (case-insensitive): sample, drug, gene,
- mutation, grade, marker_name/marker, af, dp. `grade` may be a WHO string like '1) Assoc w R' -> gn=1."""
- if not path or not os.path.exists(path):
- return None
- calls, samples = [], []
- try:
- with open(path, encoding="utf-8", errors="replace") as fh:
- header = None
- for line in fh:
- if not line.strip() or line.startswith("#"):
- continue
- parts = line.rstrip("\n").split("\t")
- if header is None:
- header = [h.strip().lower() for h in parts]
- continue
- d = dict(zip(header, parts))
- s = (d.get("sample") or (parts[0] if parts else "")).strip()
- if not s:
- continue
- grade = (d.get("grade") or "").strip()
- mg = re.match(r"\s*(\d)", grade)
- dp = d.get("dp")
- try:
- dp = int(dp) if dp not in (None, "", "NA", ".") else None
- except ValueError:
- dp = None
- calls.append({"s": s, "drug": (d.get("drug") or "").strip(), "gene": (d.get("gene") or "").strip(),
- "mutation": (d.get("mutation") or "").strip(), "grade": grade,
- "gn": (int(mg.group(1)) if mg else None),
- "marker": (d.get("marker_name") or d.get("marker") or "").strip(),
- "af": to_float(d.get("af")), "dp": dp})
- if s not in samples:
- samples.append(s)
- except OSError:
- return None
- if not calls:
- return None
- drugs = sorted({c["drug"] for c in calls if c["drug"]},
- key=lambda x: (_DR_DRUG_ORDER.index(x) if x in _DR_DRUG_ORDER else 99, x))
- return {"samples": samples, "drugs": drugs, "calls": calls}
-
-
-def parse_kraken(paths):
- """Per-sample Kraken2 `.report` files (6 columns: clade%, clade reads, taxon reads, rank, taxid,
- name) -> {'samples':[{s, unclassified, classified, primary, secondary, top}]} or None. The sample id
- is the filename up to the first '__' (the pipeline writes '__.kraken.report'). The
- 'primary' taxon is the top species-level clade; contamination shows as a low primary% / large
- secondary taxon / high unclassified%."""
- out = []
- for p in paths or []:
- if not p or not os.path.exists(p):
- continue
- sid = os.path.basename(p).split("__")[0].split(".")[0]
- unclass, species = 0.0, []
- try:
- with open(p, encoding="utf-8", errors="replace") as fh:
- for line in fh:
- c = line.rstrip("\n").split("\t")
- if len(c) < 6:
- continue
- try:
- pct = float(c[0])
- except ValueError:
- continue
- rank, name = c[3].strip(), c[5].strip()
- if rank == "U":
- unclass = pct
- elif rank == "S": # species-level clades = the interpretable composition
- species.append((pct, name))
- except OSError:
- continue
- if not species and unclass == 0.0:
- continue
- species.sort(key=lambda x: -x[0])
- top = [{"name": n, "pct": round(p, 2)} for p, n in species[:6]]
- out.append({"s": sid, "unclassified": round(unclass, 2),
- "classified": round(100.0 - unclass, 2),
- "primary": (top[0] if top else None),
- "secondary": (top[1] if len(top) > 1 else None), "top": top})
- if not out:
- return None
- return {"samples": out}
-
-
-def parse_gff(path):
- """Best-effort GFF3 parse of gene/CDS features -> [{name, start, end}] (1-based). Empty on any problem."""
- if not path or not os.path.exists(path):
- return []
- genes = {}
- try:
- op = gzip.open if str(path).endswith(".gz") else open
- with op(path, "rt", encoding="utf-8", errors="replace") as fh:
- for line in fh:
- if line.startswith("#") or "\t" not in line:
- continue
- c = line.rstrip("\n").split("\t")
- if len(c) < 9 or c[2] not in ("gene", "CDS"):
- continue
- try:
- start, end = int(c[3]), int(c[4])
- except ValueError:
- continue
- name = None
- for key in ("Name=", "gene=", "locus_tag=", "ID="):
- i = c[8].find(key)
- if i >= 0:
- name = c[8][i + len(key):].split(";")[0].strip()
- break
- name = name or f"{start}-{end}"
- # prefer a 'gene' feature over a 'CDS' of the same name
- if name not in genes or (c[2] == "gene" and genes[name].get("type") != "gene"):
- genes[name] = {"name": name, "start": start, "end": end, "type": c[2]}
- except (OSError, ValueError):
- return []
- return sorted(({"name": g["name"], "start": g["start"], "end": g["end"]} for g in genes.values()),
- key=lambda g: g["start"])
-
-
-# Curated H37Rv locus tags (Mycobrowser) for the common MTB resistance genes, used as a fallback when
-# the run's reference GFF carries no H37Rv-style locus_tag for them (e.g. a non-H37Rv reference). These
-# are the WHO-catalogue first/second-line resistance genes.
-GENE_RV = {
- "rpoB": "Rv0667", "rpoC": "Rv0668", "katG": "Rv1908c", "inhA": "Rv1484", "fabG1": "Rv1483",
- "ahpC": "Rv2428", "gyrA": "Rv0006", "gyrB": "Rv0005", "embA": "Rv3794", "embB": "Rv3795",
- "embC": "Rv3793", "pncA": "Rv2043c", "rpsL": "Rv0682", "eis": "Rv2416c", "ethA": "Rv3854c",
- "ethR": "Rv3855", "gid": "Rv3919c", "gidB": "Rv3919c", "tlyA": "Rv1694", "Rv0678": "Rv0678",
- "atpE": "Rv1305", "pepQ": "Rv2535c", "folC": "Rv2447c", "thyA": "Rv2764c", "alr": "Rv3423c",
- "ddn": "Rv3547", "fgd1": "Rv0407", "fbiA": "Rv3261", "fbiB": "Rv3262", "fbiC": "Rv1173",
- "rrs": "MTB000019", "rrl": "MTB000020", "mmpL5": "Rv0676c", "mmpS5": "Rv0677c", "whiB7": "Rv3197A",
- "clpC1": "Rv3596c", "panD": "Rv3601c", "ddlA": "Rv2981c",
-}
-_RV_LOCUS_RE = re.compile(r'^(Rv\d|MTB\d)', re.I) # H37Rv (Mycobrowser) locus-tag schemes
-
-
-def parse_gene_locus(path):
- """{gene_display_name: locus_tag} from a GFF3 (e.g. rpoB -> Rv0667). Empty on any problem."""
- out = {}
- if not path or not os.path.exists(path):
- return out
- try:
- op = gzip.open if str(path).endswith(".gz") else open
- with op(path, "rt", encoding="utf-8", errors="replace") as fh:
- for line in fh:
- if line.startswith("#") or "\t" not in line:
- continue
- c = line.rstrip("\n").split("\t")
- if len(c) < 9 or c[2] not in ("gene", "CDS"):
- continue
- attrs = c[8]
-
- def attr(key):
- i = attrs.find(key)
- return attrs[i + len(key):].split(";")[0].strip() if i >= 0 else ""
- locus = attr("locus_tag=")
- name = attr("gene=") or attr("Name=")
- if name and locus and name not in out:
- out[name] = locus
- except (OSError, ValueError):
- return out
- return out
-
-
-def build_gene_map(gff_path):
- """{gene: Mycobrowser (H37Rv) locus tag}. Curated resistance-gene map, overlaid with any H37Rv-style
- locus tags from the run's GFF (so an H37Rv reference contributes every gene; a non-H37Rv reference's
- strain-specific tags are ignored in favour of the curated H37Rv equivalents)."""
- m = dict(GENE_RV)
- for gene, locus in parse_gene_locus(gff_path).items():
- if _RV_LOCUS_RE.match(locus):
- m[gene] = locus
- return m
-
-
-def robust(vals):
- v = sorted(x for x in vals if x is not None)
- if not v:
- return None, None
- n = len(v)
- med = v[n // 2] if n % 2 else (v[n // 2 - 1] + v[n // 2]) / 2
- devs = sorted(abs(x - med) for x in v)
- mad = devs[n // 2] if n % 2 else (devs[n // 2 - 1] + devs[n // 2]) / 2
- sigma = 1.4826 * mad
- if sigma <= 0:
- mean = sum(v) / n
- sigma = (sum((x - mean) ** 2 for x in v) / n) ** 0.5 or 1.0
- return med, sigma
-
-
-MIN_HET_VARIANTS = 50 # below this variant count the het fraction is too noisy to interpret -> None
-
-
-def het_frac(m):
- """Fraction of VARIANT genotypes that are heterozygous (het/(het+hom)). Illustrative mixture signal; the
- denominator is variant-count (not callable sites), so it is entangled with divergence -> read with care.
- Returns None below a minimum variant count (too noisy to interpret)."""
- het = to_float(m.get("het_variants"))
- hom = to_float(m.get("homo_variants"))
- if het is None or hom is None or (het + hom) < MIN_HET_VARIANTS:
- return None
- return het / (het + hom)
-
-
-def lineage_counts_parsed(lc):
- """Parse 'L4.1:3200;L4.3:120;L2:40' into {main_lineage: total_count}, COLLAPSING sub-lineages (L4.1 -> L4)
- so two sub-lineages of the same main lineage do not spuriously read as a mixed infection."""
- if not lc or str(lc) in ("NA", ".", ""):
- return {}
- out = {}
- for part in str(lc).replace(",", ";").split(";"):
- if ":" in part:
- k, v = part.rsplit(":", 1)
- n = to_float(v)
- if n is not None:
- main = k.strip().split(".")[0] # collapse sub-lineage to main lineage
- out[main] = out.get(main, 0.0) + n
- return out
-
-
-def lineage_fracs(lc):
- """{main_lineage: fraction} from the collapsed counts."""
- counts = lineage_counts_parsed(lc)
- tot = sum(counts.values())
- return {k: v / tot for k, v in counts.items()} if tot > 0 else {}
-
-
-def is_ancient(m):
- """True if the summary row marks this sample ancient (aDNA). Missing column -> False (modern)."""
- return str(m.get("ancient", "")).strip().lower() in ("yes", "y", "true", "1", "ancient")
-
-
-def infer_metric(key, values):
- """Infer (label, kind, dir) for an un-registered numeric summary column. Neutral dir = no auto-flag / no claim."""
- kl = key.lower()
- if any(h in kl for h in _PCT_HINTS):
- kind = "pct"
- elif all(v is None or float(v).is_integer() for v in values) and any(h in kl for h in _INT_HINTS):
- kind = "int"
- else:
- kind = "float"
- label = key.replace("_pct", " %").replace("_", " ").strip()
- label = label[:1].upper() + label[1:]
- return label, kind, "neu"
-
-
-def discover_extra_metrics(summ):
- """Any summary column not a KNOWN_KEY and numeric for >=1 sample -> an auto 'extra' metric (future-proofing)."""
- if not summ:
- return []
- seen = set()
- for mm in summ.values():
- for k in mm.keys():
- if k not in KNOWN_KEYS and k != "sample_id":
- seen.add(k)
- extras = []
- for k in sorted(seen):
- vals = [to_float(mm.get(k)) for mm in summ.values()]
- if any(v is not None for v in vals):
- label, kind, direction = infer_metric(k, vals)
- extras.append({"key": k, "label": label, "kind": kind, "dir": direction})
- return extras
-
-
-def _read_freq_file(path, max_pos=25):
- """Parse a mapDamage 5pCtoT_freq.txt / 3pGtoA_freq.txt ('posfreq'). Returns
- {'pos1': float|None, 'profile': [[pos,freq],...] up to max_pos} or None. Robust to header + blanks."""
- prof = []
- try:
- with open(path) as fh:
- for line in fh:
- line = line.strip()
- if not line or line.lower().startswith("pos"):
- continue
- parts = line.replace(",", "\t").split("\t")
- if len(parts) < 2:
- parts = line.split()
- if len(parts) < 2:
- continue
- p = to_float(parts[0]); v = to_float(parts[1])
- if p is None or v is None:
- continue
- prof.append((int(p), v))
- except (OSError, ValueError):
- return None
- if not prof:
- return None
- prof.sort(key=lambda x: x[0])
- pos1 = next((v for p, v in prof if p == 1), prof[0][1])
- return {"pos1": pos1, "profile": [[p, round(v, 5)] for p, v in prof[:max_pos]]}
-
-
-def _mean_fraglen(path):
- """Mean fragment length from lgdistribution.txt (StdLengthOccurences, '#'-comment lines)."""
- num = den = 0.0
- try:
- with open(path) as fh:
- for line in fh:
- if not line.strip() or line.startswith("#"):
- continue
- parts = line.split()
- if len(parts) < 3 or parts[1].lower() == "length":
- continue
- L = to_float(parts[1]); occ = to_float(parts[2])
- if L is None or occ is None:
- continue
- num += L * occ; den += occ
- except (OSError, ValueError):
- return None
- return (num / den) if den > 0 else None
-
-
-def mapdamage_stats(sample_dir):
- """Read a mapDamage2 output directory for one sample. Returns a dict or None if nothing parseable.
- keys: ct1 (5' C>T pos1), ga1 (3' G>A pos1), fraglen (mean), ct_profile [[pos,freq],...]."""
- if not sample_dir or not os.path.isdir(sample_dir):
- return None
- five = _read_freq_file(os.path.join(sample_dir, "5pCtoT_freq.txt"))
- three = _read_freq_file(os.path.join(sample_dir, "3pGtoA_freq.txt"))
- fl = _mean_fraglen(os.path.join(sample_dir, "lgdistribution.txt"))
- ct1 = five["pos1"] if five else None
- if five is None and three is None and fl is None and ct1 is None:
- return None
- return {"ct1": ct1, "ga1": (three["pos1"] if three else None), "fraglen": fl,
- "ct_profile": five["profile"] if five else None}
-
-
-def flag_sample(m, thr, snp_med, snp_sig, ancient=False, anc_thr=None, dmg=None):
- flags = []
- at = anc_thr if (ancient and anc_thr) else thr # ancient samples judged against the aDNA threshold view
- g = lambda k: to_float(m.get(k))
- # produced-nothing guard: a sample with no consensus AND no mapping metrics is a FAIL, not a quiet PASS-by-absence
- if all(g(k) is None for k in ("mean_depth", "breadth_pct", "missing_pct", "mapped_pct")):
- return "FAIL", ["NO_DATA"]
- if (v := g("mean_depth")) is not None and v < at["depth_min"]:
- flags.append("LOW_DEPTH")
- if (v := g("breadth_pct")) is not None and v < at["breadth_min"]:
- flags.append("LOW_BREADTH")
- if (v := g("missing_pct")) is not None and v > at["missing_max"]:
- flags.append("HIGH_MISSING")
- if (v := g("mapped_pct")) is not None and v < at["mapping_min"]:
- flags.append("MAPPING_LOW")
- if (v := g("duplication_pct")) is not None and v > at["dup_max"]:
- flags.append("HIGH_DUP")
- if (v := g("iupac_pct")) is not None and v > at["iupac_max"]:
- flags.append("HIGH_IUPAC")
- if thr["titv_min"] > 0 and (v := g("ti_tv")) is not None and v < thr["titv_min"]:
- flags.append("TITV_LOW")
- # SNP count is dominated by phylogenetic divergence, not quality -> a cohort distributional outlier (WARN only).
- # snp_med/sig are computed over MODERN samples; ancient SNP counts are low by design, so skip the z-flag for aDNA.
- if not ancient:
- snp = g("snps")
- if snp is not None and snp_med is not None and snp_sig:
- if snp < snp_med - thr["snp_z"] * snp_sig:
- flags.append("SNP_LOW")
- elif snp > snp_med + thr["snp_z"] * snp_sig:
- flags.append("SNP_HIGH")
- # het fraction: NOT applied to aDNA (post-mortem deamination inflates apparent het -> would false-flag mummies)
- if not ancient:
- hf = het_frac(m)
- if hf is not None and hf * 100.0 > thr["het_max_frac"]:
- flags.append("HET_HIGH")
- # MIXED: >1 MAIN lineage above the fraction cut-off AND above a minimum marker count (guards single-read noise)
- counts = lineage_counts_parsed(m.get("lineage_counts"))
- tot = sum(counts.values())
- if tot > 0 and sum(1 for c in counts.values()
- if (c / tot) * 100.0 >= thr["mixed_min_frac"] and c >= 3) > 1:
- flags.append("MIXED")
- # aDNA authentication: an ANCIENT sample whose terminal 5' C>T is below the floor -> possible modern contaminant
- if ancient and anc_thr is not None and dmg and dmg.get("ct1") is not None \
- and dmg["ct1"] < anc_thr["damage_min_ct"]:
- flags.append("DAMAGE_LOW")
- verdict = "FAIL" if any(f in FAIL_FLAGS for f in flags) else ("WARN" if flags else "PASS")
- return verdict, flags
-
-
-# ============================================================================= CSS / JS / SHELL
-# ---- report front-end assets (CSS / JS / HTML shell / logo) live next to this script in
-# report_assets/, one file per language so each stays editable with its own tooling; build_html()
-# splices them into the single self-contained output HTML. ----
-_ASSET_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "report_assets")
-
-
-def _asset(name):
- with open(os.path.join(_ASSET_DIR, name), encoding="utf-8") as fh:
- return fh.read()
-
-
-# CSS split into ordered modules under report_assets/css/ (base tokens -> layout -> components ->
-# panels); concatenated in cascade order, so the result is identical to a single stylesheet.
-_CSS_MODULES = ["css/01_base.css", "css/02_layout.css", "css/03_components.css", "css/04_panels.css"]
-CSS = "".join(_asset(m) for m in _CSS_MODULES)
-
-# The report front-end is one ES5 IIFE split into ordered modules under report_assets/js/ (one per
-# area) purely for editability; they are concatenated verbatim here, so the running code is identical
-# to a single file. Order matters: 01 opens the IIFE and sets up shared state; 14 wires events + closes it.
-_JS_MODULES = [
- "js/01_prelude.js", "js/02_qcspace.js", "js/03_state.js", "js/04_helpers.js",
- "js/05_insights_qc.js", "js/06_insights_genome.js", "js/07_render_core.js", "js/08_curation.js",
- "js/09_genome_genes.js", "js/10_dynamics.js", "js/11_epistasis.js", "js/12_snpmatrix.js",
- "js/13_drug_kraken.js", "js/14_boot.js",
-]
-JS = "".join(_asset(m) for m in _JS_MODULES)
-
-# BAMpiro header logo (bampiro2.png resized to ~90x100 and embedded so the report stays
-# self-contained - no external image request). Regenerate from .github/bampiro2.png if the logo changes.
-LOGO_DATA_URI = _asset("logo.b64")
-
-SHELL = _asset("shell.html") # named shell.html (not report_*.html) so .gitignore's report*.html rule can't swallow it
-
-
-def build_html(title, payload):
- # Embed the payload gzip-compressed + base64 so the whole cohort travels in a light, self-contained
- # HTML (JSON deflates ~5-10x); the browser inflates it natively at load. base64 has no "" so it
- # cannot break out of the