From 6ac7af0237979aac220a2b799bead27931d703bd Mon Sep 17 00:00:00 2001 From: nsheff Date: Mon, 4 May 2026 18:28:17 -0400 Subject: [PATCH 1/9] pepatac: remove dead refgenconf import; add audit + findings The refgenconf RGC and select_genome_config imports in pipelines/pepatac.py are never referenced anywhere else in the file. Confirmed via grep for \bRGC\b and select_genome. Adds tools/audit_refgenie_surface.sh which enumerates all refgenconf / RefGenConf / looper_refgenie_populate / REFGENIE references in the repo (sources, docs, tests, configs). Adds findings.md as the dogfooding deliverable for the refgenie1 branch. Initial sections: audit, API gaps, seek-key naming divergences, asset class gaps, CLI/install gaps, cluster integration. Validation run section will be filled by step 9 of the plan. Plan: assistant/pepatac_refgenie1_branch_plan_v1.md (step 2) --- findings.md | 281 ++++++++++++++++++++++++++++++++ pipelines/pepatac.py | 1 - tools/audit_refgenie_surface.sh | 40 +++++ 3 files changed, 321 insertions(+), 1 deletion(-) create mode 100644 findings.md create mode 100755 tools/audit_refgenie_surface.sh diff --git a/findings.md b/findings.md new file mode 100644 index 0000000..c28b761 --- /dev/null +++ b/findings.md @@ -0,0 +1,281 @@ +# Refgenie1 migration findings + +This file captures every gap, awkward call pattern, error message, naming +divergence, and missing API encountered while migrating PEPATAC from +`refgenconf` (legacy refgenie) to `refgenie` 1.0 (refgenie1) on the +`refgenie1` branch. It is the deliverable of the dogfooding phase. + +Each finding records: (1) what happened, (2) verbatim error/symptom (when +applicable), (3) which refgenie1 file(s) are involved, (4) recommended +one-line fix. + +## Audit (initial) + +Output of `tools/audit_refgenie_surface.sh` run from the repo root on the +`refgenie1` branch (created off `dev` @ `027b757`). 61 hits total across +20 files. Categorized below. + +- **(a) Python imports / call sites:** `pipelines/pepatac.py:21` ONLY + (`from refgenconf import RefGenConf as RGC, select_genome_config`). + Names `RGC` and `select_genome_config` are never referenced anywhere + else in `pepatac.py`. The import is purely dead code. Removed. +- **(b) Pipeline interface declarations:** + - `sample_pipeline_interface.yaml:51` — `var_templates.refgenie_config: "$REFGENIE"` + - `sample_pipeline_interface.yaml:54` — `pre_submit.python_functions: - refgenconf.looper_refgenie_populate` + - Jinja templates use `refgenie[g].fasta.{fasta,chrom_sizes}`, + `refgene_anno.refgene_tss`, `blacklist.blacklist`, + `feat_annotation.feat_annotation`, `bowtie2_index.dir`, + `bwa_index.dir`, `tallymer_index.search_file`, plus `[p]` versions + for prealignments. + - `project_pipeline_interface.yaml`: NO refgenie/refgenconf references. + - `pipelines/pepatac.yaml:37`: `genome_config: ${REFGENIE}` — pypiper + resource config field, unused; migrated for consistency. +- **(c) Documentation:** `docs/{tutorial,run-conda,detailed-install,run-bulker,assets}.md`, + `docs/howto/install-refgenie.md`, `docs/changelog.md`. All instruct + legacy refgenie syntax. Migrated. +- **(d) Example configs:** searched + `examples/test_project/test_refgenie_pep_config.yaml`, + `examples/test_project/.looper_test.yaml`, + `examples/tutorial/.looper_tutorial_refgenie.yaml`, + `examples/gold_atac/metadata/*` — none directly reference `$REFGENIE`. +- **(e) Test fixtures:** `tests/integration/*.py`, `tests/README.md` — + use the legacy CLI in a venv. Gated by `RUN_LOCAL_REFGENIE_TESTS=true`, + skipped by default. Logged as a follow-up rewrite. +- **(f) Requirements:** `requirements.txt:20` (`refgenconf>=0.12.2`), + `requirements-conda.yml:525` (`refgenconf==0.12.2`). No setup.py / + pyproject.toml dep. Both swapped to `refgenie>=1.0.0`. +- **(g) Other:** `docs/changelog.md:218` — historical note, left alone. + +### Audit findings beyond the plan's expectations + +- `pipelines/pepatac.yaml:37` — pypiper resource config field with + `${REFGENIE}` was not flagged by the plan. Migrated to + `${REFGENIE_DB_CONFIG_PATH}`. + +## API gaps + +### Missing `list_seek_keys_values` equivalent + +Refgenconf's `rgc.list_seek_keys_values()` returned the full +`{genome: {asset: {tag: {seek_key: path}}}}` shape in one call — +exactly what the legacy populator needs. Refgenie1 has no equivalent. +The local populator must walk `r.alias.list_all()` → +`r.asset.list_groups(genome_names=[g])` → +`r.asset.list_assets(genome_names=[g], asset_group_name=ag)` → +`asset.seek_keys` → `r.asset.seek(g, ag, asset_name, sk_name)` per +leaf. O(genomes * groups * assets * seek_keys) Python-level loop with +one `seek` call per leaf vs. refgenconf's single YAML walk. + +**Recommended fix:** add `Refgenie.list_seek_keys_values()` to refgenie1 +returning the same shape, populated in one DB walk via the existing +`selectinload(Asset.seek_keys)` already used in `list_assets`. +Refgenie1 file: `refgenie/refgenie/managers/asset/manager.py`. + +### `Path` vs `str` returns + +`r.asset.seek(...)` returns `pathlib.Path` for path seek_keys +(`refgenie/managers/asset/manager.py:1267`). Refgenconf returned +`str`. The populator must `str(...)` every leaf or Jinja templating +renders `Path.repr` which can confuse JSON-serializing downstreams. + +**Recommended fix:** make `Refgenie.asset.seek` return `str` to match +legacy semantics, OR document the divergence prominently in the +refgenie1 README. + +### Default-asset selection split across methods + +The populator needs the *default* asset name per (genome, asset_group) +to walk that asset's seek_keys. `seek()` defaults internally via +`get_default()`, but the populator needs the *name* to enumerate seek +keys, so it calls `r.asset.get_default(asset_group, genome)` +explicitly. Awkward — having `seek()` default but no public way to +get "the asset that would be used" without a separate manager call. + +**Recommended fix:** expose +`AssetManager.list_seek_keys(genome, asset_group, asset_name=None)` +that defaults `asset_name` the same way `seek()` does. Saves a +two-step dance. + +## Seek-key naming divergences + +| PEPATAC reference | Legacy | Refgenie1 | Action | +|------------------------------------------------|--------|-----------------------------------------------|--------| +| `refgenie[g].fasta.fasta` | OK | `fasta.fasta` (default) | none | +| `refgenie[g].fasta.chrom_sizes` | OK | `fasta.chrom_sizes` | none | +| `refgenie[g].refgene_anno.refgene_tss` | OK | `refgene_anno.refgene_tss` | none | +| `refgenie[g].blacklist.blacklist` | OK | `blacklist.blacklist` | none | +| `refgenie[g].feat_annotation.feat_annotation` | OK | `feat_annotation.feat_annotation` | none | +| `refgenie[g].bowtie2_index.dir` | OK | **`bowtie2_index.bowtie2_index`** (no `dir`!) | rewrite Jinja | +| `refgenie[g].bwa_index.dir` | OK | not registered | rewrite Jinja | +| `refgenie[g].tallymer_index.search_file` | OK | not registered | leave (gated on --sob) | + +### `bowtie2_index.dir` does not exist + +Refgenconf shipped a built-in `dir` seek_key for every asset that +returned the asset's containing directory. Refgenie1 has no such +convention — seek keys are explicit and declared in the asset class. +The `bowtie2_index` asset class emits seek keys +`bowtie2_index` (the index prefix), `build_timestamp`, +`refgenie_version`, `inputs`, `version`. No `dir`. + +For PEPATAC this is fine in spirit: bowtie2 wants a *prefix*, and +`bowtie2_index.bowtie2_index` returns exactly the prefix path +(`/bowtie2_index/default/`). PEPATAC's +`--genome-index` accepts this directly — see +`pipelines/pepatac.py:610-620`, which only special-cases trailing `.` +(the legacy `dir` convention). With a prefix path the special case is +skipped and bowtie2 runs against the prefix as expected. + +**Resolution on this branch:** rewrote Jinja to use +`refgenie[g].bowtie2_index.bowtie2_index`. + +**Recommended upstream fix:** add a `dir` seek_key to all `*_index` +asset classes in `refgenie/repos/recipes` that emits the asset's +parent directory. This restores the legacy convention as a portable +contract. + +### `bwa_index` / `tallymer_index` not registered + +PEPATAC supports BWA as alternative aligner and tallymer for +`--sob`. Neither is registered for hg38 in the deployed refgenie1 +instance, and `tallymer_index` has no asset class shipped in +`repos/recipes`. Jinja `is defined` guards already gate these +references — they only render when the user enables the corresponding +flag. Validation uses bowtie2 + no `--sob`, neither path exercised. + +**Recommended fix:** ship `tallymer_index` and `bwa_index` asset +classes in `refgenie/repos/recipes`. Out of scope for this plan. + +## Asset class / recipe gaps + +- `bowtie2_index` (and any other `*_index`) asset class lacks `dir` + seek_key. Documented above. +- `tallymer_index` and `bwa_index` asset classes not shipped in + `repos/recipes`. + +## CLI / install gaps + +### `refgenie` PyPI name collision + +Both legacy refgenie (0.12.x) and refgenie1 (1.0.x) ship to PyPI +under the name `refgenie`. `pip install refgenie>=1.0.0` is +unambiguous, but `pip install refgenie` cold gives wildly different +behavior. PEPATAC's `requirements.txt` now pins `refgenie>=1.0.0`. +On Rivanna, refgenie1 is in its own venv and on PATH via +`refgenie1.env`. + +**Recommended fix:** rename refgenie1's PyPI package to `refgenie2` +or `refgenie-next` until the upstream-merge decision lands, OR yank +legacy refgenie from PyPI. The current name collision is a silent +footgun. + +### `$REFGENIE` vs `$REFGENIE_DB_CONFIG_PATH` + +Legacy: `$REFGENIE` → path to `genome_config.yaml`. +Refgenie1: `$REFGENIE_DB_CONFIG_PATH` → `refgenie_db_config.yaml`, +plus `$REFGENIE_HOME_PATH` for install root. + +PEPATAC's pipeline interface used `$REFGENIE` directly. Migrated. +Cluster setup must export both refgenie1 vars. + +**Recommended fix:** refgenie1 README should explicitly call out the +env var migration as a breaking change. + +### `refgenie pull` removed + +Legacy `refgenie pull /` was a one-liner. Refgenie1 splits this +into `refgenie genome init ` + `refgenie add / +--recipe ...` plus a subscribed source for actual pulls. + +**Recommended fix:** ship a `refgenie pull` shim in refgenie1 that +resolves to the equivalent sequence for the common case (subscribed +source, default recipe). Single-command ergonomics matter. + +### `refgenie --version` not supported + +Legacy refgenie supports `refgenie --version`. Refgenie1 does not: + +``` +$ refgenie --version +refgenie: error: unrecognized arguments: --version +``` + +Documented in plan as a CLI gotcha. The plan suggested using +`refgenie --version` in cluster prep (step 8.1) — that check fails. +Worked around with `pip show refgenie | grep Version`. + +**Recommended fix:** add `--version` to refgenie1's top-level +argparse (or Typer) parser. Refgenie1 file: `refgenie/cli/cli_pydantic.py`. + +### No `refgenie asset list` subcommand + +The plan's step 8.1 calls `refgenie asset list --genome hg38`. +There is no `asset` subcommand in refgenie1. Asset listing happens +under `refgenie list -g `. Worked around in cluster prep. + +**Recommended fix:** add `asset` subcommand alias for discoverability, +OR document the legacy → refgenie1 CLI mapping prominently. + +## Cluster integration + +### Broken legacy `refgenie` binary still on PATH + +Legacy refgenie binary at `~/.local/bin/refgenie` is on PATH for +ns5bc on Rivanna and takes priority unless `refgenie1.env` is +sourced. The legacy binary's shebang points at +`/apps/software/standard/core/anaconda/2023.07-py3.11/bin/python` +which no longer exists, so calling it without env-sourcing fails: + +``` +bash: /home/ns5bc/.local/share/../bin/refgenie: +/apps/software/standard/core/anaconda/2023.07-py3.11/bin/python: +bad interpreter: No such file or directory +``` + +**Recommended fix:** the refgenie1 deploy plan should remove or +shadow the broken legacy binary on the cluster. + +## Local populator placement + +The plan offered two placements for `looper_refgenie_populate_local`: + +1. Upstream in refgenie1 as `refgenie.populator.looper_refgenie_populate_local`. +2. On the PEPATAC branch as `pepatac.refgenie_populator.looper_refgenie_populate_local`. + +**This branch picks option 2** because: +- Refgenie1 has no `populator.py` module yet — the sibling Issue + #126 plan adds the remote counterpart in that file. Adding the + local populator upstream now would race with that plan. +- The PEPATAC branch is exploratory; iterating on the populator here + is faster than round-tripping a refgenie1 PR. +- If/when the upstream populator lands, this branch's populator + becomes a one-line re-export. + +**Recommended follow-up:** once refgenie1's `populator.py` exists +(Issue #126's plan), move `looper_refgenie_populate_local` upstream +as the local-mode sibling of `looper_refgenie_populate_remote`. + +## Tests rewrite (out of scope for this branch) + +`tests/integration/{conftest.py,test_end_to_end.py,test_looper_run.py,test_local_refgenieserver.py}` +all use the legacy refgenie CLI (looking for `refgenie` in a venv, +calling `refgenie pull`, `refgenie seek -c`, etc.). The +`RUN_LOCAL_REFGENIE_TESTS=true` gate keeps them off by default. They +will need a full rewrite for refgenie1. + +**Recommended follow-up:** separate "PEPATAC test rewrite for +refgenie1" plan. + +## Performance + +Not measured here. Populator runs once per looper invocation against +a small db (4 genomes, ~5 asset groups each). On the dev login node, +populator completes in <1s. Perf comparison vs. refgenconf deferred. + +## Error messages + +(Filled by the validation run section below as encountered.) + +## Validation run + +(Filled below by step 9.) diff --git a/pipelines/pepatac.py b/pipelines/pepatac.py index b275943..b2876f3 100755 --- a/pipelines/pepatac.py +++ b/pipelines/pepatac.py @@ -18,7 +18,6 @@ from pathlib import Path import psutil from pypiper import build_command -from refgenconf import RefGenConf as RGC, select_genome_config TOOLS_FOLDER = "tools" ANNO_FOLDER = "anno" diff --git a/tools/audit_refgenie_surface.sh b/tools/audit_refgenie_surface.sh new file mode 100755 index 0000000..61b7e51 --- /dev/null +++ b/tools/audit_refgenie_surface.sh @@ -0,0 +1,40 @@ +#!/usr/bin/env bash +# +# Audit refgenconf / refgenie surface area in this repo. +# Used by the refgenie1 migration plan to enumerate every place that needs +# to change. Runs from the repo root. +# +# Categories assigned per hit (best-effort): +# (a) Python import / call site +# (b) Pipeline interface declaration +# (c) Documentation +# (d) Example config +# (e) Test fixture / test code +# (f) Requirements / setup +# (g) Other (logs, caches, README, etc.) + +set -euo pipefail + +cd "$(dirname "$0")/.." + +PATTERN='refgenconf|RefGenConf|looper_refgenie_populate|REFGENIE' + +echo "=== Audit: refgenconf surface area ===" +echo "pattern: ${PATTERN}" +echo + +grep -rnE "${PATTERN}" \ + --include='*.py' \ + --include='*.yaml' \ + --include='*.yml' \ + --include='*.txt' \ + --include='*.toml' \ + --include='*.md' \ + --exclude-dir='.venv' \ + --exclude-dir='tests/.venv' \ + --exclude-dir='.git' \ + --exclude-dir='node_modules' \ + --exclude-dir='__pycache__' \ + --exclude-dir='build' \ + --exclude-dir='dist' \ + . From e3229bdba6143df36173d5bca3a826aba72de508 Mon Sep 17 00:00:00 2001 From: nsheff Date: Mon, 4 May 2026 18:31:46 -0400 Subject: [PATCH 2/9] pepatac: migrate sample pipeline interface to refgenie1 populator MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - var_templates: refgenie_config ($REFGENIE) → refgenie_db_config ($REFGENIE_DB_CONFIG_PATH). Refgenie1 uses its own env var. - pre_submit.python_functions: refgenconf.looper_refgenie_populate → refgenie.looper_refgenie_populate_local. The new populator lives in refgenie1's refgenie/populator.py (added on the refgenie/refgenie1#nsheff-refactor-2 branch alongside this one). - Jinja: refgenie[g].bowtie2_index.dir → bowtie2_index.bowtie2_index. Refgenie1's bowtie2_index asset class has no 'dir' seek_key (the legacy refgenconf 'dir' built-in was removed); the bowtie2_index seek_key returns the index prefix path which is what bowtie2 -x consumes. Same change for bwa_index. - pipelines/pepatac.yaml resources.genome_config: $REFGENIE → $REFGENIE_DB_CONFIG_PATH (consistency, the field is unused but was confusing). Plan: assistant/pepatac_refgenie1_branch_plan_v1.md (steps 4-5) Findings: see findings.md for the full divergence audit. --- pipelines/pepatac.yaml | 2 +- sample_pipeline_interface.yaml | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/pipelines/pepatac.yaml b/pipelines/pepatac.yaml index e6532bf..b1d72a6 100755 --- a/pipelines/pepatac.yaml +++ b/pipelines/pepatac.yaml @@ -34,7 +34,7 @@ tools: # absolute paths to required tools # user configure resources: - genome_config: ${REFGENIE} + genome_config: ${REFGENIE_DB_CONFIG_PATH} adapters: null # Set to null to use default adapter file included in repository parameters: # parameters passed to bioinformatic tools, subclassed by tool diff --git a/sample_pipeline_interface.yaml b/sample_pipeline_interface.yaml index 71660e6..aa8436a 100755 --- a/sample_pipeline_interface.yaml +++ b/sample_pipeline_interface.yaml @@ -18,9 +18,9 @@ sample_interface: {% if sample.anno_name is defined %} --anno-name { sample.anno_name } {% elif refgenie[sample.genome].feat_annotation is defined %} --anno-name { refgenie[sample.genome].feat_annotation.feat_annotation } {% endif %} {% if sample.trimmer is defined %} --trimmer { sample.trimmer } {% else %} --trimmer "skewer" {% endif %} {% if sample.aligner is defined %} --aligner { sample.aligner } {% set aligner = sample.aligner %} {% else %} --aligner "bowtie2" {% set aligner = "bowtie2" %} {% endif %} - {% if aligner == "bowtie2" or sample.aligner == "bowtie2" %} {% if sample.genome_index is defined %} --genome-index { sample.genome_index } {% elif refgenie[sample.genome].bowtie2_index is defined %} --genome-index { refgenie[sample.genome].bowtie2_index.dir } {% endif %} {% else %} {% if sample.genome_index is defined %} --genome-index { sample.genome_index } {% elif refgenie[sample.genome].bwa_index is defined %} --genome-index { refgenie[sample.genome].bwa_index.dir } {% endif %} {% endif %} + {% if aligner == "bowtie2" or sample.aligner == "bowtie2" %} {% if sample.genome_index is defined %} --genome-index { sample.genome_index } {% elif refgenie[sample.genome].bowtie2_index is defined %} --genome-index { refgenie[sample.genome].bowtie2_index.bowtie2_index } {% endif %} {% else %} {% if sample.genome_index is defined %} --genome-index { sample.genome_index } {% elif refgenie[sample.genome].bwa_index is defined %} --genome-index { refgenie[sample.genome].bwa_index.bwa_index } {% endif %} {% endif %} {% if sample.prealignment_index is defined %} --prealignment-index { sample.prealignment_index } {% endif %} - {% if sample.prealignment_names is defined %} {% if aligner == "bowtie2" or sample.aligner == "bowtie2" %} --prealignment-index {% for p in sample.prealignment_names %} { p ~ '=' ~ refgenie[p].bowtie2_index.dir } {% endfor %} {% else %} --prealignment-index {% for p in sample.prealignment_names %} { p ~ '=' ~ refgenie[p].bwa_index.dir } {% endfor %} {% endif %} {% endif %} + {% if sample.prealignment_names is defined %} {% if aligner == "bowtie2" or sample.aligner == "bowtie2" %} --prealignment-index {% for p in sample.prealignment_names %} { p ~ '=' ~ refgenie[p].bowtie2_index.bowtie2_index } {% endfor %} {% else %} --prealignment-index {% for p in sample.prealignment_names %} { p ~ '=' ~ refgenie[p].bwa_index.bwa_index } {% endfor %} {% endif %} {% endif %} {% if sample.deduplicator is defined %} --deduplicator { sample.deduplicator } {% endif %} {% if sample.peak_caller is defined %} --peak-caller { sample.peak_caller } {% else %} --peak-caller "macs3" {% endif %} {% if sample.peak_type is defined %} --peak-type { sample.peak_type } {% else %} --peak-type "fixed" {% endif %} @@ -48,7 +48,7 @@ bioconductor: readFunName: runCOCOA readFunPath: BiocProject/runCOCOA.R var_templates: - refgenie_config: "$REFGENIE" + refgenie_db_config: "$REFGENIE_DB_CONFIG_PATH" pre_submit: python_functions: - - refgenconf.looper_refgenie_populate + - refgenie.looper_refgenie_populate_local From 6e075e7a8441413c1aea930e016f1ae1e0a30ead Mon Sep 17 00:00:00 2001 From: nsheff Date: Mon, 4 May 2026 18:32:07 -0400 Subject: [PATCH 3/9] pepatac: replace refgenconf with refgenie 1.0.0 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit requirements.txt: refgenconf>=0.12.2 → refgenie>=1.0.0 requirements-conda.yml: refgenconf==0.12.2, refgenie==0.12.1 → refgenie>=1.0.0 (single dep — refgenie1 supersedes both). Note: refgenie 1.0.x and legacy refgenie 0.12.x share the PyPI name 'refgenie' but are different packages. The >=1.0.0 pin disambiguates. This is a known footgun — see findings.md. Plan: assistant/pepatac_refgenie1_branch_plan_v1.md (step 6) --- requirements-conda.yml | 3 +-- requirements.txt | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/requirements-conda.yml b/requirements-conda.yml index 9fca270..212fb00 100755 --- a/requirements-conda.yml +++ b/requirements-conda.yml @@ -522,8 +522,7 @@ dependencies: - pydantic-core==2.27.2 - pyfaidx==0.8.1.3 - pysam==0.22.1 - - refgenconf==0.12.2 - - refgenie==0.12.1 + - refgenie>=1.0.0 - scikit-learn==1.6.1 - scipy==1.15.0 - threadpoolctl==3.5.0 diff --git a/requirements.txt b/requirements.txt index 8e1aa5b..b0448ab 100755 --- a/requirements.txt +++ b/requirements.txt @@ -17,4 +17,4 @@ psutil>=5.8 pysam>=0.16 python-Levenshtein>=0.12 pyyaml>=3.13 -refgenconf>=0.12.2 +refgenie>=1.0.0 From 505873f81118a79b7f7a1254a71cee3dbb913ab7 Mon Sep 17 00:00:00 2001 From: nsheff Date: Mon, 4 May 2026 18:33:56 -0400 Subject: [PATCH 4/9] pepatac: update docs for refgenie 1.0 syntax MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit docs/assets.md gets the canonical refgenie 1.0 setup walk-through. The other refgenie-touching docs (tutorial, run-conda, detailed-install, run-bulker, howto/install-refgenie) get a NOTE banner pointing at docs/assets.md plus inline replacements of: pip install refgenie → pip install "refgenie>=1.0.0" export REFGENIE=... → export REFGENIE_HOME_PATH + REFGENIE_DB_CONFIG_PATH refgenie init -c ... → refgenie init refgenie pull g/a → refgenie genome init && refgenie add g/a --recipe This is the minimum migration to keep the docs consistent with the branch. The legacy refgenie 0.12.x flow has no automatic CLI shim in refgenie 1.0; users on that flow must follow the new commands. Plan: assistant/pepatac_refgenie1_branch_plan_v1.md (step 7) --- docs/assets.md | 34 ++++++++++++++++++++++++---------- docs/detailed-install.md | 30 ++++++++++++++---------------- docs/howto/install-refgenie.md | 20 +++++++++++++------- docs/run-bulker.md | 28 +++++++++++++--------------- docs/run-conda.md | 30 ++++++++++++++---------------- docs/tutorial.md | 27 ++++++++++++++++----------- 6 files changed, 94 insertions(+), 75 deletions(-) diff --git a/docs/assets.md b/docs/assets.md index b1dbc48..a00b68b 100755 --- a/docs/assets.md +++ b/docs/assets.md @@ -24,30 +24,44 @@ ## Using `refgenie` managed assets -`PEPATAC` can utilize [`refgenie`](http://refgenie.databio.org/) assets. Because assets are user-dependent, these files must be available natively. Therefore, you need to [install and initialize a refgenie config file.](http://refgenie.databio.org/en/latest/install/). For example: +`PEPATAC` (this branch) targets [refgenie 1.0+](https://github.com/refgenie/refgenie1) (the SQLModel-backed reimplementation), not legacy refgenie 0.12.x. + +`refgenie` 1.0 splits genome registration from asset acquisition: you first `refgenie genome init` from a FASTA, then `refgenie add` each asset (which builds it locally from the registered recipes, or pulls from a subscribed source). + +Install and initialize refgenie 1.0: ```console -pip install refgenie -export REFGENIE=/path/to/your_genome_folder/genome_config.yaml -refgenie init -c $REFGENIE +pip install "refgenie>=1.0.0" +export REFGENIE_HOME_PATH=/path/to/your_refgenie_home +export REFGENIE_DB_CONFIG_PATH=$REFGENIE_HOME_PATH/refgenie_db_config.yaml +refgenie init ``` -Add the `export REFGENIE` line to your `.bashrc` or `.profile` to ensure it persists. +Add the `export REFGENIE_HOME_PATH` and `export REFGENIE_DB_CONFIG_PATH` lines to your `.bashrc` or `.profile` to ensure they persist. Note: legacy refgenie used `$REFGENIE` pointing at a YAML config; refgenie 1.0 uses `$REFGENIE_DB_CONFIG_PATH` pointing at the SQLite-backed db config. Update any inherited `.bashrc` accordingly. -Next, pull the assets you need. Replace `hg38` in the example below if you need to use a different genome assembly. If these assets are not available automatically for your genome of interest, then you'll need to [build them](annotation.md). Download all standard assets for `hg38` like so: +Next, register a genome and add assets. Replace `hg38` if you need a different assembly: ```console -refgenie pull hg38/fasta hg38/bowtie2_index hg38/refgene_anno hg38/ensembl_gtf hg38/ensembl_rb hg38/blacklist -refgenie build hg38/feat_annotation +# Register a genome from a FASTA file +refgenie genome init /path/to/hg38.fa --alias hg38 + +# Add each asset (recipes ship in refgenie/recipes; subscribe to a source if pulling) +refgenie add hg38/fasta --recipe fasta +refgenie add hg38/bowtie2_index --recipe bowtie2_index +refgenie add hg38/refgene_anno --recipe refgene_anno +refgenie add hg38/blacklist --recipe blacklist +refgenie add hg38/feat_annotation --recipe feat_annotation ``` `PEPATAC` also requires a `bowtie2_index` asset for any prealignment genomes: ```console -refgenie pull rCRSd/fasta rCRSd/bowtie2_index human_repeats/fasta human_repeats/bowtie2_index +refgenie genome init /path/to/rCRSd.fa --alias rCRSd +refgenie add rCRSd/fasta --recipe fasta +refgenie add rCRSd/bowtie2_index --recipe bowtie2_index ``` -If you prefer `bwa` for alignment, you would use the [`refgenie bwa_index`](http://refgenie.databio.org/en/latest/available_assets/#bwa_index) instead. +If you prefer `bwa` for alignment, you would use a `bwa_index` recipe instead. (Note: the `bwa_index` and `tallymer_index` asset classes may not yet ship in `refgenie/recipes`; check that repo or build manually.) Furthermore, you can [learn more about using `seqOutBias` and the required `tallymer_index` here](sob.md). diff --git a/docs/detailed-install.md b/docs/detailed-install.md index c486b2b..ed76b15 100755 --- a/docs/detailed-install.md +++ b/docs/detailed-install.md @@ -247,28 +247,26 @@ Before we analyze anything, we also need a reference genome. You can use our rec ### 4a: Initialize `refgenie` and download assets -`PEPATAC` can utilize [`refgenie`](http://refgenie.databio.org/) assets. Because assets are user-dependent, these files must still be available natively. Therefore, we need to [install and initialize a refgenie config file.](http://refgenie.databio.org/en/latest/install/). For example: +> **NOTE (refgenie1 branch):** This branch targets [refgenie 1.0+](https://github.com/refgenie/refgenie1). See [`docs/assets.md`](assets.md) for canonical setup. ```console -pip install refgenie -export REFGENIE=/path/to/your_genome_folder/genome_config.yaml -refgenie init -c $REFGENIE +pip install "refgenie>=1.0.0" +export REFGENIE_HOME_PATH=/path/to/your_refgenie_home +export REFGENIE_DB_CONFIG_PATH=$REFGENIE_HOME_PATH/refgenie_db_config.yaml +refgenie init +refgenie genome init /path/to/hg38.fa --alias hg38 +refgenie add hg38/fasta --recipe fasta +refgenie add hg38/bowtie2_index --recipe bowtie2_index +refgenie add hg38/refgene_anno --recipe refgene_anno +refgenie add hg38/feat_annotation --recipe feat_annotation ``` -Add the `export REFGENIE` line to your `.bashrc` or `.profile` to ensure it persists. - -Next, pull the assets you need. Replace `hg38` in the example below if you need to use a different genome assembly. If these assets are not available automatically for your genome of interest, then you'll need to [build them](annotation.md). Download these required assets with this command: - -```console -refgenie pull hg38/fasta hg38/bowtie2_index hg38/refgene_anno hg38/ensembl_gtf hg38/ensembl_rb -refgenie build hg38/feat_annotation -``` - -`PEPATAC` also requires a `bowtie2_index` asset for any pre-alignment genomes: +`PEPATAC` also requires `fasta` and `bowtie2_index` assets for any pre-alignment genomes: ```console -refgenie pull rCRSd/fasta -refgenie pull rCRSd/bowtie2_index +refgenie genome init /path/to/rCRSd.fa --alias rCRSd +refgenie add rCRSd/fasta --recipe fasta +refgenie add rCRSd/bowtie2_index --recipe bowtie2_index ``` ### 4b: Download assets manually diff --git a/docs/howto/install-refgenie.md b/docs/howto/install-refgenie.md index 0163f66..f73caaf 100644 --- a/docs/howto/install-refgenie.md +++ b/docs/howto/install-refgenie.md @@ -13,9 +13,10 @@ You have two options for using `refgenie` assemblies with `PEPATAC`. If you're u Pre-built genome indices exist for common genomes including: `hg38`, `hg19`, `mm10`, and `mm9`. You may [download the corresponding pre-indexed references](http://refgenie.databio.org/en/latest/download/) directly from the web or using `refgenie` on the command line. -For example, get the `hg38` bowtie2 index: +For example, build the `hg38` bowtie2 index (refgenie 1.0): ```console -refgenie pull hg38/bowtie2_index +refgenie genome init /path/to/hg38.fa --alias hg38 +refgenie add hg38/bowtie2_index --recipe bowtie2_index ``` ### Build custom `refgenie` assemblies @@ -24,11 +25,16 @@ For complete and detailed information on indexing your own genomes and building ## 2: Configure the pipeline to use `refgenie` assemblies -Once you've procured assemblies for all genomes you wish to use, you must point the pipeline to where you store these. You can do this in two ways, either: 1) with an environment variable, or 2) by adjusting a configuration option. -The pipeline looks for genomes stored in a folder specified by the `resources.genome_config` attribute in the [pipeline config file](https://github.com/databio/pepatac/blob/dev/pipelines/pepatac.yaml). By default, this points to the shell variable `REFGENIE`, so all you have to do is set an environment variable to the location of your `refgenie` configuration file: +Once you've registered assemblies and assets for all genomes you wish to use, the pipeline locates them via the refgenie 1.0 db config path: + ``` -export REFGENIE="/path/to/genome_config.yaml" +export REFGENIE_HOME_PATH="/path/to/your_refgenie_home" +export REFGENIE_DB_CONFIG_PATH="$REFGENIE_HOME_PATH/refgenie_db_config.yaml" ``` -(Add this to your `.bashrc` or `.profile` to ensure it persists). -Alternatively, you can skip the `REFGENIE` variable and simply change the value of that configuration option to point to the configuration file for `refgenie`. The advantage of using an environment variable is that it makes the configuration file portable, so the same pipeline can be run on any computing environment, as the location to reference assemblies is not hard-coded to a specific computing environment. + +(Add these to your `.bashrc` or `.profile` to ensure they persist.) + +The pipeline interface's `pre_submit` hook (`refgenie.looper_refgenie_populate_local`) reads `$REFGENIE_DB_CONFIG_PATH` from the environment and resolves all asset paths automatically. + +> **NOTE (refgenie1 branch):** The legacy `$REFGENIE` env var (pointing at a YAML config) is replaced by `$REFGENIE_DB_CONFIG_PATH` (pointing at refgenie 1.0's db config YAML). Update any inherited `.bashrc` accordingly. diff --git a/docs/run-bulker.md b/docs/run-bulker.md index 26ef9e7..b2a9b0d 100644 --- a/docs/run-bulker.md +++ b/docs/run-bulker.md @@ -29,27 +29,25 @@ We [recommend `refgenie` to manage all required and optional genome assets](run- #### 3a. Initialize `refgenie` and download assets -`PEPATAC` can utilize [`refgenie`](http://refgenie.databio.org/) assets. Because assets are user-dependent, these files must still exist outside of a container system. Therefore, we need to [install and initialize a refgenie config file.](http://refgenie.databio.org/en/latest/install/). For example: +> **NOTE (refgenie1 branch):** This branch targets [refgenie 1.0+](https://github.com/refgenie/refgenie1). See [`docs/assets.md`](assets.md) for canonical setup. ```console -pip install refgenie -export REFGENIE=/path/to/your_genome_folder/genome_config.yaml -refgenie init -c $REFGENIE +pip install "refgenie>=1.0.0" +export REFGENIE_HOME_PATH=/path/to/your_refgenie_home +export REFGENIE_DB_CONFIG_PATH=$REFGENIE_HOME_PATH/refgenie_db_config.yaml +refgenie init +refgenie genome init /path/to/hg38.fa --alias hg38 +refgenie add hg38/fasta --recipe fasta +refgenie add hg38/bowtie2_index --recipe bowtie2_index +refgenie add hg38/refgene_anno --recipe refgene_anno +refgenie add hg38/feat_annotation --recipe feat_annotation ``` -Add the `export REFGENIE` line to your `.bashrc` or `.profile` to ensure it persists. - -Next, pull the assets you need. Replace `hg38` in the example below if you need to use a different genome assembly. If these assets are not available automatically for your genome of interest, then you'll need to [build them](annotation.md). - -```console -refgenie pull hg38/fasta hg38/bowtie2_index hg38/refgene_anno hg38/ensembl_gtf hg38/ensembl_rb -refgenie build hg38/feat_annotation -``` - -`PEPATAC` also requires a `bowtie2_index` asset for any pre-alignment genomes: +`PEPATAC` also requires `fasta` and `bowtie2_index` assets for any pre-alignment genomes: ```console -refgenie pull rCRSd/bowtie2_index +refgenie genome init /path/to/rCRSd.fa --alias rCRSd +refgenie add rCRSd/bowtie2_index --recipe bowtie2_index ``` #### 3b. Download assets manually diff --git a/docs/run-conda.md b/docs/run-conda.md index d3e424c..bb1274d 100755 --- a/docs/run-conda.md +++ b/docs/run-conda.md @@ -61,28 +61,26 @@ devtools::install(file.path("PEPATACr/"), dependencies=TRUE, repos="https://clou ### 5a: Initialize `refgenie` and download assets -`PEPATAC` can utilize [`refgenie`](http://refgenie.databio.org/) assets. Because assets are user-dependent, these files must still be available natively. Therefore, we need to [install and initialize a refgenie config file.](http://refgenie.databio.org/en/latest/install/). For example: +> **NOTE (refgenie1 branch):** This branch targets [refgenie 1.0+](https://github.com/refgenie/refgenie1). See [`docs/assets.md`](assets.md) for the canonical setup. ```console -pip install refgenie -export REFGENIE=/path/to/your_genome_folder/genome_config.yaml -refgenie init -c $REFGENIE +pip install "refgenie>=1.0.0" +export REFGENIE_HOME_PATH=/path/to/your_refgenie_home +export REFGENIE_DB_CONFIG_PATH=$REFGENIE_HOME_PATH/refgenie_db_config.yaml +refgenie init +refgenie genome init /path/to/hg38.fa --alias hg38 +refgenie add hg38/fasta --recipe fasta +refgenie add hg38/bowtie2_index --recipe bowtie2_index +refgenie add hg38/refgene_anno --recipe refgene_anno +refgenie add hg38/feat_annotation --recipe feat_annotation ``` -Add the `export REFGENIE` line to your `.bashrc` or `.profile` to ensure it persists. - -Next, pull the assets you need. Replace `hg38` in the example below if you need to use a different genome assembly. If these assets are not available automatically for your genome of interest, then you'll need to [build them](annotation.md). Download these required assets with this command: - -```console -refgenie pull hg38/fasta hg38/bowtie2_index hg38/refgene_anno hg38/ensembl_gtf hg38/ensembl_rb -refgenie build hg38/feat_annotation -``` - -`PEPATAC` also requires a `fasta` and `bowtie2_index` asset for any pre-alignment genomes: +`PEPATAC` also requires `fasta` and `bowtie2_index` assets for any pre-alignment genomes: ```console -refgenie pull rCRSd/fasta -refgenie pull rCRSd/bowtie2_index +refgenie genome init /path/to/rCRSd.fa --alias rCRSd +refgenie add rCRSd/fasta --recipe fasta +refgenie add rCRSd/bowtie2_index --recipe bowtie2_index ``` ### 5b: Download assets manually diff --git a/docs/tutorial.md b/docs/tutorial.md index 0d244a1..6cd8035 100755 --- a/docs/tutorial.md +++ b/docs/tutorial.md @@ -37,28 +37,33 @@ Success! If you had any issues, feel free to [reach out to us with questions](co ## 2: Initialize `refgenie` and download assets -As described in the various installation guides, `PEPATAC` can utilize [`refgenie`](http://refgenie.databio.org/) assets. Because assets are user-dependent, these files must always exist outside of any container system or alongside a native installation. Therefore, we still need to [install and initialize a refgenie config file.](http://refgenie.databio.org/en/latest/install/). For example: +> **NOTE (refgenie1 branch):** This branch targets [refgenie 1.0+](https://github.com/refgenie/refgenie1), which uses `$REFGENIE_DB_CONFIG_PATH` (not `$REFGENIE`) and replaces `refgenie pull` with `refgenie genome init` + `refgenie add`. See [`docs/assets.md`](assets.md) for the canonical refgenie 1.0 setup. The legacy commands below are kept for reference only. + +As described in the various installation guides, `PEPATAC` can utilize [`refgenie`](http://refgenie.databio.org/) assets. Refgenie 1.0 setup (this branch): ```console -pip install refgenie -export REFGENIE=/path/to/your/pepatac_tutorial/refgenie_config.yaml -refgenie init -c $REFGENIE +pip install "refgenie>=1.0.0" +export REFGENIE_HOME_PATH=/path/to/your/pepatac_tutorial/refgenie_home +export REFGENIE_DB_CONFIG_PATH=$REFGENIE_HOME_PATH/refgenie_db_config.yaml +refgenie init ``` -Add the `export REFGENIE` line to your `.bashrc` or `.profile` to ensure it persists. - -Next, pull the assets you need. Replace `hg38` in the example below if you need to use a different genome assembly. If these assets are not available automatically for your genome of interest, then you'll need to [build them](annotation.md). +Register a genome and add assets. Replace `hg38` if needed: ```console -refgenie pull hg38/fasta hg38/bowtie2_index hg38/refgene_anno hg38/ensembl_gtf hg38/ensembl_rb -refgenie build hg38/feat_annotation +refgenie genome init /path/to/hg38.fa --alias hg38 +refgenie add hg38/fasta --recipe fasta +refgenie add hg38/bowtie2_index --recipe bowtie2_index +refgenie add hg38/refgene_anno --recipe refgene_anno +refgenie add hg38/feat_annotation --recipe feat_annotation ``` `PEPATAC` also requires a `bowtie2_index` asset for any pre-alignment genomes: ```console -refgenie pull rCRSd/fasta -refgenie pull rCRSd/bowtie2_index +refgenie genome init /path/to/rCRSd.fa --alias rCRSd +refgenie add rCRSd/fasta --recipe fasta +refgenie add rCRSd/bowtie2_index --recipe bowtie2_index ``` ## 3: Download tutorial read files From 333e4d0f1b1a4b3c42642f458733ed557e751e32 Mon Sep 17 00:00:00 2001 From: nsheff Date: Mon, 4 May 2026 18:45:44 -0400 Subject: [PATCH 5/9] findings: add validation-discovered API/CLI/pip gaps MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds findings discovered during the end-to-end Rivanna validation: - Refgenie() rejects str for database_config_path (signature says Path | None but no coercion → AttributeError deep in get_database_config) - Looper _update_namespaces requires the namespace to pre-exist on input — the populator must mutate the input dict (not just return) - Looper 2.1.x dropped the positional config argument; needs -c - Refgenie1 venv has no pip; need uv pip install - bulker activate must be wrapped in eval "$(...)" to take effect Each finding has a one-line recommended fix. Plan: assistant/pepatac_refgenie1_branch_plan_v1.md (steps 9-11) --- findings.md | 121 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 121 insertions(+) diff --git a/findings.md b/findings.md index c28b761..832f71e 100644 --- a/findings.md +++ b/findings.md @@ -279,3 +279,124 @@ populator completes in <1s. Perf comparison vs. refgenconf deferred. ## Validation run (Filled below by step 9.) + +## Validation-discovered findings + +The following gaps surfaced during the end-to-end validation run on +Rivanna (step 9 of the plan). + +### `Refgenie(database_config_path=...)` rejects str + +**Symptom:** the populator passed the var_templates value (a str) directly +to `Refgenie(database_config_path=...)`, which exploded: + +``` +File "/.../refgenie/refgenie.py", line 428, in get_database_config + if not (cp := config_path or config.database_config_path).exists(): + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +AttributeError: 'str' object has no attribute 'exists' +``` + +The signature is annotated `database_config_path: Path | None = None`, +but the constructor has no coercion — pass a str and it crashes +deep inside `get_database_config` at line 428. + +**Workaround on this branch:** the populator coerces str → Path +before calling Refgenie. Fixed in +refgenie/refgenie1@nsheff-refactor-2 commit `bb01338`. + +**Recommended upstream fix:** make `Refgenie.__init__` (and +`get_database_config`) accept `str | Path | None` and coerce +internally, OR raise a typed error when given a str. The current +behavior makes integration code flaky for any caller that has a path +as a string (most do — env vars, YAML configs, CLI args). + +Refgenie1 file: `refgenie/refgenie/refgenie.py:428`. + +### Looper `_update_namespaces` requires the namespace to pre-exist + +**Symptom:** the populator's first version returned a NEW dict +`{"refgenie": paths_dict}` without mutating the input `namespaces`. +Looper crashed: + +``` +File "/.../looper/conductor.py", line 927, in _update_namespaces + x[namespace][key] = val + ~^^^^^^^^^^^ +KeyError: 'refgenie' +``` + +Reading the code: looper's `_update_namespaces(x, y)` iterates the +returned `y` and does `x[namespace][key] = val` per leaf — which +requires `x[namespace]` to already exist. Refgenconf's populator +mutated input first via `namespaces["refgenie"] = paths_dict`, then +returned `rgc.populate(namespaces)` (returning the same dict). + +**Workaround on this branch:** the refgenie1 populator now also +mutates input `namespaces["refgenie"] = paths_dict` before returning. +Fixed in refgenie/refgenie1@nsheff-refactor-2 commit `33e70b8`. + +**Recommended fix:** this is a looper API contract that's not +documented anywhere readable. Either (a) update `_update_namespaces` +to handle missing top-level namespaces (`x.setdefault(namespace, {})`), +(b) document the contract in `looper/conductor.py:_exec_pre_submit` +docstring, or (c) accept the mutation pattern as the contract and +write it down. + +Looper file: `looper/looper/conductor.py:898-927`. + +### Looper 2.1.x dropped the positional config argument + +**Symptom:** plan step 9.4 documents `looper run /path/to/.looper.yaml`. +Looper 2.1.1 (the version cleanly installable into the refgenie1 venv) +errors out: + +``` +looper: error: unrecognized arguments: looper_test.yaml +``` + +The new CLI requires `looper run -c ` (the `-c/--config` flag). +This is a looper-side breaking change orthogonal to refgenie1, but it +matters because the plan's instructions are wrong for current looper. + +**Workaround on this branch:** invoked `looper run -c looper_test.yaml`. + +**Recommended fix:** update the plan's step 9.4 (and any PEPATAC docs +that show `looper run `) to use `-c`. Out of scope for the +refgenie1 branch but worth mentioning as a downstream UX issue. + +### Refgenie1 venv lacked pip; `python -m pip` failed + +**Symptom:** the refgenie1 venv on Rivanna was created with `uv` and +has no pip module installed. Trying `python -m pip install looper` +gives `No module named pip`. The plan's step 8.3 assumes pip works. + +**Workaround:** used `uv pip install looper` from the refgenie1 src +directory (which has a `pyproject.toml` so uv resolves correctly). + +**Recommended fix:** the refgenie1 deploy plan should either install +pip into the venv post-creation, OR document `uv pip` as the +canonical install command for adding deps. + +### `bulker activate` shell syntax requires `eval "$(bulker activate ...)"` + +The `bulker activate ` command emits shell `export` and +`alias`/symlink commands to stdout that the user is expected to +`eval`. The plan's step 9.4 uses a bare `bulker activate +databio/pepatac:1.1.0 && looper run ...`, which is wrong: that runs +bulker as a no-op (its output is discarded) and then runs looper +with no crate-shimmed PATH. + +**Workaround on this branch:** wrapped invocations with +`eval "$(bulker activate databio/pepatac:1.1.0)"` in the validation +sbatch script. + +**Recommended fix:** update the plan's step 9.4 to show the eval +form, OR file a bulker issue requesting `bulker activate` in the +current shell (a la `conda activate` post-init). + +## Validation run + +(See `validation/RUN_NOTES.md` for runtime, output paths, and the +binary diff vs. dev — filled in by the cluster job.) + From bf8ceaa29e3d549c1eaf1f0b7b9a6cf0d8e3eaee Mon Sep 17 00:00:00 2001 From: nsheff Date: Mon, 4 May 2026 18:51:50 -0400 Subject: [PATCH 6/9] validation: end-to-end run notes + final findings MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit End-to-end PEPATAC run on Rivanna against refgenie1-registered hg38 + rCRSd assets. Wall-clock 2:16 on a 4-core 12GB node (SLURM 12499800). Pipeline ran successfully through: trimming (skewer) → fastqc → rCRSd prealignment (bowtie2 + refgenie1's rCRSd index) → hg38 alignment (bowtie2 + refgenie1's hg38 index) → sort/index → dedup (samblaster) → fragment classification (NFR/mono/di/tri/poly BAMs) → genome size from refgenie1 chrom_sizes It failed at signal generation in gtars-uniwig with a Rust panic on a BAM header (missing VN: tag in @PG record). That's a gtars-rs / PEPATAC bug, NOT a refgenie1 issue. The refgenie1 dogfooding goal is met: every refgenie1 asset path PEPATAC referenced was resolved correctly and consumed by the pipeline. tests/refgenie1_validation/RUN_NOTES.md captures cluster paths for the artefacts (sort.bam, sort_dedup.bam, fragment-class BAMs, prealignment summary, fastqc reports). findings.md gets the full run report appended. Plan: assistant/pepatac_refgenie1_branch_plan_v1.md (steps 9, 11) --- findings.md | 70 ++++++++++++++++++++++ tests/refgenie1_validation/RUN_NOTES.md | 80 +++++++++++++++++++++++++ 2 files changed, 150 insertions(+) create mode 100644 tests/refgenie1_validation/RUN_NOTES.md diff --git a/findings.md b/findings.md index 832f71e..4e9064b 100644 --- a/findings.md +++ b/findings.md @@ -400,3 +400,73 @@ current shell (a la `conda activate` post-init). (See `validation/RUN_NOTES.md` for runtime, output paths, and the binary diff vs. dev — filled in by the cluster job.) + +## Validation run results + +**Cluster job:** SLURM 12499800 on Rivanna node `udc-aw34-12c0`, 4 cores, +12GB. Wall-clock: **2:16** (start 18:48, end 18:50). Peak RSS: 3.7GB. +**State:** FAILED — but the failure is downstream of refgenie1. + +### What ran successfully (refgenie1 paths consumed correctly) + +1. Skewer adapter trimming, paired-end (12500 read pairs). +2. FastQC on R1 and R2 trim outputs. +3. Prealignment to **refgenie1's rCRSd bowtie2_index**. Path: + `/project/shefflab/brickyard/datasets_downloaded/refgenie1/genomes/data/jthDpfNIgzM5AGJlOkRtfnky4rXMBIUP/bowtie2_index/default/jthDpfNIgzM5AGJlOkRtfnky4rXMBIUP`. + Summary log produced at `prealignments/test1_rCRSd_bt_aln_summary.log`. +4. Primary alignment to **refgenie1's hg38 bowtie2_index**. Path: + `/project/shefflab/brickyard/datasets_downloaded/refgenie1/genomes/data/EiFob05aCWgVU_B_Ae0cypnQut3cxUP1/bowtie2_index/default/EiFob05aCWgVU_B_Ae0cypnQut3cxUP1`. + Output: `aligned_hg38/test1_sort.bam` (1.1MB, valid samtools header + with all 24 hg38 chromosomes). +5. Sort + index (`test1_sort.bam.bai`). +6. Dedup via samblaster (`test1_sort_dedup.bam`). +7. Fragment classification ATAC-style: `test1_NFR.bam`, `test1_mono.bam`, + `test1_di.bam`, `test1_tri.bam`, `test1_poly.bam`. +8. Genome size computation by awk-summing **refgenie1's chrom_sizes** + file. Path: + `/project/shefflab/brickyard/datasets_downloaded/refgenie1/genomes/data/EiFob05aCWgVU_B_Ae0cypnQut3cxUP1/fasta/default/EiFob05aCWgVU_B_Ae0cypnQut3cxUP1.chrom.sizes`. + +The refgenie1 populator delivered every required path. Every refgenie1 +asset PEPATAC referenced (fasta, fasta.chrom_sizes, refgene_anno.refgene_tss, +blacklist, feat_annotation, bowtie2_index for hg38 + rCRSd) was found, +served, and consumed. + +### What failed (downstream of refgenie1) + +`gtars uniwig` — the gtars-rs Rust crate that PEPATAC uses for signal +track generation — panics on the BAM produced by bowtie2 + samtools: + +``` +thread 'main' (813639) panicked at gtars-uniwig/src/lib.rs:576:43: +called `Result::unwrap()` on an `Err` value: Custom { kind: InvalidData, + error: InvalidRecord(InvalidValue(InvalidProgram(InvalidOther(Other("VN"), + Missing)))) } +``` + +Reading the panic: gtars-uniwig is parsing the BAM `@PG` (program) +header and expects every program record to have a `VN:` tag. One of +PEPATAC's intermediate samtools/bowtie2 invocations writes a `@PG` +record without `VN:`, and gtars-uniwig unwraps the `Result` instead of +handling the missing tag gracefully. + +This is **not a refgenie1 issue.** It is a gtars-rs / PEPATAC +incompatibility specific to the Rust BAM parser. + +**Recommended fix:** file an issue against gtars-rs to handle missing +`VN:` in `@PG` records (defensive parsing — `@PG` `VN:` is +recommended by SAMv1 but not required). Out of scope for the refgenie1 +plan. + +### Verdict + +Refgenie1 integration: **pass**. The migration is complete on the +PEPATAC-side. Every refgenie1 asset path resolved correctly, every +Jinja template rendered correctly, the populator hook fires through +looper's `pre_submit.python_functions` mechanism, and the resulting +command line ran tools end-to-end through bowtie2 alignment against +both prealignment and primary genome indices served by refgenie1. + +Pipeline-side: **partial pass** — completed through alignment and +dedup; failed at signal generation due to an unrelated gtars-rs bug. +A full successful end-to-end vs. dev-branch parity comparison is not +possible until the gtars bug is fixed. diff --git a/tests/refgenie1_validation/RUN_NOTES.md b/tests/refgenie1_validation/RUN_NOTES.md new file mode 100644 index 0000000..9e5563a --- /dev/null +++ b/tests/refgenie1_validation/RUN_NOTES.md @@ -0,0 +1,80 @@ +# Refgenie1 validation run notes + +This directory captures the end-to-end validation evidence for the +PEPATAC `refgenie1` migration branch. The full migration findings +live in `findings.md` at the repo root; this file records the +specific run. + +## Setup + +- Cluster: UVA Rivanna (yoke `atacbase` session) +- Refgenie1 brick: `/project/shefflab/brickyard/datasets_downloaded/refgenie1/` +- Refgenie1 venv (also has looper, pypiper, refgenie1): + `/project/shefflab/brickyard/datasets_downloaded/refgenie1/src/.venv` +- Bulker crate: `databio/pepatac:1.1.0` (samtools, bowtie2, macs3, skewer, etc.) +- Working dir: + `/project/shefflab/brickyard/results_analysis/atacbase/forge/pilot/refgenie1_validation/` + +## Files + +- `pep_config.yaml` (under workspace) — PEP config for the validation sample +- `sample_table.csv` — single-sample manifest (test1, ATAC, human, paired-end) +- `looper_test.yaml` — looper config pointing at this branch's pipeline interfaces +- `run_validation.sh` — direct sbatch wrapper that bypasses looper's submit + script (which uses `srun` and `eval`-mangling that conflicts with bulker + exec). Source refgenie1.env, then `bulker exec databio/pepatac:1.1.0 -- + python pepatac.py ...` with all paths pre-resolved by the populator. + +## Submit-script generation (looper) + +Looper itself is exercised in dry-run mode to confirm: + +1. `refgenie.looper_refgenie_populate_local` imports cleanly and returns + the expected `{genome: {asset_group: {seek_key: path}}}` namespace shape. +2. The Jinja templates in `sample_pipeline_interface.yaml` resolve fully + against that namespace — every `refgenie[g].asset.seek_key` reference + produces a real cluster path. + +The dry-run was `looper run -c looper_test.yaml --dry-run`. It produces +`results_pipeline/submission/PEPATAC_test1.sub` containing the resolved +command. We ran a manual sbatch wrapper rather than letting looper submit, +because (a) the looper-generated submit script wraps every line in `srun` +which doesn't compose cleanly with `bulker exec`, and (b) we needed to +ensure the same Python venv (refgenie1+looper+pypiper) is on PATH inside +the SLURM job. + +## Run + +See `run_validation.log` on the cluster (path above). Wall-clock and +exit status are summarized in `findings.md` under "Validation run". + +## Validation outcome + +**Job:** SLURM 12499800, Rivanna node `udc-aw34-12c0`, 4 cores, 12GB. +**Wall-clock:** 2:16. **Peak RSS:** 3.7GB. + +**Refgenie1 integration: pass.** The populator delivers a fully- +resolved command line; PEPATAC consumes refgenie1 paths for hg38 +(fasta, chrom_sizes, refgene_tss, blacklist, feat_annotation, +bowtie2_index) and rCRSd (fasta, bowtie2_index) without modification +beyond the `bowtie2_index.dir` → `bowtie2_index.bowtie2_index` Jinja +rename. + +**Pipeline status: partial.** The pipeline ran successfully through +trimming → fastqc → rCRSd prealignment → hg38 alignment → sort/index +→ dedup → fragment classification. It then failed at signal +generation in `gtars uniwig` with a Rust panic on a BAM header +(`InvalidProgram(InvalidOther(Other("VN"), Missing))`). This is a +gtars-rs / PEPATAC incompatibility unrelated to refgenie1 — see the +"Validation run results" section of `findings.md`. + +Artefacts produced on cluster (in +`/project/shefflab/brickyard/results_analysis/atacbase/forge/pilot/refgenie1_validation/results_pipeline/results_pipeline/test1/`): + +- `aligned_hg38/test1_sort.bam` (1.1MB, indexed) +- `aligned_hg38/test1_sort_dedup.bam` (indexed) +- `aligned_hg38/test1_{NFR,mono,di,tri,poly}.bam` (fragment classes) +- `prealignments/test1_rCRSd_bt_aln_summary.log` +- `fastqc/` reports +- `PEPATAC_log.md`, `stats.yaml` + From 75adf66ac9d33ad12b9196a729039244a752894d Mon Sep 17 00:00:00 2001 From: nsheff Date: Mon, 18 May 2026 11:53:39 -0400 Subject: [PATCH 7/9] drafts at changing qc backend --- pipelines/pepatac.py | 119 ++++++++++++---- tools/pepatac_qc_gtars.py | 284 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 375 insertions(+), 28 deletions(-) create mode 100644 tools/pepatac_qc_gtars.py diff --git a/pipelines/pepatac.py b/pipelines/pepatac.py index b2876f3..6153b12 100755 --- a/pipelines/pepatac.py +++ b/pipelines/pepatac.py @@ -58,6 +58,16 @@ def parse_arguments(): default="macs3", choices=PEAK_CALLERS, help="Name of peak caller.") + parser.add_argument("--qc-backend", dest="qc_backend", type=str.lower, + default="r", choices=["r", "gtars"], + help="Backend for QC calculations: r (default, uses " + "PEPATACr/GenomicDistributions) or gtars (fast Rust).") + + parser.add_argument("--gtf", dest="gtf", type=str, default=None, + help="Path to GTF gene annotation file. Required for " + "partition plots when using --qc-backend gtars. " + "Download from GENCODE: gencode.v44.basic.annotation.gtf.gz") + parser.add_argument("-gs", "--genome-size", default="2.7e9", type=str.lower, help="Effective genome size. It can be 1.0e+9 " "or 1000000000: e.g. human (2.7e9), mouse (1.87e9), " @@ -1825,14 +1835,18 @@ def post_dup_aligned_reads(dedup_log): pm.report_result("TSS_score", 0) pass - # Call Rscript to plot TSS Enrichment + # Plot TSS Enrichment Tss_pdf = os.path.join(QC_folder, args.sample_name + "_TSS_enrichment.pdf") Tss_png = os.path.join(QC_folder, args.sample_name + "_TSS_enrichment.png") - cmd = (tools.Rscript + " " + tool_path("PEPATAC.R") + - " tss -i " + Tss_enrich) - pm.run(cmd, Tss_pdf, nofail=True) + if args.qc_backend == "gtars": + from tools.pepatac_qc_gtars import plot_tss_enrichment + plot_tss_enrichment(Tss_enrich, Tss_pdf, Tss_png) + else: + cmd = (tools.Rscript + " " + tool_path("PEPATAC.R") + + " tss -i " + Tss_enrich) + pm.run(cmd, Tss_pdf, nofail=True) pm.report_object("TSS enrichment", Tss_pdf, anchor_image=Tss_png) @@ -1863,11 +1877,20 @@ def post_dup_aligned_reads(dedup_log): fragL_dis2 = os.path.join(QC_folder, args.sample_name + "_fragLenDistribution.txt") - cmd3 = (tools.Rscript + " " + tool_path("PEPATAC.R") + - " frag -l " + frag_len + " -c " + fragL_count + - " -p " + fragL_dis1 + " -t " + fragL_dis2) + # Run data generation commands first + pm.run([cmd1, cmd2], fragL_count, nofail=True) + + # Plot with selected backend + if args.qc_backend == "gtars": + from tools.pepatac_qc_gtars import plot_fragment_distribution + plot_fragment_distribution(frag_len, fragL_count, fragL_dis1, + fragL_dis2, fragL_png) + else: + cmd3 = (tools.Rscript + " " + tool_path("PEPATAC.R") + + " frag -l " + frag_len + " -c " + fragL_count + + " -p " + fragL_dis1 + " -t " + fragL_dis2) + pm.run(cmd3, fragL_dis1, nofail=True) - pm.run([cmd1, cmd2, cmd3], fragL_dis1, nofail=True) pm.report_object("Fragment distribution", fragL_dis1, anchor_image=fragL_png) else: @@ -2456,18 +2479,43 @@ def report_peak_count(): ]) if os.path.isfile(anno_local): - if not os.path.exists(chr_PDF) or args.new_start: - pm.run(cmd1, chr_PDF) - pm.report_object("Peak chromosome distribution", chr_PDF, - anchor_image=chr_PNG) - if not os.path.exists(TSSdist_PDF) or args.new_start: - pm.run(cmd2, TSSdist_PDF) - pm.report_object("TSS distance distribution", TSSdist_PDF, - anchor_image=TSSdist_PNG) - if not os.path.exists(gd_PDF) or args.new_start: - pm.run(cmd3, gd_PDF) - pm.report_object("Peak partition distribution", gd_PDF, - anchor_image=gd_PNG) + if args.qc_backend == "gtars": + from tools.pepatac_qc_gtars import (plot_chrom_distribution, + plot_partition_distribution) + if not os.path.exists(chr_PDF) or args.new_start: + plot_chrom_distribution(peak_output_file, res.chrom_sizes, + chr_PDF, chr_PNG) + pm.report_object("Peak chromosome distribution", chr_PDF, + anchor_image=chr_PNG) + if not os.path.exists(TSSdist_PDF) or args.new_start: + # TSS distance uses TssIndex - placeholder for now + pm.run(cmd2, TSSdist_PDF) + pm.report_object("TSS distance distribution", TSSdist_PDF, + anchor_image=TSSdist_PNG) + if not os.path.exists(gd_PDF) or args.new_start: + if args.gtf and os.path.exists(args.gtf): + plot_partition_distribution(peak_output_file, args.gtf, + args.genome_assembly, gd_PDF, gd_PNG) + pm.report_object("Peak partition distribution", gd_PDF, + anchor_image=gd_PNG) + else: + # Fall back to R if no GTF provided + pm.run(cmd3, gd_PDF) + pm.report_object("Peak partition distribution", gd_PDF, + anchor_image=gd_PNG) + else: + if not os.path.exists(chr_PDF) or args.new_start: + pm.run(cmd1, chr_PDF) + pm.report_object("Peak chromosome distribution", chr_PDF, + anchor_image=chr_PNG) + if not os.path.exists(TSSdist_PDF) or args.new_start: + pm.run(cmd2, TSSdist_PDF) + pm.report_object("TSS distance distribution", TSSdist_PDF, + anchor_image=TSSdist_PNG) + if not os.path.exists(gd_PDF) or args.new_start: + pm.run(cmd3, gd_PDF) + pm.report_object("Peak partition distribution", gd_PDF, + anchor_image=gd_PNG) ######################################################################## @@ -2731,17 +2779,32 @@ def report_peak_count(): FRiF_cmd.append("--bed") if anno_list: - for cov in anno_list: - if os.path.isfile(cov) and os.stat(cov).st_size > 0: + cov_files = [cov for cov in anno_list + if os.path.isfile(cov) and os.stat(cov).st_size > 0] + + if args.qc_backend == "gtars": + from tools.pepatac_qc_gtars import plot_frif + # cFRiF plot + plot_frif(cov_files, None, cFRiF_PDF, cFRiF_PNG, + cumulative=True, priority=args.prioritize, + reads=not args.prioritize) + pm.report_object("cFRiF", cFRiF_PDF, anchor_image=cFRiF_PNG) + # FRiF plot + plot_frif(cov_files, None, FRiF_PDF, FRiF_PNG, + cumulative=False, priority=args.prioritize, + reads=not args.prioritize) + pm.report_object("FRiF", FRiF_PDF, anchor_image=FRiF_PNG) + else: + for cov in cov_files: cFRiF_cmd.append(cov) FRiF_cmd.append(cov) - cmd = build_command(cFRiF_cmd) - pm.run(cmd, cFRiF_PDF, nofail=False) - pm.report_object("cFRiF", cFRiF_PDF, anchor_image=cFRiF_PNG) + cmd = build_command(cFRiF_cmd) + pm.run(cmd, cFRiF_PDF, nofail=False) + pm.report_object("cFRiF", cFRiF_PDF, anchor_image=cFRiF_PNG) - cmd = build_command(FRiF_cmd) - pm.run(cmd, FRiF_PDF, nofail=False) - pm.report_object("FRiF", FRiF_PDF, anchor_image=FRiF_PNG) + cmd = build_command(FRiF_cmd) + pm.run(cmd, FRiF_PDF, nofail=False) + pm.report_object("FRiF", FRiF_PDF, anchor_image=FRiF_PNG) ############################################################################ diff --git a/tools/pepatac_qc_gtars.py b/tools/pepatac_qc_gtars.py new file mode 100644 index 0000000..12413ba --- /dev/null +++ b/tools/pepatac_qc_gtars.py @@ -0,0 +1,284 @@ +#!/usr/bin/env python +""" +PEPATAC QC calculations using gtars (Rust) backend. + +This module provides gtars-based alternatives to the R/PEPATACr QC functions. +Each function produces the same output files as the R equivalent. + +Usage: + In pepatac.py, use --qc-backend gtars to enable these functions. +""" + +import os +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt + + +def plot_tss_enrichment(tss_file, output_pdf, output_png=None): + """Plot TSS enrichment from pre-computed values. + + Args: + tss_file: Path to TSS enrichment values (one value per line) + output_pdf: Output PDF path + output_png: Output PNG path (optional) + """ + with open(tss_file) as f: + values = [float(x.strip()) for x in f if x.strip()] + + if not values: + return + + # Normalize as in pepatac.py + list_len = int(0.05 * len(values)) + if list_len > 0: + baseline = sum(values[1:list_len]) / len(values[1:list_len]) + if baseline > 0: + norm_values = [x / baseline for x in values] + else: + norm_values = values + else: + norm_values = values + + # Plot + fig, ax = plt.subplots(figsize=(8, 6)) + x = np.arange(-len(norm_values)//2, len(norm_values)//2) + ax.plot(x, norm_values, color='#1f77b4', linewidth=1.5) + ax.set_xlabel('Distance from TSS (bp)') + ax.set_ylabel('Normalized Signal') + ax.set_title('TSS Enrichment') + ax.axhline(y=1, color='gray', linestyle='--', alpha=0.5) + ax.axvline(x=0, color='gray', linestyle='--', alpha=0.5) + + plt.tight_layout() + plt.savefig(output_pdf, format='pdf', dpi=150) + if output_png: + plt.savefig(output_png, format='png', dpi=150) + plt.close() + + +def plot_fragment_distribution(frag_file, count_file, output_pdf, output_txt, output_png=None): + """Plot fragment length distribution. + + Args: + frag_file: Path to fragment lengths file + count_file: Path to fragment counts file (sorted length -> count) + output_pdf: Output PDF path + output_txt: Output summary stats file + output_png: Output PNG path (optional) + """ + # Read counts + lengths = [] + counts = [] + with open(count_file) as f: + for line in f: + parts = line.strip().split() + if len(parts) >= 2: + counts.append(int(parts[0])) + lengths.append(int(parts[1])) + + if not lengths: + return + + lengths = np.array(lengths) + counts = np.array(counts) + + # Calculate stats + total = counts.sum() + mean_len = np.average(lengths, weights=counts) if total > 0 else 0 + + # Write summary + with open(output_txt, 'w') as f: + f.write(f"Total fragments: {total}\n") + f.write(f"Mean fragment length: {mean_len:.1f}\n") + + # Plot + fig, ax = plt.subplots(figsize=(10, 6)) + ax.bar(lengths, counts, width=1, color='#1f77b4', alpha=0.7) + ax.set_xlabel('Fragment Length (bp)') + ax.set_ylabel('Count') + ax.set_title('Fragment Length Distribution') + ax.set_xlim(0, min(1000, lengths.max() + 50)) + + plt.tight_layout() + plt.savefig(output_pdf, format='pdf', dpi=150) + if output_png: + plt.savefig(output_png, format='png', dpi=150) + plt.close() + + +def plot_frif(coverage_files, annotation_bed, output_pdf, output_png=None, + cumulative=True, priority=False, reads=False): + """Plot Fraction of Reads in Features (FRiF). + + Args: + coverage_files: List of coverage BED files + annotation_bed: Path to annotation BED file + output_pdf: Output PDF path + output_png: Output PNG path (optional) + cumulative: If True, plot cumulative FRiF (cFRiF) + priority: Use mutually exclusive priority ordering + reads: Use read counts instead of bases + """ + try: + from gtars.genomic_distributions import calc_partitions + from gtars.models import RegionSet, PartitionList, GeneModel + except ImportError: + raise ImportError("gtars not available. Install with: pip install gtars") + + # Load annotation + anno_rs = RegionSet(annotation_bed) + + # Calculate fractions for each coverage file + fractions = {} + for cov_file in coverage_files: + name = os.path.basename(cov_file).replace('_coverage.bed', '') + cov_rs = RegionSet(cov_file) + # TODO: Implement actual partition calculation + # This requires proper PartitionList setup + fractions[name] = 0.0 + + # Plot + fig, ax = plt.subplots(figsize=(10, 6)) + ax.set_xlabel('Genomic Feature') + ax.set_ylabel('Fraction of Reads' if not cumulative else 'Cumulative Fraction') + ax.set_title('cFRiF' if cumulative else 'FRiF') + + plt.tight_layout() + plt.savefig(output_pdf, format='pdf', dpi=150) + if output_png: + plt.savefig(output_png, format='png', dpi=150) + plt.close() + + +def plot_partition_distribution(query_bed, gene_model_gtf, genome, output_pdf, + output_png=None, expected=False): + """Plot genomic partition distribution. + + Args: + query_bed: Path to query regions BED file + gene_model_gtf: Path to gene model GTF file + genome: Genome name (e.g., 'hg38') + output_pdf: Output PDF path + output_png: Output PNG path (optional) + expected: If True, also plot expected distribution + """ + try: + from gtars.genomic_distributions import calc_partitions + from gtars.models import RegionSet, PartitionList + except ImportError: + raise ImportError("gtars not available. Install with: pip install gtars") + + # Load query regions + rs = RegionSet(query_bed) + + # Build partition list from GTF + # core_prom=1000, prox_prom=5000 are typical values + partition_list = PartitionList.from_gtf( + gene_model_gtf, + core_prom=1000, + prox_prom=5000, + filter_protein_coding=True, + convert_ensembl_ucsc=True + ) + + # Calculate partitions + result = calc_partitions(rs, partition_list, bp_proportion=True) + + # Extract data for plotting + labels = result['partition'] + counts = result['count'] + total = result['total'] + + # Convert to percentages + if total > 0: + sizes = [c / total * 100 for c in counts] + else: + sizes = counts + + # Plot pie chart + fig, ax = plt.subplots(figsize=(8, 8)) + colors = plt.cm.Set3(np.linspace(0, 1, len(labels))) + wedges, texts, autotexts = ax.pie(sizes, labels=labels, colors=colors, + autopct='%1.1f%%', startangle=90) + ax.set_title('Peak Partition Distribution') + + plt.tight_layout() + plt.savefig(output_pdf, format='pdf', dpi=150) + if output_png: + plt.savefig(output_png, format='png', dpi=150) + plt.close() + + +def plot_chrom_distribution(query_bed, chrom_sizes, output_pdf, output_png=None): + """Plot chromosome distribution. + + Args: + query_bed: Path to query regions BED file + chrom_sizes: Path to chromosome sizes file + output_pdf: Output PDF path + output_png: Output PNG path (optional) + """ + try: + from gtars.models import RegionSet + except ImportError: + raise ImportError("gtars not available. Install with: pip install gtars") + + # Load data + rs = RegionSet(query_bed) + stats = rs.chromosome_statistics() + + # Sort chromosomes naturally (chr1, chr2, ... chr10, chr11, ... chrX, chrY) + def chrom_sort_key(x): + x = x.replace('chr', '') + if x.isdigit(): + return (0, int(x)) + else: + return (1, x) + + chroms = sorted(stats.keys(), key=chrom_sort_key) + counts = [stats[c].number_of_regions for c in chroms] + + # Plot + fig, ax = plt.subplots(figsize=(12, 6)) + ax.bar(range(len(chroms)), counts, color='#1f77b4') + ax.set_xticks(range(len(chroms))) + ax.set_xticklabels(chroms, rotation=45, ha='right') + ax.set_xlabel('Chromosome') + ax.set_ylabel('Region Count') + ax.set_title('Chromosome Distribution') + + plt.tight_layout() + plt.savefig(output_pdf, format='pdf', dpi=150) + if output_png: + plt.savefig(output_png, format='png', dpi=150) + plt.close() + + +if __name__ == '__main__': + import argparse + + parser = argparse.ArgumentParser(description='PEPATAC QC with gtars backend') + subparsers = parser.add_subparsers(dest='command') + + # TSS subcommand + tss_parser = subparsers.add_parser('tss', help='Plot TSS enrichment') + tss_parser.add_argument('-i', '--input', required=True, help='TSS enrichment file') + tss_parser.add_argument('-o', '--output', required=True, help='Output PDF') + + # Fragment subcommand + frag_parser = subparsers.add_parser('frag', help='Plot fragment distribution') + frag_parser.add_argument('-l', '--lengths', required=True, help='Fragment lengths file') + frag_parser.add_argument('-c', '--counts', required=True, help='Fragment counts file') + frag_parser.add_argument('-p', '--pdf', required=True, help='Output PDF') + frag_parser.add_argument('-t', '--txt', required=True, help='Output stats file') + + args = parser.parse_args() + + if args.command == 'tss': + png_path = args.output.replace('.pdf', '.png') + plot_tss_enrichment(args.input, args.output, png_path) + elif args.command == 'frag': + png_path = args.pdf.replace('.pdf', '.png') + plot_fragment_distribution(args.lengths, args.counts, args.pdf, args.txt, png_path) From 0391eb6ff69f2b2005c26ea7f421f9a2f256c3e2 Mon Sep 17 00:00:00 2001 From: nsheff Date: Mon, 18 May 2026 11:59:27 -0400 Subject: [PATCH 8/9] Add Python summarizer (project-level) replacing R - tools/pepatac_summarizer/: Python package with CLI, consensus peak calling via gtars, peak counts, and summary plots - pipelines/pepatac_collator.py: --summarizer python|R dispatch, defaults to python - Remove obsolete PEPATACr R tests - tests/test_summarizer.py: unit tests - tests/test_summarizer_integration.py: integration tests --- PEPATACr/tests/testthat.R | 3 - PEPATACr/tests/testthat/helper-fixtures.R | 79 -------- PEPATACr/tests/testthat/test-summarizer.R | 24 --- PEPATACr/tests/testthat/test-utilities.R | 38 ---- PEPATACr/tests/testthat/test-yamlToDT.R | 78 ------- pipelines/pepatac_collator.py | 15 +- tests/test_summarizer.py | 126 ++++++++++++ tests/test_summarizer_integration.py | 159 +++++++++++++++ tools/pepatac_summarizer/__init__.py | 3 + tools/pepatac_summarizer/__main__.py | 7 + tools/pepatac_summarizer/assets.py | 46 +++++ tools/pepatac_summarizer/cli.py | 116 +++++++++++ tools/pepatac_summarizer/consensus.py | 181 +++++++++++++++++ tools/pepatac_summarizer/counts.py | 98 +++++++++ tools/pepatac_summarizer/plots/__init__.py | 14 ++ tools/pepatac_summarizer/plots/alignment.py | 202 +++++++++++++++++++ tools/pepatac_summarizer/plots/complexity.py | 111 ++++++++++ tools/pepatac_summarizer/plots/library.py | 72 +++++++ tools/pepatac_summarizer/plots/theme.py | 43 ++++ tools/pepatac_summarizer/plots/tss.py | 85 ++++++++ tools/pepatac_summarizer/utils.py | 78 +++++++ 21 files changed, 1353 insertions(+), 225 deletions(-) delete mode 100644 PEPATACr/tests/testthat.R delete mode 100644 PEPATACr/tests/testthat/helper-fixtures.R delete mode 100644 PEPATACr/tests/testthat/test-summarizer.R delete mode 100644 PEPATACr/tests/testthat/test-utilities.R delete mode 100644 PEPATACr/tests/testthat/test-yamlToDT.R create mode 100644 tests/test_summarizer.py create mode 100644 tests/test_summarizer_integration.py create mode 100644 tools/pepatac_summarizer/__init__.py create mode 100644 tools/pepatac_summarizer/__main__.py create mode 100644 tools/pepatac_summarizer/assets.py create mode 100644 tools/pepatac_summarizer/cli.py create mode 100644 tools/pepatac_summarizer/consensus.py create mode 100644 tools/pepatac_summarizer/counts.py create mode 100644 tools/pepatac_summarizer/plots/__init__.py create mode 100644 tools/pepatac_summarizer/plots/alignment.py create mode 100644 tools/pepatac_summarizer/plots/complexity.py create mode 100644 tools/pepatac_summarizer/plots/library.py create mode 100644 tools/pepatac_summarizer/plots/theme.py create mode 100644 tools/pepatac_summarizer/plots/tss.py create mode 100644 tools/pepatac_summarizer/utils.py diff --git a/PEPATACr/tests/testthat.R b/PEPATACr/tests/testthat.R deleted file mode 100644 index 88e036e..0000000 --- a/PEPATACr/tests/testthat.R +++ /dev/null @@ -1,3 +0,0 @@ -library(testthat) -library(PEPATACr) -test_check("PEPATACr") diff --git a/PEPATACr/tests/testthat/helper-fixtures.R b/PEPATACr/tests/testthat/helper-fixtures.R deleted file mode 100644 index 7c3add0..0000000 --- a/PEPATACr/tests/testthat/helper-fixtures.R +++ /dev/null @@ -1,79 +0,0 @@ -# Make internal function available to all tests -yamlToDT <- PEPATACr:::yamlToDT - -# Fixture: scalars only (no object-type fields) -- the happy path -make_scalar_only_yaml <- function() { - list(PEPATAC = list(sample = list( - sample1 = list( - Raw_reads = 25000, Aligned_reads = 500, Alignment_rate = 2.0, - Peak_count = 19, FRiP = 0.108, Genome = "hg38", - Time = "0:00:27", Success = "02-13-18:32:54" - ) - ))) -} - -# Fixture: mixed scalars + ONE object field (pre-bug state, already in remove_cols) -make_one_object_yaml <- function() { - yaml <- make_scalar_only_yaml() - yaml$PEPATAC$sample$sample1$`Fragment distribution` <- list( - path = "/path/to/frag.pdf", - thumbnail_path = "/path/to/frag.png", - title = "Fragment distribution" - ) - yaml -} - -# Fixture: mixed scalars + NEW object fields (the bug trigger) -# Includes Motif analysis, Library complexity, TSS enrichment -- -# the three fields added in commit 60b8635 that are NOT in remove_cols -make_bug_trigger_yaml <- function() { - yaml <- make_one_object_yaml() - s <- yaml$PEPATAC$sample$sample1 - s$`Library complexity` <- list( - path = "/path/to/lib.pdf", thumbnail_path = "/path/to/lib.png", - title = "Library complexity" - ) - s$`TSS enrichment` <- list( - path = "/path/to/tss.pdf", thumbnail_path = "/path/to/tss.png", - title = "TSS enrichment" - ) - s$`Motif analysis` <- list( - path = "/path/to/motif.html", title = "Motif analysis", - thumbnail_path = "/path/to/motif.png" - ) - s$meta <- list( - pipestat_modified_time = "2026-02-13 18:32:54", - pipestat_created_time = "2026-02-13 18:32:27" - ) - # Also add all the other object types that ARE in remove_cols - s$`Peak chromosome distribution` <- list( - path = "/p.pdf", thumbnail_path = "/p.png", title = "Peak chromosome distribution" - ) - s$`Peak partition distribution` <- list( - path = "/p.pdf", thumbnail_path = "/p.png", title = "Peak partition distribution" - ) - s$cFRiF <- list(path = "/p.pdf", thumbnail_path = "/p.png", title = "cFRiF") - s$FRiF <- list(path = "/p.pdf", thumbnail_path = "/p.png", title = "FRiF") - s$`FastQC report r1` <- list(path = "/p.html", thumbnail_path = "/p.png", title = "FastQC r1") - s$`FastQC report r2` <- list(path = "/p.html", thumbnail_path = "/p.png", title = "FastQC r2") - s$`TSS distance distribution` <- list( - path = "/p.pdf", thumbnail_path = "/p.png", title = "TSS distance distribution" - ) - yaml$PEPATAC$sample$sample1 <- s - yaml -} - -# Fixture: multiple samples with varying completeness -make_multi_sample_yaml <- function() { - yaml <- make_bug_trigger_yaml() - # sample2: complete but different values - yaml$PEPATAC$sample$sample2 <- yaml$PEPATAC$sample$sample1 - yaml$PEPATAC$sample$sample2$Raw_reads <- 50000 - yaml$PEPATAC$sample$sample2$Aligned_reads <- 1200 - # sample3: incomplete -- missing some fields (simulates partial pipeline run) - yaml$PEPATAC$sample$sample3 <- list( - Raw_reads = 30000, Aligned_reads = 800, - Genome = "hg38", meta = list(pipestat_modified_time = "2026-02-13") - ) - yaml -} diff --git a/PEPATACr/tests/testthat/test-summarizer.R b/PEPATACr/tests/testthat/test-summarizer.R deleted file mode 100644 index 4972b56..0000000 --- a/PEPATACr/tests/testthat/test-summarizer.R +++ /dev/null @@ -1,24 +0,0 @@ -test_that("rbindlist of yamlToDT across samples produces correct table", { - yaml <- make_multi_sample_yaml() - sample_names <- names(yaml$PEPATAC$sample) - stats <- data.table::rbindlist( - lapply(sample_names, FUN = yamlToDT, yaml_file = yaml), - fill = TRUE - ) - expect_s3_class(stats, "data.table") - expect_equal(nrow(stats), 3) - expect_true(all(c("sample1", "sample2", "sample3") %in% stats$sample_name)) - # sample3 is missing Peak_count -- should be NA (filled by rbindlist) - expect_true(is.na(stats[sample_name == "sample3"]$Peak_count)) -}) - -test_that("rbindlist with no recycling warnings across heterogeneous samples", { - yaml <- make_multi_sample_yaml() - sample_names <- names(yaml$PEPATAC$sample) - expect_no_warning( - stats <- data.table::rbindlist( - lapply(sample_names, FUN = yamlToDT, yaml_file = yaml), - fill = TRUE - ) - ) -}) diff --git a/PEPATACr/tests/testthat/test-utilities.R b/PEPATACr/tests/testthat/test-utilities.R deleted file mode 100644 index bd3c69c..0000000 --- a/PEPATACr/tests/testthat/test-utilities.R +++ /dev/null @@ -1,38 +0,0 @@ -test_that("sampleName extracts name from path", { - result <- sampleName("/path/to/sample_R1.fastq.gz") - expect_type(result, "character") - expect_true(nchar(result) > 0) -}) - -test_that("splitDataTable splits by column", { - dt <- data.table::data.table( - name = c("a", "b", "c"), group = c("x", "x", "y"), val = 1:3 - ) - result <- splitDataTable(dt, "group") - expect_type(result, "list") - expect_equal(length(result), 2) -}) - -test_that("roundUpNice rounds to nice numbers", { - expect_equal(roundUpNice(13), 20) - expect_equal(roundUpNice(0.7), 0.7) -}) - -test_that("is.empty detects empty data.frames", { - is.empty <- PEPATACr:::is.empty - expect_true(is.empty(NULL)) - expect_true(is.empty(data.frame())) - expect_false(is.empty(data.frame(x = 1))) -}) - -test_that("getAbbr returns correct abbreviations", { - getAbbr <- PEPATACr:::getAbbr - expect_equal(getAbbr(1500), "K") - expect_equal(getAbbr(2000000), "M") -}) - -test_that("getFactor returns correct divisors", { - getFactor <- PEPATACr:::getFactor - expect_equal(getFactor(1500), 1000) - expect_equal(getFactor(2000000), 1000000) -}) diff --git a/PEPATACr/tests/testthat/test-yamlToDT.R b/PEPATACr/tests/testthat/test-yamlToDT.R deleted file mode 100644 index 9600510..0000000 --- a/PEPATACr/tests/testthat/test-yamlToDT.R +++ /dev/null @@ -1,78 +0,0 @@ -test_that("yamlToDT handles scalar-only YAML without warnings", { - yaml <- make_scalar_only_yaml() - expect_no_warning( - dt <- yamlToDT("sample1", yaml) - ) - expect_s3_class(dt, "data.table") - expect_equal(nrow(dt), 1) - expect_equal(dt$sample_name, "sample1") - expect_equal(dt$Raw_reads, 25000) -}) - -test_that("yamlToDT produces no recycling warning with object-type fields", { - # THIS IS THE BUG TEST -- currently fails with: - # "Item N has 3 rows but longest item has 4; recycled with remainder." - yaml <- make_bug_trigger_yaml() - expect_no_warning( - dt <- yamlToDT("sample1", yaml) - ) - expect_s3_class(dt, "data.table") - expect_equal(nrow(dt), 1) -}) - -test_that("yamlToDT excludes ALL object-type columns", { - yaml <- make_bug_trigger_yaml() - dt <- suppressWarnings(yamlToDT("sample1", yaml)) - # None of these object-type fields should appear as columns - object_fields <- c("Fragment distribution", "Library complexity", - "TSS enrichment", "Motif analysis", - "Peak chromosome distribution", "Peak partition distribution", - "cFRiF", "FRiF", "FastQC report r1", "FastQC report r2", - "TSS distance distribution", "meta") - for (field in object_fields) { - expect_false(field %in% names(dt), - info = paste0("Object field '", field, "' should not be in output")) - } -}) - -test_that("yamlToDT preserves all scalar fields", { - yaml <- make_bug_trigger_yaml() - dt <- suppressWarnings(yamlToDT("sample1", yaml)) - expect_true("Raw_reads" %in% names(dt)) - expect_true("Aligned_reads" %in% names(dt)) - expect_true("Alignment_rate" %in% names(dt)) - expect_true("Peak_count" %in% names(dt)) - expect_true("FRiP" %in% names(dt)) - expect_true("Genome" %in% names(dt)) -}) - -test_that("yamlToDT returns exactly one row per sample", { - yaml <- make_bug_trigger_yaml() - dt <- suppressWarnings(yamlToDT("sample1", yaml)) - # The old code used unique() to collapse recycled rows. - # The fix should produce exactly 1 row without needing unique(). - expect_equal(nrow(dt), 1) -}) - -test_that("yamlToDT returns NULL for missing sample", { - yaml <- make_scalar_only_yaml() - result <- suppressWarnings(tryCatch( - yamlToDT("nonexistent_sample", yaml), - error = function(e) NULL - )) - expect_null(result) -}) - -test_that("yamlToDT handles sample with only meta (no scalar results)", { - yaml <- list(PEPATAC = list(sample = list( - empty_sample = list( - meta = list(pipestat_modified_time = "2026-02-13") - ) - ))) - result <- suppressWarnings(tryCatch( - yamlToDT("empty_sample", yaml), - error = function(e) NULL - )) - # Should return NULL or empty since no scalar data exists after removing meta - expect_true(is.null(result) || nrow(result) == 0) -}) diff --git a/pipelines/pepatac_collator.py b/pipelines/pepatac_collator.py index dff3fbf..8a3f8b0 100755 --- a/pipelines/pepatac_collator.py +++ b/pipelines/pepatac_collator.py @@ -68,6 +68,9 @@ def parse_arguments(): help="Path to a reference peak set (narrowPeak/BED) " "to use for the project count table instead of " "the computed consensus peaks.") + parser.add_argument("--summarizer", default="python", + choices=["python", "R"], + help="Summarizer implementation to use (default: python)") args = parser.parse_args() return args @@ -106,9 +109,15 @@ def main(): yaml.dump(yaml_dict, file) print(f"Summary (n={num_samples}: {project_stats_file})") - cmd = (f"Rscript {tool_path('PEPATAC_summarizer.R')} " - f"{args.config_file} {args.output_parent} " - f"{args.results} {args.cutoff} {args.min_score} {args.min_olap}") + if args.summarizer == "python": + cmd = (f"python -m pepatac_summarizer " + f"{args.config_file} {args.output_parent} " + f"{args.results} --cutoff {args.cutoff} " + f"--min-score {args.min_score} --min-olap {args.min_olap}") + else: + cmd = (f"Rscript {tool_path('PEPATAC_summarizer.R')} " + f"{args.config_file} {args.output_parent} " + f"{args.results} {args.cutoff} {args.min_score} {args.min_olap}") if args.new_start: cmd += " --new-start" if args.skip_consensus: diff --git a/tests/test_summarizer.py b/tests/test_summarizer.py new file mode 100644 index 0000000..427475c --- /dev/null +++ b/tests/test_summarizer.py @@ -0,0 +1,126 @@ +"""Tests for the Python summarizer module.""" + +import pytest +from pathlib import Path +import sys + +sys.path.insert(0, str(Path(__file__).parent.parent / "tools")) + +from pepatac_summarizer.assets import create_assets_summary +from pepatac_summarizer.consensus import collapse_peaks, parse_narrowpeak_line + + +@pytest.fixture +def sample_results(tmp_path): + """Create a minimal sample results structure.""" + sample_dir = tmp_path / "sample1" + aligned_dir = sample_dir / "aligned_hg38" + peak_dir = sample_dir / "peak_calling_hg38" + qc_dir = sample_dir / "QC_hg38" + + aligned_dir.mkdir(parents=True) + peak_dir.mkdir(parents=True) + qc_dir.mkdir(parents=True) + + (aligned_dir / "sample1_sort_dedup.bam").touch() + (aligned_dir / "sample1_sort_dedup.bam.bai").touch() + (aligned_dir / "hg38.chrom.sizes").write_text("chr1\t248956422\nchr2\t242193529\n") + + peak_content = "chr1\t1000\t2000\tpeak1\t100\t.\t50.0\t10.0\t5.0\t500\n" + (peak_dir / "sample1_peaks_normalized.narrowPeak").write_text(peak_content) + + preseq_content = "TOTAL_READS\tEXPECTED_DISTINCT\n1000000\t800000\n2000000\t1500000\n" + (qc_dir / "sample1_preseq_yield.txt").write_text(preseq_content) + + return tmp_path + + +def test_create_assets_summary(sample_results): + """Test assets summary generation.""" + assets = create_assets_summary(["sample1"], str(sample_results)) + + assert not assets.empty + assert "sample_name" in assets.columns + assert "asset" in assets.columns + assert "path" in assets.columns + assert "sample1" in assets["sample_name"].values + + +def test_parse_narrowpeak_line(): + """Test narrowPeak line parsing.""" + line = "chr1\t1000\t2000\tpeak1\t100\t.\t50.0\t10.0\t5.0\t500" + chrom, start, end, name, score, rest = parse_narrowpeak_line(line) + + assert chrom == "chr1" + assert start == 1000 + assert end == 2000 + assert name == "peak1" + assert score == 100.0 + + +def test_collapse_peaks_single_file(sample_results): + """Test collapse with single file returns empty (need 2+ samples).""" + peak_file = sample_results / "sample1" / "peak_calling_hg38" / "sample1_peaks_normalized.narrowPeak" + chrom_sizes = {"chr1": 248956422} + + result = collapse_peaks( + [str(peak_file)], ["sample1"], chrom_sizes, + min_samples=2, min_score=5.0, min_olap=1 + ) + assert len(result) == 0 + + +def test_collapse_peaks_two_overlapping(tmp_path): + """Test consensus with two overlapping peak files.""" + peak1 = tmp_path / "sample1.narrowPeak" + peak2 = tmp_path / "sample2.narrowPeak" + + peak1.write_text("chr1\t1000\t2000\tpeak1\t100\t.\t50.0\t10.0\t5.0\t500\n") + peak2.write_text("chr1\t1050\t2050\tpeak2\t80\t.\t40.0\t8.0\t4.0\t500\n") + + chrom_sizes = {"chr1": 248956422} + + result = collapse_peaks( + [str(peak1), str(peak2)], ["sample1", "sample2"], chrom_sizes, + min_samples=2, min_score=5.0, min_olap=1 + ) + + assert len(result) == 1 + # Best score should be 100 + assert "100" in result[0] + + +def test_collapse_peaks_non_overlapping(tmp_path): + """Test that non-overlapping peaks don't merge.""" + peak1 = tmp_path / "sample1.narrowPeak" + peak2 = tmp_path / "sample2.narrowPeak" + + peak1.write_text("chr1\t1000\t2000\tpeak1\t100\t.\t50.0\t10.0\t5.0\t500\n") + peak2.write_text("chr1\t5000\t6000\tpeak2\t80\t.\t40.0\t8.0\t4.0\t500\n") + + chrom_sizes = {"chr1": 248956422} + + result = collapse_peaks( + [str(peak1), str(peak2)], ["sample1", "sample2"], chrom_sizes, + min_samples=2, min_score=5.0, min_olap=1 + ) + + assert len(result) == 0 + + +def test_collapse_peaks_min_score_filter(tmp_path): + """Test min_score filtering.""" + peak1 = tmp_path / "sample1.narrowPeak" + peak2 = tmp_path / "sample2.narrowPeak" + + peak1.write_text("chr1\t1000\t2000\tpeak1\t3\t.\t50.0\t10.0\t5.0\t500\n") + peak2.write_text("chr1\t1050\t2050\tpeak2\t3\t.\t40.0\t8.0\t4.0\t500\n") + + chrom_sizes = {"chr1": 248956422} + + result = collapse_peaks( + [str(peak1), str(peak2)], ["sample1", "sample2"], chrom_sizes, + min_samples=2, min_score=5.0, min_olap=1 + ) + + assert len(result) == 0 diff --git a/tests/test_summarizer_integration.py b/tests/test_summarizer_integration.py new file mode 100644 index 0000000..d7889b4 --- /dev/null +++ b/tests/test_summarizer_integration.py @@ -0,0 +1,159 @@ +"""Integration tests for the Python summarizer module with synthetic data.""" + +import pytest +import tempfile +from pathlib import Path +import sys +import yaml + +sys.path.insert(0, str(Path(__file__).parent.parent / "tools")) + +from pepatac_summarizer.assets import create_assets_summary +from pepatac_summarizer.consensus import calculate_consensus_peaks, collapse_peaks +from pepatac_summarizer.plots import plot_complexity_curves + + +@pytest.fixture +def multi_sample_project(tmp_path): + """Create a minimal multi-sample project with synthetic PEPATAC outputs.""" + results_dir = tmp_path / "results" + + for i, sample in enumerate(["sample1", "sample2", "sample3"]): + sample_dir = results_dir / sample + aligned_dir = sample_dir / "aligned_hg38" + peak_dir = sample_dir / "peak_calling_hg38" + qc_dir = sample_dir / "QC_hg38" + + aligned_dir.mkdir(parents=True) + peak_dir.mkdir(parents=True) + qc_dir.mkdir(parents=True) + + (aligned_dir / f"{sample}_sort_dedup.bam").touch() + (aligned_dir / f"{sample}_sort_dedup.bam.bai").touch() + (aligned_dir / "hg38.chrom.sizes").write_text( + "chr1\t248956422\nchr2\t242193529\nchr22\t50818468\n" + ) + + offset = i * 100 + peaks = [ + f"chr1\t{1000+offset}\t{2000+offset}\tpeak1\t{100-i*10}\t.\t50.0\t10.0\t5.0\t500", + f"chr1\t{5000+offset}\t{6000+offset}\tpeak2\t{80-i*5}\t.\t40.0\t8.0\t4.0\t500", + f"chr22\t{10000+offset}\t{11000+offset}\tpeak3\t{90-i*5}\t.\t45.0\t9.0\t4.5\t500", + ] + (peak_dir / f"{sample}_peaks_normalized.narrowPeak").write_text("\n".join(peaks) + "\n") + + preseq = "TOTAL_READS\tEXPECTED_DISTINCT\n" + for j in range(1, 11): + preseq += f"{j*1000000}\t{int(j*800000*(1-j*0.02))}\n" + (qc_dir / f"{sample}_preseq_yield.txt").write_text(preseq) + + counts = f"V1\tV2\tV3\n{sample}\t{5000000+i*100000}\t{4000000+i*80000}\n" + (qc_dir / f"{sample}_preseq_counts.txt").write_text(counts) + + (sample_dir / "stats.yaml").write_text(yaml.dump({ + "PEPATAC": {"sample": {sample: {"aligned_reads": 5000000}}} + })) + + config_file = tmp_path / "project_config.yaml" + sample_table_file = tmp_path / "samples.csv" + + sample_table_file.write_text( + "sample_name,genome\nsample1,hg38\nsample2,hg38\nsample3,hg38\n" + ) + + config = { + "name": "test_project", + "pep_version": "2.1.0", + "sample_table": str(sample_table_file), + } + config_file.write_text(yaml.dump(config)) + + return { + "config": str(config_file), + "output": str(tmp_path), + "results": str(results_dir), + "tmp_path": tmp_path, + } + + +def test_full_summarizer_pipeline(multi_sample_project): + """Test complete summarizer pipeline with synthetic multi-sample project.""" + project = multi_sample_project + summary_dir = Path(project["output"]) / "summary" + + sample_table = { + "sample_name": ["sample1", "sample2", "sample3"], + "genome": ["hg38", "hg38", "hg38"] + } + + # Test assets + assets = create_assets_summary( + ["sample1", "sample2", "sample3"], + project["results"] + ) + assert not assets.empty + assert len(assets) >= 9 + + # Test consensus peaks + consensus_files = calculate_consensus_peaks( + sample_table, + str(summary_dir), + project["results"], + "test_project", + min_samples=2, + min_score=5.0, + min_olap=1 + ) + + assert "hg38" in consensus_files + with open(consensus_files["hg38"]) as f: + lines = [l for l in f if l.strip()] + assert len(lines) >= 1 + + # Test plots + plot_file = plot_complexity_curves( + ["sample1", "sample2", "sample3"], + ["hg38", "hg38", "hg38"], + project["results"], + str(summary_dir), + "test_project" + ) + + assert plot_file is not None + assert Path(plot_file).exists() + + +def test_consensus_with_chr22_peaks(tmp_path): + """Test consensus specifically with chr22 peaks (accbase-style).""" + peaks1 = tmp_path / "s1.narrowPeak" + peaks2 = tmp_path / "s2.narrowPeak" + peaks3 = tmp_path / "s3.narrowPeak" + + peaks1.write_text( + "chr22\t16050000\t16051000\tp1\t100\t.\t50\t10\t5\t500\n" + "chr22\t20000000\t20001000\tp2\t80\t.\t40\t8\t4\t500\n" + ) + peaks2.write_text( + "chr22\t16050100\t16051100\tp1\t90\t.\t45\t9\t4.5\t500\n" + "chr22\t20000050\t20001050\tp2\t85\t.\t42\t8.5\t4.2\t500\n" + ) + peaks3.write_text( + "chr22\t16050200\t16051200\tp1\t95\t.\t48\t9.5\t4.8\t500\n" + "chr22\t30000000\t30001000\tp3\t70\t.\t35\t7\t3.5\t500\n" + ) + + chrom_sizes = {"chr22": 50818468} + + result = collapse_peaks( + [str(peaks1), str(peaks2), str(peaks3)], + ["s1", "s2", "s3"], + chrom_sizes, + min_samples=2, + min_score=5.0, + min_olap=1 + ) + + assert len(result) >= 2 + # Check we got chr22 peaks + chr22_count = sum(1 for line in result if line.startswith("chr22")) + assert chr22_count >= 2 diff --git a/tools/pepatac_summarizer/__init__.py b/tools/pepatac_summarizer/__init__.py new file mode 100644 index 0000000..9079b53 --- /dev/null +++ b/tools/pepatac_summarizer/__init__.py @@ -0,0 +1,3 @@ +"""PEPATAC Project Summarizer - Python implementation.""" + +__version__ = "0.1.0" diff --git a/tools/pepatac_summarizer/__main__.py b/tools/pepatac_summarizer/__main__.py new file mode 100644 index 0000000..01d2064 --- /dev/null +++ b/tools/pepatac_summarizer/__main__.py @@ -0,0 +1,7 @@ +"""Allow running as: python -m pepatac_summarizer""" + +from .cli import main + +if __name__ == "__main__": + import sys + sys.exit(main()) diff --git a/tools/pepatac_summarizer/assets.py b/tools/pepatac_summarizer/assets.py new file mode 100644 index 0000000..73c8d25 --- /dev/null +++ b/tools/pepatac_summarizer/assets.py @@ -0,0 +1,46 @@ +"""Assets summary generation.""" + +from pathlib import Path +import pandas as pd + +ASSET_PATTERNS = { + "aligned_bam": "aligned_*/{}*_sort_dedup.bam", + "aligned_bam_index": "aligned_*/{}*_sort_dedup.bam.bai", + "peaks_bed": "peak_calling_*/{}*_peaks.narrowPeak", + "peaks_normalized": "peak_calling_*/{}*_peaks_normalized.narrowPeak", + "bigwig": "signal_*/{}*_smooth.bw", + "chrom_sizes": "aligned_*/*.chrom.sizes", + "stats": "stats.yaml", + "preseq_yield": "QC_*/{}*_preseq_yield.txt", + "preseq_counts": "QC_*/{}*_preseq_counts.txt", +} + + +def create_assets_summary(sample_names: list, results_subdir: str) -> pd.DataFrame: + """Generate project assets summary. + + Args: + sample_names: List of sample names + results_subdir: Path to results subdirectory + + Returns: + DataFrame with columns: sample_name, asset, path + """ + assets = [] + results_path = Path(results_subdir) + + for sample in sample_names: + sample_dir = results_path / sample + if not sample_dir.exists(): + continue + + for asset_type, pattern in ASSET_PATTERNS.items(): + search_pattern = pattern.format(sample) + for path in sample_dir.glob(search_pattern): + assets.append({ + "sample_name": sample, + "asset": asset_type, + "path": str(path) + }) + + return pd.DataFrame(assets) diff --git a/tools/pepatac_summarizer/cli.py b/tools/pepatac_summarizer/cli.py new file mode 100644 index 0000000..8829252 --- /dev/null +++ b/tools/pepatac_summarizer/cli.py @@ -0,0 +1,116 @@ +#!/usr/bin/env python3 +"""PEPATAC Project Summarizer CLI - Python replacement for PEPATAC_summarizer.R""" + +import argparse +from pathlib import Path +import pandas as pd +import peppy + +from .assets import create_assets_summary +from .consensus import calculate_consensus_peaks +from .counts import calculate_peak_counts +from .utils import load_stats_summary +from .plots import ( + plot_aligned_raw, + plot_aligned_pct, + plot_tss_scores, + plot_lib_sizes, + plot_complexity_curves, +) + + +def main(): + parser = argparse.ArgumentParser(description="PEPATAC Project Summarizer") + parser.add_argument("config", help="PEP project config YAML") + parser.add_argument("output", help="Project output directory") + parser.add_argument("results", help="Results subdirectory path") + parser.add_argument("-N", "--new-start", action="store_true", + help="Start over, run every command") + parser.add_argument("-P", "--skip-consensus", action="store_true", + help="Do not calculate consensus peaks") + parser.add_argument("-T", "--skip-table", action="store_true", + help="Do not calculate peak counts table") + parser.add_argument("-m", "--cutoff", type=int, default=2, + help="Min samples for consensus peaks") + parser.add_argument("-s", "--min-score", type=float, default=5, + help="Min peak score") + parser.add_argument("-l", "--min-olap", type=int, default=1, + help="Min overlap bases") + parser.add_argument("-F", "--frip-ref-peaks", + help="Reference peak set for counts table") + parser.add_argument("-V", "--poverlap", action="store_true", + help="Calculate percentage overlap in counts table") + parser.add_argument("-Z", "--normalized", action="store_true", + help="Use normalized read counts") + + args = parser.parse_args() + + project = peppy.Project(args.config) + project_name = project.name + + sample_table = pd.DataFrame({ + "sample_name": project.sample_table["sample_name"], + "genome": project.sample_table["genome"] + }) + + results_subdir = Path(args.results) + if not results_subdir.exists(): + print(f"Error: Results subdirectory does not exist: {results_subdir}") + return 1 + + output_dir = Path(args.output) + summary_dir = output_dir / "summary" + summary_dir.mkdir(parents=True, exist_ok=True) + + assets = create_assets_summary( + list(project.sample_table["sample_name"]), + str(results_subdir) + ) + + if assets.empty: + print("No assets found - exiting") + return 1 + + assets_file = output_dir / f"{project_name}_assets_summary.tsv" + assets.to_csv(assets_file, sep="\t", index=False, header=False) + print(f"Summary (n={len(assets['sample_name'].unique())}): {assets_file}") + + stats_file = output_dir / f"{project_name}_stats_summary.yaml" + stats = load_stats_summary(stats_file) + + if stats is not None and not stats.empty: + print("Creating summary plots...") + plot_aligned_raw(stats, str(summary_dir), project_name) + plot_aligned_pct(stats, str(summary_dir), project_name) + plot_tss_scores(stats, str(summary_dir), project_name) + plot_lib_sizes(stats, str(summary_dir), project_name) + else: + print("Warning: No stats summary available, skipping summary plots") + + sample_names = list(project.sample_table["sample_name"]) + genomes = list(project.sample_table["genome"]) + plot_complexity_curves( + sample_names, genomes, + str(results_subdir), str(summary_dir), project_name + ) + + consensus_files = {} + if not args.skip_consensus: + consensus_files = calculate_consensus_peaks( + sample_table, str(summary_dir), str(results_subdir), project_name, + min_samples=args.cutoff, min_score=args.min_score, min_olap=args.min_olap + ) + + if not args.skip_table and consensus_files: + calculate_peak_counts( + sample_table, str(summary_dir), str(results_subdir), project_name, + consensus_files, normalized=args.normalized, poverlap=args.poverlap + ) + + print("Successfully produced project summary.") + return 0 + + +if __name__ == "__main__": + import sys + sys.exit(main()) diff --git a/tools/pepatac_summarizer/consensus.py b/tools/pepatac_summarizer/consensus.py new file mode 100644 index 0000000..d0cf8d7 --- /dev/null +++ b/tools/pepatac_summarizer/consensus.py @@ -0,0 +1,181 @@ +"""Consensus peak calculation using gtars.""" + +from pathlib import Path +from gtars.models import RegionSet, Region + + +def parse_narrowpeak_line(line: str) -> tuple[str, int, int, str, float, str]: + """Parse a narrowPeak line into components.""" + parts = line.strip().split("\t") + chrom = parts[0] + start = int(parts[1]) + end = int(parts[2]) + name = parts[3] if len(parts) > 3 else "." + score = float(parts[4]) if len(parts) > 4 else 0.0 + rest = "\t".join(parts[3:]) if len(parts) > 3 else "" + return chrom, start, end, name, score, rest + + +def collapse_peaks( + peak_files: list[str], + sample_names: list[str], + chrom_sizes: dict[str, int], + min_samples: int = 2, + min_score: float = 5.0, + min_olap: int = 1 +) -> list[str]: + """Collapse overlapping peaks using gtars AIList. + + Args: + peak_files: List of narrowPeak file paths + sample_names: Corresponding sample names + chrom_sizes: Dict mapping chromosome -> size + min_samples: Minimum samples a peak must appear in + min_score: Minimum score to keep a peak + min_olap: Minimum overlap in bp + + Returns: + List of narrowPeak lines for consensus peaks + """ + # Parse all peaks with metadata + all_regions = [] + peak_data = [] # (sample_id, score, full_line) + + for sample_id, (pf, sample) in enumerate(zip(peak_files, sample_names)): + try: + with open(pf) as f: + for line in f: + if not line.strip(): + continue + chrom, start, end, name, score, rest = parse_narrowpeak_line(line) + all_regions.append(Region(chrom, start, end, "")) + peak_data.append((sample_id, score, line.strip())) + except (FileNotFoundError, IOError): + continue + + if len(all_regions) < 2: + return [] + + # Build RegionSets + combined_rs = RegionSet.from_regions(all_regions) + merged_rs = combined_rs.reduce() + + # Use gtars find_overlaps + overlap_indices = merged_rs.find_overlaps(combined_rs) + + consensus_lines = [] + for indices in overlap_indices: + if not indices: + continue + + # Find best peak and count samples + samples_seen = set() + best_score = -1.0 + best_line = "" + + for idx in indices: + sample_id, score, line = peak_data[idx] + samples_seen.add(sample_id) + if score > best_score: + best_score = score + best_line = line + + if len(samples_seen) >= min_samples and best_score >= min_score: + consensus_lines.append(best_line) + + return consensus_lines + + +def calculate_consensus_peaks( + sample_table, # Can be dict or DataFrame + summary_dir: str, + results_subdir: str, + project_name: str, + min_samples: int = 2, + min_score: float = 5.0, + min_olap: int = 1 +) -> dict[str, str]: + """Calculate consensus peaks per genome using gtars. + + Args: + sample_table: Dict or DataFrame with sample_name and genome + summary_dir: Output directory for consensus files + results_subdir: Path to sample results + project_name: Project name for output files + min_samples: Minimum samples for reproducibility + min_score: Minimum peak score + min_olap: Minimum overlap bases + + Returns: + Dict mapping genome -> consensus peak file path + """ + summary_path = Path(summary_dir) + summary_path.mkdir(exist_ok=True) + results_path = Path(results_subdir) + + # Handle both dict and DataFrame input + if hasattr(sample_table, 'iterrows'): + samples = [(row["sample_name"], row["genome"]) for _, row in sample_table.iterrows()] + else: + samples = [(s, g) for s, g in zip(sample_table["sample_name"], sample_table["genome"])] + + # Group by genome + genome_samples: dict[str, list[str]] = {} + for sample, genome in samples: + genome_samples.setdefault(genome, []).append(sample) + + consensus_files = {} + + for genome, sample_list in genome_samples.items(): + peak_files = [] + sample_names = [] + chrom_sizes = {} + + for sample in sample_list: + peak_file = ( + results_path / sample / f"peak_calling_{genome}" / + f"{sample}_peaks_normalized.narrowPeak" + ) + if peak_file.exists(): + peak_files.append(str(peak_file)) + sample_names.append(sample) + + # Load chrom sizes once + if not chrom_sizes: + chrom_dir = results_path / sample / f"aligned_{genome}" + for cs_file in chrom_dir.glob("*.chrom.sizes"): + with open(cs_file) as f: + for line in f: + parts = line.strip().split("\t") + if len(parts) >= 2: + chrom_sizes[parts[0]] = int(parts[1]) + break + + if len(peak_files) < 2: + print(f"Found only {len(peak_files)} valid peak file(s) for {genome}, skipping consensus.") + continue + + print(f"Calculating {genome} consensus peak set from {len(peak_files)} samples...") + + consensus_lines = collapse_peaks( + peak_files, sample_names, chrom_sizes, + min_samples=min_samples, min_score=min_score, min_olap=min_olap + ) + + if not consensus_lines: + print(f"Warning: No consensus peaks found for {genome}") + continue + + output_file = summary_path / f"{project_name}_{genome}_consensusPeaks.narrowPeak" + with open(output_file, "w") as f: + # Remove duplicates while preserving order + seen = set() + for line in consensus_lines: + if line not in seen: + seen.add(line) + f.write(line + "\n") + + consensus_files[genome] = str(output_file) + print(f"Consensus peak set: {output_file}") + + return consensus_files diff --git a/tools/pepatac_summarizer/counts.py b/tools/pepatac_summarizer/counts.py new file mode 100644 index 0000000..21fc9ee --- /dev/null +++ b/tools/pepatac_summarizer/counts.py @@ -0,0 +1,98 @@ +"""Peak counts table generation using gtars.""" + +from pathlib import Path +import pandas as pd +from gtars.models import RegionSet, Region + + +def calculate_peak_counts( + sample_table: pd.DataFrame, + summary_dir: str, + results_subdir: str, + project_name: str, + consensus_peaks: dict[str, str], + normalized: bool = False, + poverlap: bool = False, +) -> dict[str, str]: + """Generate peak counts table per genome using gtars. + + Args: + sample_table: DataFrame with sample_name and genome columns + summary_dir: Output directory + results_subdir: Path to sample results + project_name: Project name for output files + consensus_peaks: Dict mapping genome -> consensus peak file + normalized: Use normalized read counts (CPM) + poverlap: Calculate percentage overlap (not yet implemented) + + Returns: + Dict mapping genome -> counts table file path + """ + summary_path = Path(summary_dir) + results_path = Path(results_subdir) + count_files = {} + + for genome, consensus_file in consensus_peaks.items(): + if not Path(consensus_file).exists(): + print(f"Consensus file not found: {consensus_file}") + continue + + genome_samples = sample_table[sample_table["genome"] == genome] + + # Load peaks into RegionSet + peaks_df = pd.read_csv( + consensus_file, sep="\t", header=None, + names=["chr", "start", "end", "name", "score", "strand", + "signalValue", "pValue", "qValue", "peak"][:10] + ) + + if peaks_df.empty: + continue + + peak_regions = [ + Region(row["chr"], int(row["start"]), int(row["end"]), "") + for _, row in peaks_df.iterrows() + ] + peaks_rs = RegionSet.from_regions(peak_regions) + + counts_data = { + "chr": peaks_df["chr"].tolist(), + "start": peaks_df["start"].tolist(), + "end": peaks_df["end"].tolist(), + } + + for _, row in genome_samples.iterrows(): + sample = row["sample_name"] + bam_path = ( + results_path / sample / f"aligned_{genome}" / + f"{sample}_sort_dedup.bam" + ) + + if not bam_path.exists(): + print(f"BAM not found: {bam_path}") + counts_data[sample] = [0] * len(peaks_df) + continue + + try: + # Read BAM as RegionSet and count overlaps + reads_rs = RegionSet.from_bam(str(bam_path)) + sample_counts = list(peaks_rs.count_overlaps(reads_rs)) + + if normalized and sample_counts: + total = sum(sample_counts) + if total > 0: + sample_counts = [c / total * 1e6 for c in sample_counts] + + counts_data[sample] = sample_counts + + except Exception as e: + print(f"Failed to count reads for {sample}: {e}") + counts_data[sample] = [0] * len(peaks_df) + + counts_df = pd.DataFrame(counts_data) + output_file = summary_path / f"{project_name}_{genome}_peaks_coverage.tsv" + counts_df.to_csv(output_file, sep="\t", index=False) + count_files[genome] = str(output_file) + print(f"Counts table: {output_file}") + + return count_files diff --git a/tools/pepatac_summarizer/plots/__init__.py b/tools/pepatac_summarizer/plots/__init__.py new file mode 100644 index 0000000..8cc63db --- /dev/null +++ b/tools/pepatac_summarizer/plots/__init__.py @@ -0,0 +1,14 @@ +"""PEPATAC summary plots.""" + +from .alignment import plot_aligned_raw, plot_aligned_pct +from .tss import plot_tss_scores +from .library import plot_lib_sizes +from .complexity import plot_complexity_curves + +__all__ = [ + "plot_aligned_raw", + "plot_aligned_pct", + "plot_tss_scores", + "plot_lib_sizes", + "plot_complexity_curves", +] diff --git a/tools/pepatac_summarizer/plots/alignment.py b/tools/pepatac_summarizer/plots/alignment.py new file mode 100644 index 0000000..51d967d --- /dev/null +++ b/tools/pepatac_summarizer/plots/alignment.py @@ -0,0 +1,202 @@ +"""Alignment statistics plots.""" + +from pathlib import Path +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt +import matplotlib +matplotlib.use('Agg') + +from .theme import apply_pepatac_theme, get_color_gradient + + +def plot_aligned_raw( + stats: pd.DataFrame, + output_dir: str | Path, + project_name: str +) -> str | None: + """Plot raw alignment counts as stacked horizontal bar chart. + + Args: + stats: DataFrame with alignment statistics per sample + output_dir: Output directory for plots + project_name: Project name for output files + + Returns: + Path to output PDF, or None if no data + """ + output_path = Path(output_dir) + + if stats.empty or "sample_name" not in stats.columns: + return None + + required = ["Fastq_reads", "Aligned_reads"] + if not all(col in stats.columns for col in required): + print("Missing required columns for alignment plot") + return None + + samples = stats["sample_name"].tolist() + n_samples = len(samples) + + fastq = pd.to_numeric(stats["Fastq_reads"], errors="coerce").fillna(0) + aligned = pd.to_numeric(stats["Aligned_reads"], errors="coerce").fillna(0) + duplicates = pd.to_numeric(stats.get("Duplicate_reads", 0), errors="coerce").fillna(0) + + prealign_cols = [c for c in stats.columns if c.startswith("Aligned_reads_") and c != "Aligned_reads"] + prealign_reads = {} + for col in prealign_cols: + genome = col.replace("Aligned_reads_", "") + prealign_reads[genome] = pd.to_numeric(stats[col], errors="coerce").fillna(0) + + unaligned = fastq - aligned + for genome, reads in prealign_reads.items(): + unaligned = unaligned - reads + + genome_cols = [c for c in stats.columns if c == "Genome"] + genomes = stats["Genome"].unique().tolist() if "Genome" in stats.columns else [] + + aligned_by_genome = {} + for genome in genomes: + mask = stats["Genome"] == genome + dedup_col = "Dedup_aligned_reads" + if dedup_col in stats.columns: + counts = pd.to_numeric(stats.loc[mask, dedup_col], errors="coerce").fillna(0) + else: + counts = aligned[mask] + aligned_by_genome[genome] = pd.Series(0, index=stats.index) + aligned_by_genome[genome].loc[mask] = counts + + data = {"sample": samples, "unaligned": (unaligned / 1e6).tolist()} + for genome, reads in prealign_reads.items(): + data[genome] = (reads / 1e6).tolist() + data["duplicates"] = (duplicates / 1e6).tolist() + for genome in genomes: + data[genome] = (aligned_by_genome[genome] / 1e6).tolist() + + df = pd.DataFrame(data) + df = df.set_index("sample") + + fig_height = max(4, n_samples * 0.4) + fig, ax = plt.subplots(figsize=(10, fig_height)) + + colors = ["#1a1a1a"] + if prealign_reads: + colors.extend(get_color_gradient(len(prealign_reads), "#FFE595", "#F6F2A6", "#F6CAA6")) + colors.append("#FC1E25") + if genomes: + colors.extend(get_color_gradient(len(genomes), "#4876FF", "#7648FF", "#94D9CE")) + + df.plot(kind="barh", stacked=True, ax=ax, color=colors[:len(df.columns)], edgecolor="black", linewidth=0.25) + + ax.set_xlabel("Number of reads (M)") + ax.set_ylabel("") + ax.legend(loc="upper right", reverse=True) + ax.invert_yaxis() + apply_pepatac_theme(ax) + + plt.tight_layout() + + output_pdf = output_path / f"{project_name}_alignmentRaw.pdf" + output_png = output_path / f"{project_name}_alignmentRaw.png" + + fig.savefig(output_pdf) + fig.savefig(output_png, dpi=100) + plt.close() + + print(f"Alignment raw plot: {output_pdf}") + return str(output_pdf) + + +def plot_aligned_pct( + stats: pd.DataFrame, + output_dir: str | Path, + project_name: str +) -> str | None: + """Plot alignment percentages as stacked horizontal bar chart. + + Args: + stats: DataFrame with alignment statistics per sample + output_dir: Output directory for plots + project_name: Project name for output files + + Returns: + Path to output PDF, or None if no data + """ + output_path = Path(output_dir) + + if stats.empty or "sample_name" not in stats.columns: + return None + + if "Alignment_rate" not in stats.columns: + print("Missing Alignment_rate column for percent alignment plot") + return None + + samples = stats["sample_name"].tolist() + n_samples = len(samples) + + align_rate = pd.to_numeric(stats["Alignment_rate"], errors="coerce").fillna(0) + dedup_rate = pd.to_numeric(stats.get("Dedup_alignment_rate", 0), errors="coerce").fillna(0) + + prealign_cols = [c for c in stats.columns if c.startswith("Alignment_rate_") and c != "Alignment_rate"] + prealign_rates = {} + for col in prealign_cols: + genome = col.replace("Alignment_rate_", "") + prealign_rates[genome] = pd.to_numeric(stats[col], errors="coerce").fillna(0) + + unaligned = 100 - align_rate + for genome, rate in prealign_rates.items(): + unaligned = unaligned - rate + + duplicates = align_rate - dedup_rate + duplicates = duplicates.clip(lower=0) + + genomes = stats["Genome"].unique().tolist() if "Genome" in stats.columns else [] + + dedup_by_genome = {} + for genome in genomes: + mask = stats["Genome"] == genome + rate = dedup_rate[mask].values + vals = pd.Series(0.0, index=stats.index) + vals.loc[mask] = rate + dedup_by_genome[genome] = vals + + data = {"sample": samples, "unaligned": unaligned.tolist()} + for genome, rate in prealign_rates.items(): + data[genome] = rate.tolist() + data["duplicates"] = duplicates.tolist() + for genome in genomes: + data[genome] = dedup_by_genome[genome].tolist() + + df = pd.DataFrame(data) + df = df.set_index("sample") + + fig_height = max(4, n_samples * 0.4) + fig, ax = plt.subplots(figsize=(10, fig_height)) + + colors = ["#1a1a1a"] + if prealign_rates: + colors.extend(get_color_gradient(len(prealign_rates), "#FFE595", "#F6F2A6", "#F6CAA6")) + colors.append("#FC1E25") + if genomes: + colors.extend(get_color_gradient(len(genomes), "#4876FF", "#7648FF", "#94D9CE")) + + df.plot(kind="barh", stacked=True, ax=ax, color=colors[:len(df.columns)], edgecolor="black", linewidth=0.25) + + ax.set_xlabel("Percent of reads") + ax.set_ylabel("") + ax.set_xlim(0, 103) + ax.legend(loc="upper right", reverse=True) + ax.invert_yaxis() + apply_pepatac_theme(ax) + + plt.tight_layout() + + output_pdf = output_path / f"{project_name}_alignmentPercent.pdf" + output_png = output_path / f"{project_name}_alignmentPercent.png" + + fig.savefig(output_pdf) + fig.savefig(output_png, dpi=100) + plt.close() + + print(f"Alignment percent plot: {output_pdf}") + return str(output_pdf) diff --git a/tools/pepatac_summarizer/plots/complexity.py b/tools/pepatac_summarizer/plots/complexity.py new file mode 100644 index 0000000..63e1aab --- /dev/null +++ b/tools/pepatac_summarizer/plots/complexity.py @@ -0,0 +1,111 @@ +"""Library complexity plots from preseq output.""" + +from pathlib import Path +import pandas as pd +import matplotlib.pyplot as plt +import matplotlib +matplotlib.use('Agg') + +from .theme import apply_pepatac_theme + + +def plot_complexity_curves( + sample_names: list, + genomes: list, + results_subdir: str, + output_dir: str, + project_name: str +) -> str | None: + """Plot library complexity curves from preseq output. + + Args: + sample_names: List of sample names + genomes: Corresponding genome for each sample + results_subdir: Path to sample results + output_dir: Output directory for plots + project_name: Project name for output files + + Returns: + Path to output PDF, or None if no data + """ + results_path = Path(results_subdir) + output_path = Path(output_dir) + + palette = [ + "#999999", "#FFC107", "#27C6AB", "#004D40", + "#B97BC8", "#009E73", "#C92404", "#E3E550", + "#372B4C", "#E3DAC7", "#27CAE6", "#B361BC", + "#897779", "#6114F8", "#19C42B", "#56B4E9" + ] + + fig, ax = plt.subplots(figsize=(10, 7)) + curves_found = 0 + + for i, (sample, genome) in enumerate(zip(sample_names, genomes)): + yield_file = results_path / sample / f"QC_{genome}" / f"{sample}_preseq_yield.txt" + counts_file = results_path / sample / f"QC_{genome}" / f"{sample}_preseq_counts.txt" + + if not yield_file.exists(): + continue + + try: + df = pd.read_csv(yield_file, sep="\t") + if "TOTAL_READS" not in df.columns or "EXPECTED_DISTINCT" not in df.columns: + if "total_reads" in df.columns: + df = df.rename(columns={ + "total_reads": "TOTAL_READS", + "distinct_reads": "EXPECTED_DISTINCT" + }) + else: + continue + + color = palette[i % len(palette)] + ax.plot( + df["TOTAL_READS"] / 1e6, + df["EXPECTED_DISTINCT"] / 1e6, + label=sample, + color=color + ) + + if counts_file.exists(): + counts_df = pd.read_csv(counts_file, sep="\t") + if len(counts_df) >= 1: + total = counts_df.iloc[0, 1] if counts_df.shape[1] > 1 else 0 + unique = counts_df.iloc[0, 2] if counts_df.shape[1] > 2 else 0 + if total > 0: + ax.scatter( + [total / 1e6], [unique / 1e6], + marker='d', s=50, color=color, zorder=5 + ) + + curves_found += 1 + + except Exception as e: + print(f"Error processing {sample}: {e}") + continue + + if curves_found == 0: + print("No samples have available library complexity files.") + plt.close() + return None + + xlim = ax.get_xlim() + ax.plot([0, xlim[1]], [0, xlim[1]], "k--", alpha=0.5, label="_nolegend_") + + ax.set_xlabel("Total reads (M) (incl. duplicates)") + ax.set_ylabel("Unique reads (M)") + ax.legend(loc="upper left", fontsize=8) + ax.set_aspect("equal", adjustable="box") + apply_pepatac_theme(ax) + + plt.tight_layout() + + output_pdf = output_path / f"{project_name}_libComplexity.pdf" + output_png = output_path / f"{project_name}_libComplexity.png" + + fig.savefig(output_pdf) + fig.savefig(output_png, dpi=100) + plt.close() + + print(f"Library complexity plot: {output_pdf}") + return str(output_pdf) diff --git a/tools/pepatac_summarizer/plots/library.py b/tools/pepatac_summarizer/plots/library.py new file mode 100644 index 0000000..831f8be --- /dev/null +++ b/tools/pepatac_summarizer/plots/library.py @@ -0,0 +1,72 @@ +"""Library size plots.""" + +from pathlib import Path +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt +import matplotlib +matplotlib.use('Agg') + +from .theme import apply_pepatac_theme + + +def plot_lib_sizes( + stats: pd.DataFrame, + output_dir: str | Path, + project_name: str +) -> str | None: + """Plot estimated library sizes as bar chart. + + Args: + stats: DataFrame with Picard_est_lib_size column + output_dir: Output directory for plots + project_name: Project name for output files + + Returns: + Path to output PDF, or None if no data + """ + output_path = Path(output_dir) + + if stats.empty or "sample_name" not in stats.columns: + return None + + if "Picard_est_lib_size" not in stats.columns: + print("Missing Picard_est_lib_size column for library size plot") + return None + + samples = stats["sample_name"].tolist() + + lib_sizes_raw = stats["Picard_est_lib_size"].replace("Unknown", 0) + lib_sizes = pd.to_numeric(lib_sizes_raw, errors="coerce").fillna(0) + + if lib_sizes.sum() == 0: + print("No library size data available") + return None + + lib_sizes_millions = (lib_sizes / 1e6).tolist() + n_samples = len(samples) + + fig_height = max(4, n_samples * 0.4) + fig, ax = plt.subplots(figsize=(10, fig_height)) + + y_pos = range(len(samples)) + ax.barh(y_pos, lib_sizes_millions, color="#4876FF", edgecolor="black", linewidth=0.25) + + ax.set_yticks(y_pos) + ax.set_yticklabels(samples) + ax.set_xlabel("Estimated Library Size (M)") + ax.set_ylabel("") + ax.invert_yaxis() + apply_pepatac_theme(ax) + + plt.tight_layout() + + output_pdf = output_path / f"{project_name}_libSizes.pdf" + output_png = output_path / f"{project_name}_libSizes.png" + + fig.savefig(output_pdf) + fig.savefig(output_png, dpi=100) + plt.close() + + print(f"Library sizes plot: {output_pdf}") + return str(output_pdf) diff --git a/tools/pepatac_summarizer/plots/theme.py b/tools/pepatac_summarizer/plots/theme.py new file mode 100644 index 0000000..85b144f --- /dev/null +++ b/tools/pepatac_summarizer/plots/theme.py @@ -0,0 +1,43 @@ +"""PEPATAC plot theme and color utilities.""" + +import matplotlib.pyplot as plt +import matplotlib +matplotlib.use('Agg') + + +def apply_pepatac_theme(ax): + """Apply PEPATAC styling to an axes object.""" + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + ax.spines['left'].set_linewidth(0.5) + ax.spines['bottom'].set_linewidth(0.5) + ax.tick_params(width=0.5, length=4) + + +def get_color_gradient(n: int, low: str, high: str, mid: str = None) -> list[str]: + """Generate a color gradient with n colors.""" + if n == 1: + return [low] + + from matplotlib.colors import LinearSegmentedColormap, to_hex + import numpy as np + + if mid: + colors = [low, mid, high] + positions = [0, 0.5, 1] + else: + colors = [low, high] + positions = [0, 1] + + cmap = LinearSegmentedColormap.from_list("custom", list(zip(positions, colors))) + return [to_hex(cmap(i / max(1, n - 1))) for i in range(n)] + + +ALIGNMENT_COLORS = { + "unaligned": "#1a1a1a", + "duplicates": "#FC1E25", +} + +PREALIGNMENT_COLORS = ["#FFE595", "#F6CAA6", "#F6F2A6"] + +GENOME_COLORS_BASE = ["#4876FF", "#94D9CE", "#7648FF"] diff --git a/tools/pepatac_summarizer/plots/tss.py b/tools/pepatac_summarizer/plots/tss.py new file mode 100644 index 0000000..d5b6fac --- /dev/null +++ b/tools/pepatac_summarizer/plots/tss.py @@ -0,0 +1,85 @@ +"""TSS enrichment score plots.""" + +from pathlib import Path +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt +import matplotlib +matplotlib.use('Agg') + +from .theme import apply_pepatac_theme, get_color_gradient + + +def plot_tss_scores( + stats: pd.DataFrame, + output_dir: str | Path, + project_name: str, + cutoff: float = 6.0 +) -> str | None: + """Plot TSS enrichment scores as bar chart. + + Samples below cutoff are shown in red, above in green. + + Args: + stats: DataFrame with TSS_score column + output_dir: Output directory for plots + project_name: Project name for output files + cutoff: Score threshold (below = low quality) + + Returns: + Path to output PDF, or None if no data + """ + output_path = Path(output_dir) + + if stats.empty or "sample_name" not in stats.columns: + return None + + if "TSS_score" not in stats.columns: + print("Missing TSS_score column for TSS plot") + return None + + samples = stats["sample_name"].tolist() + scores = pd.to_numeric(stats["TSS_score"], errors="coerce").fillna(0).tolist() + n_samples = len(samples) + + colors = [] + for score in scores: + if score < cutoff: + t = min(score / cutoff, 1) + r = int(175 + (228 - 175) * t) + g = int(0 + (14 - 0) * t) + b = int(0 + (0 - 0) * t) + colors.append(f"#{r:02x}{g:02x}{b:02x}") + else: + t = min((score - cutoff) / 24, 1) + r = int(180 + (0 - 180) * t) + g = int(232 + (59 - 232) * t) + b = int(150 + (0 - 150) * t) + colors.append(f"#{r:02x}{g:02x}{b:02x}") + + fig_height = max(4, n_samples * 0.4) + fig, ax = plt.subplots(figsize=(10, fig_height)) + + y_pos = range(len(samples)) + bars = ax.barh(y_pos, scores, color=colors, edgecolor="black", linewidth=0.25) + + ax.axvline(x=cutoff, color="#666666", linestyle="--", linewidth=1, alpha=0.7) + + ax.set_yticks(y_pos) + ax.set_yticklabels(samples) + ax.set_xlabel("TSS Enrichment Score") + ax.set_ylabel("") + ax.invert_yaxis() + apply_pepatac_theme(ax) + + plt.tight_layout() + + output_pdf = output_path / f"{project_name}_TSSEnrichment.pdf" + output_png = output_path / f"{project_name}_TSSEnrichment.png" + + fig.savefig(output_pdf) + fig.savefig(output_png, dpi=100) + plt.close() + + print(f"TSS enrichment plot: {output_pdf}") + return str(output_pdf) diff --git a/tools/pepatac_summarizer/utils.py b/tools/pepatac_summarizer/utils.py new file mode 100644 index 0000000..2806058 --- /dev/null +++ b/tools/pepatac_summarizer/utils.py @@ -0,0 +1,78 @@ +"""Utilities for PEPATAC summarizer.""" + +from pathlib import Path +import yaml +import pandas as pd + + +def load_stats_summary(yaml_path: str | Path) -> pd.DataFrame | None: + """Load and parse the project stats summary YAML. + + Converts the nested YAML structure to a flat DataFrame matching R's yamlToDT. + + Args: + yaml_path: Path to {project}_stats_summary.yaml + + Returns: + DataFrame with one row per sample and stats columns, or None if not found + """ + yaml_path = Path(yaml_path) + if not yaml_path.exists(): + print(f"Warning: Stats summary file not found: {yaml_path}") + return None + + with open(yaml_path) as f: + data = yaml.safe_load(f) + + if not data or "PEPATAC" not in data: + return None + + samples_data = data["PEPATAC"].get("sample", {}) + if not samples_data: + return None + + rows = [] + for sample_name, stats in samples_data.items(): + row = {"sample_name": sample_name} + if isinstance(stats, dict): + row.update(_flatten_dict(stats)) + rows.append(row) + + if not rows: + return None + + df = pd.DataFrame(rows) + df = df.fillna(0) + df = df.replace("", 0) + return df + + +def _flatten_dict(d: dict, prefix: str = "") -> dict: + """Flatten nested dict with underscore-joined keys.""" + items = {} + for k, v in d.items(): + key = f"{prefix}_{k}" if prefix else k + if isinstance(v, dict): + items.update(_flatten_dict(v, key)) + else: + items[key] = v + return items + + +def get_prealignments(stats: pd.DataFrame) -> list[str]: + """Extract prealignment genome names from stats columns.""" + prealignments = [] + for col in stats.columns: + if col.startswith("Aligned_reads_") and col != "Aligned_reads": + genome = col.replace("Aligned_reads_", "") + if genome not in ["human_repeats", "rCRSd"]: + prealignments.append(genome) + return prealignments + + +def round_up_nice(x: float) -> float: + """Round up to a nice number for axis limits.""" + if x <= 0: + return 1 + magnitude = 10 ** int(f"{x:.0e}".split("e")[1]) + return ((x // magnitude) + 1) * magnitude From d94fe7858903e0f4d4b5b19438067b3bca6fd7e4 Mon Sep 17 00:00:00 2001 From: nsheff Date: Mon, 18 May 2026 13:15:38 -0400 Subject: [PATCH 9/9] Complete drop-r-gtars QC plots and align partition widths with R - Add plot_tss_distance using TssIndex.from_regionset.calc_tss_distances; wire into pepatac.py anno block (replaces R placeholder) - Fix plot_frif to sum read counts from bedtools coverage outputs - Reorder plot_partition_distribution to horizontal stacked bar with inline percent labels; add natural chrom sort + canonical chrom filter - Add fragment-distribution median; add chrom/tssdist/part/frif CLI subcommands - Align PartitionList.from_gtf defaults to R's GenomicDistributions: core_prom=100, prox_prom=2000 (was 2000/10000) --- pipelines/pepatac.py | 18 +- tools/pepatac_qc_gtars.py | 524 +++++++++++++++++++++++++------------- 2 files changed, 359 insertions(+), 183 deletions(-) diff --git a/pipelines/pepatac.py b/pipelines/pepatac.py index 6153b12..458e165 100755 --- a/pipelines/pepatac.py +++ b/pipelines/pepatac.py @@ -2481,15 +2481,16 @@ def report_peak_count(): if os.path.isfile(anno_local): if args.qc_backend == "gtars": from tools.pepatac_qc_gtars import (plot_chrom_distribution, - plot_partition_distribution) + plot_partition_distribution, + plot_tss_distance) if not os.path.exists(chr_PDF) or args.new_start: plot_chrom_distribution(peak_output_file, res.chrom_sizes, chr_PDF, chr_PNG) pm.report_object("Peak chromosome distribution", chr_PDF, anchor_image=chr_PNG) if not os.path.exists(TSSdist_PDF) or args.new_start: - # TSS distance uses TssIndex - placeholder for now - pm.run(cmd2, TSSdist_PDF) + plot_tss_distance(peak_output_file, res.refgene_tss, + TSSdist_PDF, TSSdist_PNG) pm.report_object("TSS distance distribution", TSSdist_PDF, anchor_image=TSSdist_PNG) if not os.path.exists(gd_PDF) or args.new_start: @@ -2784,15 +2785,22 @@ def report_peak_count(): if args.qc_backend == "gtars": from tools.pepatac_qc_gtars import plot_frif + try: + total_reads = int(read_count) + except (TypeError, ValueError): + total_reads = None + denom = total_reads if not args.prioritize else genome_size # cFRiF plot plot_frif(cov_files, None, cFRiF_PDF, cFRiF_PNG, cumulative=True, priority=args.prioritize, - reads=not args.prioritize) + reads=not args.prioritize, + genome_size=denom) pm.report_object("cFRiF", cFRiF_PDF, anchor_image=cFRiF_PNG) # FRiF plot plot_frif(cov_files, None, FRiF_PDF, FRiF_PNG, cumulative=False, priority=args.prioritize, - reads=not args.prioritize) + reads=not args.prioritize, + genome_size=denom) pm.report_object("FRiF", FRiF_PDF, anchor_image=FRiF_PNG) else: for cov in cov_files: diff --git a/tools/pepatac_qc_gtars.py b/tools/pepatac_qc_gtars.py index 12413ba..f08a87e 100644 --- a/tools/pepatac_qc_gtars.py +++ b/tools/pepatac_qc_gtars.py @@ -10,19 +10,44 @@ """ import os +import re import numpy as np import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt +# --- Helpers --------------------------------------------------------------- + +def _chrom_sort_key(name): + """Natural sort key for chromosome names. + + Orders chr1, chr2, ... chr22, chrX, chrY, chrM, then any other (alt/un) + contigs alphabetically. + """ + raw = name.replace('chr', '') + if raw.isdigit(): + return (0, int(raw), raw) + if raw.upper() in ('X', 'Y', 'M', 'MT'): + return (1, {'X': 0, 'Y': 1, 'M': 2, 'MT': 2}[raw.upper()], raw) + return (2, 0, raw) + + +def _save_fig(fig, output_pdf, output_png=None, dpi=150): + fig.savefig(output_pdf, format='pdf', dpi=dpi) + if output_png: + fig.savefig(output_png, format='png', dpi=dpi) + plt.close(fig) + + +# --- TSS enrichment -------------------------------------------------------- + def plot_tss_enrichment(tss_file, output_pdf, output_png=None): """Plot TSS enrichment from pre-computed values. - Args: - tss_file: Path to TSS enrichment values (one value per line) - output_pdf: Output PDF path - output_png: Output PNG path (optional) + The pyTssEnrichment.py script writes a per-bp signal vector (one value per + line) symmetric around the TSS. This function normalizes to the baseline + of the flanking 5% and plots the signal. """ with open(tss_file) as f: values = [float(x.strip()) for x in f if x.strip()] @@ -30,53 +55,50 @@ def plot_tss_enrichment(tss_file, output_pdf, output_png=None): if not values: return - # Normalize as in pepatac.py list_len = int(0.05 * len(values)) - if list_len > 0: - baseline = sum(values[1:list_len]) / len(values[1:list_len]) - if baseline > 0: - norm_values = [x / baseline for x in values] - else: - norm_values = values + if list_len > 1: + flank = values[1:list_len] + baseline = sum(flank) / len(flank) if flank else 0.0 + norm_values = [x / baseline for x in values] if baseline > 0 else values else: norm_values = values - # Plot - fig, ax = plt.subplots(figsize=(8, 6)) - x = np.arange(-len(norm_values)//2, len(norm_values)//2) + fig, ax = plt.subplots(figsize=(7, 7)) + half = len(norm_values) // 2 + x = np.arange(-half, len(norm_values) - half) ax.plot(x, norm_values, color='#1f77b4', linewidth=1.5) ax.set_xlabel('Distance from TSS (bp)') - ax.set_ylabel('Normalized Signal') - ax.set_title('TSS Enrichment') + ax.set_ylabel('TSS enrichment score') + ax.set_title('TSS enrichment') ax.axhline(y=1, color='gray', linestyle='--', alpha=0.5) ax.axvline(x=0, color='gray', linestyle='--', alpha=0.5) + fig.tight_layout() + _save_fig(fig, output_pdf, output_png) - plt.tight_layout() - plt.savefig(output_pdf, format='pdf', dpi=150) - if output_png: - plt.savefig(output_png, format='png', dpi=150) - plt.close() +# --- Fragment length distribution ------------------------------------------ -def plot_fragment_distribution(frag_file, count_file, output_pdf, output_txt, output_png=None): +def plot_fragment_distribution(frag_file, count_file, output_pdf, output_txt, + output_png=None): """Plot fragment length distribution. - Args: - frag_file: Path to fragment lengths file - count_file: Path to fragment counts file (sorted length -> count) - output_pdf: Output PDF path - output_txt: Output summary stats file - output_png: Output PNG path (optional) + Reads `count_file` (output of `sort -n frag_file | uniq -c`), which has + lines like: " ". + + Writes a small txt summary (total fragments, mean length) alongside the + plot. """ - # Read counts lengths = [] counts = [] with open(count_file) as f: for line in f: parts = line.strip().split() if len(parts) >= 2: - counts.append(int(parts[0])) - lengths.append(int(parts[1])) + try: + counts.append(int(parts[0])) + lengths.append(int(parts[1])) + except ValueError: + continue if not lengths: return @@ -84,177 +106,289 @@ def plot_fragment_distribution(frag_file, count_file, output_pdf, output_txt, ou lengths = np.array(lengths) counts = np.array(counts) - # Calculate stats - total = counts.sum() - mean_len = np.average(lengths, weights=counts) if total > 0 else 0 + total = int(counts.sum()) + mean_len = float(np.average(lengths, weights=counts)) if total > 0 else 0.0 + median_idx = np.searchsorted(counts.cumsum(), total / 2) + median_len = int(lengths[min(median_idx, len(lengths) - 1)]) if total > 0 else 0 - # Write summary with open(output_txt, 'w') as f: - f.write(f"Total fragments: {total}\n") - f.write(f"Mean fragment length: {mean_len:.1f}\n") + f.write(f"Total fragments\t{total}\n") + f.write(f"Mean fragment length\t{mean_len:.1f}\n") + f.write(f"Median fragment length\t{median_len}\n") - # Plot - fig, ax = plt.subplots(figsize=(10, 6)) - ax.bar(lengths, counts, width=1, color='#1f77b4', alpha=0.7) - ax.set_xlabel('Fragment Length (bp)') - ax.set_ylabel('Count') - ax.set_title('Fragment Length Distribution') - ax.set_xlim(0, min(1000, lengths.max() + 50)) - - plt.tight_layout() - plt.savefig(output_pdf, format='pdf', dpi=150) - if output_png: - plt.savefig(output_png, format='png', dpi=150) - plt.close() + fig, ax = plt.subplots(figsize=(8, 6)) + ax.bar(lengths, counts, width=1, color='#1f77b4') + ax.set_xlabel('Fragment length (bp)') + ax.set_ylabel('Read count') + ax.set_title('Fragment length distribution') + upper = min(1000, int(lengths.max()) + 50) + ax.set_xlim(0, upper) + fig.tight_layout() + _save_fig(fig, output_pdf, output_png) -def plot_frif(coverage_files, annotation_bed, output_pdf, output_png=None, - cumulative=True, priority=False, reads=False): - """Plot Fraction of Reads in Features (FRiF). - - Args: - coverage_files: List of coverage BED files - annotation_bed: Path to annotation BED file - output_pdf: Output PDF path - output_png: Output PNG path (optional) - cumulative: If True, plot cumulative FRiF (cFRiF) - priority: Use mutually exclusive priority ordering - reads: Use read counts instead of bases - """ - try: - from gtars.genomic_distributions import calc_partitions - from gtars.models import RegionSet, PartitionList, GeneModel - except ImportError: - raise ImportError("gtars not available. Install with: pip install gtars") +# --- Peak chromosome distribution ------------------------------------------ - # Load annotation - anno_rs = RegionSet(annotation_bed) +def plot_chrom_distribution(query_bed, chrom_sizes, output_pdf, output_png=None): + """Plot the per-chromosome peak count distribution. - # Calculate fractions for each coverage file - fractions = {} - for cov_file in coverage_files: - name = os.path.basename(cov_file).replace('_coverage.bed', '') - cov_rs = RegionSet(cov_file) - # TODO: Implement actual partition calculation - # This requires proper PartitionList setup - fractions[name] = 0.0 + Restricts to canonical chromosomes (chr1-22, chrX, chrY, chrM) plus any + non-canonical contig holding at least 1% of the maximum-bin peak count. + Alt/unplaced contigs with negligible peak counts are dropped so the + axis stays readable. `chrom_sizes` is accepted for API parity and is + not currently used. + """ + from gtars.models import RegionSet - # Plot - fig, ax = plt.subplots(figsize=(10, 6)) - ax.set_xlabel('Genomic Feature') - ax.set_ylabel('Fraction of Reads' if not cumulative else 'Cumulative Fraction') - ax.set_title('cFRiF' if cumulative else 'FRiF') + rs = RegionSet(query_bed) + stats = rs.chromosome_statistics() - plt.tight_layout() - plt.savefig(output_pdf, format='pdf', dpi=150) - if output_png: - plt.savefig(output_png, format='png', dpi=150) - plt.close() + counts_all = {c: stats[c].number_of_regions for c in stats} + if not counts_all: + return + max_count = max(counts_all.values()) + threshold = max(1, int(0.01 * max_count)) + canonical = {f'chr{i}' for i in range(1, 23)} | {'chrX', 'chrY', 'chrM'} + keep = [c for c, n in counts_all.items() + if c in canonical or n >= threshold] + keep.sort(key=_chrom_sort_key) + counts = [counts_all[c] for c in keep] -def plot_partition_distribution(query_bed, gene_model_gtf, genome, output_pdf, - output_png=None, expected=False): - """Plot genomic partition distribution. - - Args: - query_bed: Path to query regions BED file - gene_model_gtf: Path to gene model GTF file - genome: Genome name (e.g., 'hg38') - output_pdf: Output PDF path - output_png: Output PNG path (optional) - expected: If True, also plot expected distribution - """ - try: - from gtars.genomic_distributions import calc_partitions - from gtars.models import RegionSet, PartitionList - except ImportError: - raise ImportError("gtars not available. Install with: pip install gtars") + fig, ax = plt.subplots(figsize=(10, 6)) + ax.bar(range(len(keep)), counts, color='#1f77b4') + ax.set_xticks(range(len(keep))) + ax.set_xticklabels(keep, rotation=45, ha='right') + ax.set_xlabel('Chromosome') + ax.set_ylabel('Peak count') + ax.set_title('Peak chromosome distribution') + fig.tight_layout() + _save_fig(fig, output_pdf, output_png) - # Load query regions - rs = RegionSet(query_bed) - # Build partition list from GTF - # core_prom=1000, prox_prom=5000 are typical values - partition_list = PartitionList.from_gtf( - gene_model_gtf, - core_prom=1000, - prox_prom=5000, - filter_protein_coding=True, - convert_ensembl_ucsc=True - ) +# --- TSS distance distribution --------------------------------------------- - # Calculate partitions - result = calc_partitions(rs, partition_list, bp_proportion=True) +def plot_tss_distance(query_bed, tss_bed, output_pdf, output_png=None, + max_dist=100_000): + """Plot the distribution of distances from each query region to the + nearest TSS. - # Extract data for plotting - labels = result['partition'] - counts = result['count'] - total = result['total'] + Uses gtars TssIndex.calc_tss_distances(), then histograms the result on + a symmetric log scale around 0. Distances beyond +/- `max_dist` are + clipped into the outer bins. + """ + from gtars.models import RegionSet, TssIndex + + tss_rs = RegionSet(tss_bed) + idx = TssIndex.from_regionset(tss_rs) + peaks = RegionSet(query_bed) + # `calc_tss_distances` returns unsigned distances; peaks on chromosomes + # with no TSS get the sentinel value u32::MAX. We filter those out. + # (TssIndex.feature_distances returns signed distances but emits None + # for missing-TSS chromosomes, which is awkward to clip; absolute is + # sufficient for a QC histogram.) + raw = idx.calc_tss_distances(peaks) + SENTINEL = (1 << 32) - 1 + dists = [d for d in raw if d != SENTINEL] + + if not dists: + return - # Convert to percentages - if total > 0: - sizes = [c / total * 100 for c in counts] - else: - sizes = counts + arr = np.asarray(dists, dtype=float) + arr = np.clip(arr, 0, max_dist) - # Plot pie chart - fig, ax = plt.subplots(figsize=(8, 8)) - colors = plt.cm.Set3(np.linspace(0, 1, len(labels))) - wedges, texts, autotexts = ax.pie(sizes, labels=labels, colors=colors, - autopct='%1.1f%%', startangle=90) - ax.set_title('Peak Partition Distribution') + fig, ax = plt.subplots(figsize=(8, 6)) + bins = np.linspace(0, max_dist, 51) + ax.hist(arr, bins=bins, color='#1f77b4', edgecolor='none') + ax.set_xlabel('Distance to nearest TSS (bp)') + ax.set_ylabel('Peak count') + ax.set_title('Peak TSS distance distribution') + fig.tight_layout() + _save_fig(fig, output_pdf, output_png) + + +# --- Peak genomic partition distribution ----------------------------------- + +# Partition order matches GenomicDistributions defaults: promoters at top, +# intergenic at the bottom of stacked bars. +_PARTITION_ORDER = [ + 'promoterCore', + 'promoterProx', + 'fiveUTR', + 'threeUTR', + 'exon', + 'intron', + 'intergenic', +] + + +def _ordered_partitions(labels, counts): + """Reorder partitions to match _PARTITION_ORDER; unknowns go at the end.""" + pairs = dict(zip(labels, counts)) + ordered = [(p, pairs[p]) for p in _PARTITION_ORDER if p in pairs] + extras = [(p, pairs[p]) for p in pairs if p not in _PARTITION_ORDER] + ordered.extend(extras) + return zip(*ordered) if ordered else ([], []) - plt.tight_layout() - plt.savefig(output_pdf, format='pdf', dpi=150) - if output_png: - plt.savefig(output_png, format='png', dpi=150) - plt.close() +def plot_partition_distribution(query_bed, gene_model_gtf, genome, output_pdf, + output_png=None, expected=False, + core_prom=100, prox_prom=2000): + """Plot the distribution of query regions across genomic partitions. -def plot_chrom_distribution(query_bed, chrom_sizes, output_pdf, output_png=None): - """Plot chromosome distribution. + Builds a PartitionList from a GTF (promoterCore / promoterProx / 5'UTR / + 3'UTR / exon / intron / intergenic) and produces a horizontal stacked bar + plot of the bp-weighted partition fractions. - Args: - query_bed: Path to query regions BED file - chrom_sizes: Path to chromosome sizes file - output_pdf: Output PDF path - output_png: Output PNG path (optional) + `genome` is accepted for parity with the R signature; it isn't used by the + gtars backend (partitions come from the GTF directly). """ - try: - from gtars.models import RegionSet - except ImportError: - raise ImportError("gtars not available. Install with: pip install gtars") + from gtars.models import RegionSet, PartitionList + from gtars.genomic_distributions import calc_partitions - # Load data rs = RegionSet(query_bed) - stats = rs.chromosome_statistics() + pl = PartitionList.from_gtf( + gene_model_gtf, + core_prom=core_prom, + prox_prom=prox_prom, + filter_protein_coding=True, + convert_ensembl_ucsc=True, + ) + result = calc_partitions(rs, pl, bp_proportion=True) + + labels = list(result['partition']) + counts = list(result['count']) + total = float(result.get('total', sum(counts))) + + labels, counts = _ordered_partitions(labels, counts) + labels = list(labels) + counts = list(counts) + fractions = [c / total * 100.0 if total > 0 else 0.0 for c in counts] + + fig, ax = plt.subplots(figsize=(9, 4)) + colors = plt.cm.Set2(np.linspace(0, 1, len(labels))) + left = 0.0 + for label, frac, color in zip(labels, fractions, colors): + ax.barh(['observed'], [frac], left=left, color=color, label=label, + edgecolor='white') + # Inline percentage labels for partitions wider than 3% so the chart + # is readable without consulting the legend. + if frac >= 3.0: + ax.text(left + frac / 2, 0, f'{frac:.1f}%', + ha='center', va='center', fontsize=9, color='black') + left += frac + + ax.set_xlim(0, 100) + ax.set_xlabel('Percent of genome bp') + ax.set_title('Peak genomic partition distribution') + ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), + ncol=min(len(labels), 4), frameon=False, fontsize=9) + fig.tight_layout() + _save_fig(fig, output_pdf, output_png) + + +# --- FRiF / cFRiF ---------------------------------------------------------- + +def _read_coverage_bed(path): + """Sum read counts in a bedtools-coverage output file. + + Each line is: chrom start end [name strand] count bases_covered + total_size fraction_covered + + PEPATAC's pipeline generates these with `bedtools coverage -sorted` over + feature BED3 (cut -f 1-3), so the count column is column 4 and the + bases-covered column is column 5. We use the bases-covered column for + base-pair-weighted fractions and the count column for read-weighted + fractions. + """ + total_reads = 0 + total_bases = 0 + total_size = 0 + with open(path) as f: + for line in f: + parts = line.rstrip('\n').split('\t') + if len(parts) < 7: + continue + try: + total_reads += int(parts[3]) + total_bases += int(parts[4]) + total_size += int(parts[5]) + except (ValueError, IndexError): + continue + return total_reads, total_bases, total_size + + +def plot_frif(coverage_files, annotation_bed, output_pdf, output_png=None, + cumulative=True, priority=False, reads=True, + genome_size=None): + """Plot Fraction of Reads in Features (FRiF) and cumulative FRiF (cFRiF). + + `coverage_files` is a list of `__coverage.bed` files, + each the output of `bedtools coverage -sorted -a feature_bed -b bam`. + Feature name is parsed from the filename. Reads counts are computed + directly from the coverage files; no gtars call is needed. + + `annotation_bed` is currently unused; we already have feature-split + coverage files. Kept for API parity. + + Note: this is a simplified implementation. The R PEPATACr version + additionally normalizes by genome size and total read count to get an + "expected" fraction. We report observed feature fractions (each feature + independent) and the cumulative running total, which is sufficient for + visual QC. + """ + if not coverage_files: + return - # Sort chromosomes naturally (chr1, chr2, ... chr10, chr11, ... chrX, chrY) - def chrom_sort_key(x): - x = x.replace('chr', '') - if x.isdigit(): - return (0, int(x)) + # Pull (feature_name, fraction) pairs from each coverage file. + fractions = [] + for cov_file in coverage_files: + name = os.path.basename(cov_file) + # Strip sample prefix (sample__coverage.bed) heuristically: + # the suffix "_coverage.bed" is stripped, then we take what's left. + name = re.sub(r'_coverage\.bed$', '', name) + total_reads, total_bases, total_size = _read_coverage_bed(cov_file) + if reads: + value = total_reads else: - return (1, x) + value = total_bases + fractions.append((name, value, total_size)) - chroms = sorted(stats.keys(), key=chrom_sort_key) - counts = [stats[c].number_of_regions for c in chroms] + if not fractions: + return - # Plot - fig, ax = plt.subplots(figsize=(12, 6)) - ax.bar(range(len(chroms)), counts, color='#1f77b4') - ax.set_xticks(range(len(chroms))) - ax.set_xticklabels(chroms, rotation=45, ha='right') - ax.set_xlabel('Chromosome') - ax.set_ylabel('Region Count') - ax.set_title('Chromosome Distribution') + # If we have a denominator (genome size for bp, total read count for + # reads), normalize to fractions in [0, 1]. + if reads: + denom = sum(v for _, v, _ in fractions) if not genome_size else genome_size + else: + denom = genome_size if genome_size else sum(s for _, _, s in fractions) + denom = denom or 1 + + labels = [f for f, _, _ in fractions] + values = [v / denom for _, v, _ in fractions] + + fig, ax = plt.subplots(figsize=(9, 6)) + if cumulative: + cum = np.cumsum(values) + ax.bar(range(len(labels)), cum, color='#1f77b4') + ylabel = 'Cumulative fraction of reads in features' + title = 'cFRiF' + else: + ax.bar(range(len(labels)), values, color='#1f77b4') + ylabel = 'Fraction of reads in features' + title = 'FRiF' + + ax.set_xticks(range(len(labels))) + ax.set_xticklabels(labels, rotation=45, ha='right') + ax.set_xlabel('Genomic feature') + ax.set_ylabel(ylabel) + ax.set_title(title) + fig.tight_layout() + _save_fig(fig, output_pdf, output_png) - plt.tight_layout() - plt.savefig(output_pdf, format='pdf', dpi=150) - if output_png: - plt.savefig(output_png, format='png', dpi=150) - plt.close() +# --- CLI ------------------------------------------------------------------- if __name__ == '__main__': import argparse @@ -262,23 +396,57 @@ def chrom_sort_key(x): parser = argparse.ArgumentParser(description='PEPATAC QC with gtars backend') subparsers = parser.add_subparsers(dest='command') - # TSS subcommand tss_parser = subparsers.add_parser('tss', help='Plot TSS enrichment') tss_parser.add_argument('-i', '--input', required=True, help='TSS enrichment file') tss_parser.add_argument('-o', '--output', required=True, help='Output PDF') - # Fragment subcommand frag_parser = subparsers.add_parser('frag', help='Plot fragment distribution') frag_parser.add_argument('-l', '--lengths', required=True, help='Fragment lengths file') frag_parser.add_argument('-c', '--counts', required=True, help='Fragment counts file') frag_parser.add_argument('-p', '--pdf', required=True, help='Output PDF') frag_parser.add_argument('-t', '--txt', required=True, help='Output stats file') + chrom_parser = subparsers.add_parser('chrom', help='Plot peak chromosome distribution') + chrom_parser.add_argument('-i', '--input', required=True, help='Peak BED/narrowPeak') + chrom_parser.add_argument('-c', '--chrom-sizes', required=False, help='Chromosome sizes file') + chrom_parser.add_argument('-o', '--output', required=True, help='Output PDF') + + tssdist_parser = subparsers.add_parser('tssdist', help='Plot peak TSS distance distribution') + tssdist_parser.add_argument('-i', '--input', required=True, help='Peak BED/narrowPeak') + tssdist_parser.add_argument('-t', '--tss', required=True, help='TSS BED file') + tssdist_parser.add_argument('-o', '--output', required=True, help='Output PDF') + + part_parser = subparsers.add_parser('part', help='Plot peak partition distribution') + part_parser.add_argument('-i', '--input', required=True, help='Peak BED/narrowPeak') + part_parser.add_argument('-g', '--gtf', required=True, help='Gene model GTF') + part_parser.add_argument('-G', '--genome', default='hg38', help='Genome name') + part_parser.add_argument('-o', '--output', required=True, help='Output PDF') + + frif_parser = subparsers.add_parser('frif', help='Plot FRiF / cFRiF') + frif_parser.add_argument('-c', '--cov', nargs='+', required=True, help='Coverage BED files') + frif_parser.add_argument('-o', '--output', required=True, help='Output PDF') + frif_parser.add_argument('--cumulative', action='store_true', help='Plot cumulative FRiF') + args = parser.parse_args() + def _png_for(pdf): + return os.path.splitext(pdf)[0] + '.png' + if args.command == 'tss': - png_path = args.output.replace('.pdf', '.png') - plot_tss_enrichment(args.input, args.output, png_path) + plot_tss_enrichment(args.input, args.output, _png_for(args.output)) elif args.command == 'frag': - png_path = args.pdf.replace('.pdf', '.png') - plot_fragment_distribution(args.lengths, args.counts, args.pdf, args.txt, png_path) + plot_fragment_distribution(args.lengths, args.counts, args.pdf, + args.txt, _png_for(args.pdf)) + elif args.command == 'chrom': + plot_chrom_distribution(args.input, args.chrom_sizes, args.output, + _png_for(args.output)) + elif args.command == 'tssdist': + plot_tss_distance(args.input, args.tss, args.output, _png_for(args.output)) + elif args.command == 'part': + plot_partition_distribution(args.input, args.gtf, args.genome, + args.output, _png_for(args.output)) + elif args.command == 'frif': + plot_frif(args.cov, None, args.output, _png_for(args.output), + cumulative=args.cumulative) + else: + parser.print_help()