diff --git a/docs/assets.md b/docs/assets.md index 60a2fa81..a632ea48 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 85eb3fb4..b9ac941f 100755 --- a/docs/detailed-install.md +++ b/docs/detailed-install.md @@ -248,28 +248,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 0163f669..f73caaf6 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-conda.md b/docs/run-conda.md index 53fa07d6..24fef922 100755 --- a/docs/run-conda.md +++ b/docs/run-conda.md @@ -64,28 +64,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 0d244a11..6cd80350 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 diff --git a/findings.md b/findings.md new file mode 100644 index 00000000..4e9064b2 --- /dev/null +++ b/findings.md @@ -0,0 +1,472 @@ +# 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.) + +## 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.) + + +## 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/pipelines/pepatac.py b/pipelines/pepatac.py index 1906e18c..8ef6d943 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 # Make sibling `tools/` importable when this script is invoked directly # (`python pepatac.py ...`), since sys.path[0] is the script's directory @@ -1211,10 +1210,21 @@ def check_alignment_genome(): ngstk.make_dir(QC_folder) bamQC = os.path.join(QC_folder, args.sample_name + "_bamQC.tsv") - cmd = tool_path("bamQC.py") - cmd += " -i " + mapping_genome_bam - cmd += " -c " + str(pm.cores) - cmd += " -o " + bamQC + if args.qc_backend == "gtars": + gtars_cmd_callable = ngstk.check_command("gtars") + if gtars_cmd_callable: + cmd = "gtars uniwig bamqc" + cmd += " -i " + mapping_genome_bam + cmd += " -o " + bamQC + else: + pm.fail_pipeline(RuntimeError( + "Could not call 'gtars'. " + "Confirm the required gtars tool is in your PATH.")) + else: + cmd = tool_path("bamQC.py") + cmd += " -i " + mapping_genome_bam + cmd += " -c " + str(pm.cores) + cmd += " -o " + bamQC def report_bam_qc(bamqc_log): # Reported BAM QC metrics via the bamQC metrics file @@ -1319,7 +1329,7 @@ def post_dup_aligned_reads(dedup_log): rr = float(pm.get_stat("Raw_reads")) tr = float(pm.get_stat("Trimmed_reads")) - if not dr and not dr.strip(): + if not dr or not dr.strip(): pm.info("DEBUG: dr didn't work correctly") dr = ar dr = float(dr) @@ -2592,7 +2602,7 @@ def report_peak_count(): pm.report_object("Peak chromosome distribution", chr_PDF, anchor_image=chr_PNG) if not os.path.exists(TSSdist_PDF) or args.new_start: - if res.refgene_tss and os.path.exists(res.refgene_tss): + if hasattr(res, 'refgene_tss') and res.refgene_tss and os.path.exists(res.refgene_tss): plot_tss_distance(peak_output_file, res.refgene_tss, TSSdist_PDF, TSSdist_PNG) pm.report_object("TSS distance distribution", TSSdist_PDF, diff --git a/pipelines/pepatac.yaml b/pipelines/pepatac.yaml index 713f7e47..0ae381e5 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/requirements-conda.yml b/requirements-conda.yml index c972cc26..62cc5318 100755 --- a/requirements-conda.yml +++ b/requirements-conda.yml @@ -526,8 +526,7 @@ dependencies: - pydantic-core==2.27.2 - pyfaidx==0.8.1.3 - pysam==0.22.1 - - refgenconf==0.13.1 - - 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 80a9d5f8..a927282e 100755 --- a/requirements.txt +++ b/requirements.txt @@ -12,16 +12,17 @@ logmuse>=0.2.7 numpy>=1.21 pandas>=1.4 pararead>=0.7.0 -piper>=0.14.3 +piper>=0.15.1 psutil>=5.8 pysam>=0.16 python-Levenshtein>=0.12 pyyaml>=3.13 -refgenconf>=0.13.1 -# gtars Python binding: used by the --qc-backend gtars path (qc_gtars.py); -# consensus.py no longer uses it. Pinned to 0.9.x to stay aligned with the -# rust-gtars CLI/container and conda env (all 0.9.x). Note: PyPI skipped 0.9.0 -# (the container's version) -- the lowest 0.9.x wheel on PyPI is 0.9.1. +refgenie>=1.0.0 +# gtars Python binding: pinned to 0.9.x to stay aligned with the rust-gtars +# CLI/container and conda env (all 0.9.x), and to carry the 4-arg +# Region(chr, start, end, rest) API the --qc-backend gtars path (qc_gtars.py) +# targets. Note: PyPI skipped 0.9.0 (the container's version) -- the lowest +# 0.9.x wheel on PyPI is 0.9.1. gtars>=0.9.1 # pyBigWig: read per-sample coverage tracks for the reproducible consensus # confidence annotations (signal-centeredness, bimodality) and --recenter. diff --git a/sample_pipeline_interface.yaml b/sample_pipeline_interface.yaml index b816f66d..c2c1b285 100755 --- a/sample_pipeline_interface.yaml +++ b/sample_pipeline_interface.yaml @@ -12,15 +12,15 @@ sample_interface: {% if sample.read2 is defined %} --input2 { sample.read2 } {% endif %} --single-or-paired { sample.read_type } --genome { sample.genome } - {% if sample.chrom_sizes is defined %} --chrom-sizes { sample.chrom_sizes } {% elif sample.genome in refgenie and refgenie[sample.genome].fasta is defined %} --chrom-sizes { refgenie[sample.genome].fasta.chrom_sizes } {% endif %} - {% if sample.TSS_name is defined %} --TSS-name { sample.TSS_name } {% elif sample.genome in refgenie and refgenie[sample.genome].refgene_anno is defined %} --TSS-name { refgenie[sample.genome].refgene_anno.refgene_tss } {% endif %} - {% if sample.blacklist is defined %} --blacklist { sample.blacklist } {% elif sample.genome in refgenie and refgenie[sample.genome].blacklist is defined %} --blacklist { refgenie[sample.genome].blacklist.blacklist } {% endif %} - {% if sample.anno_name is defined %} --anno-name { sample.anno_name } {% elif sample.genome in refgenie and refgenie[sample.genome].feat_annotation is defined %} --anno-name { refgenie[sample.genome].feat_annotation.feat_annotation } {% endif %} + {% if sample.chrom_sizes is defined %} --chrom-sizes { sample.chrom_sizes } {% elif refgenie[sample.genome].fasta is defined %} --chrom-sizes { refgenie[sample.genome].fasta.chrom_sizes } {% endif %} + {% if sample.TSS_name is defined %} --TSS-name { sample.TSS_name } {% elif refgenie[sample.genome].refgene_anno is defined %} --TSS-name { refgenie[sample.genome].refgene_anno.refgene_tss } {% endif %} + {% if sample.blacklist is defined %} --blacklist { sample.blacklist } {% elif refgenie[sample.genome].blacklist is defined %} --blacklist { refgenie[sample.genome].blacklist.blacklist } {% endif %} + {% 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 sample.genome in refgenie and 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 sample.genome in refgenie and 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 %} {% if p in refgenie %} { p ~ '=' ~ refgenie[p].bowtie2_index.dir } {% endif %} {% endfor %} {% else %} --prealignment-index {% for p in sample.prealignment_names %} {% if p in refgenie %} { p ~ '=' ~ refgenie[p].bwa_index.dir } {% endif %} {% 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 %} @@ -29,16 +29,15 @@ sample_interface: {% if sample.frip_ref_peaks is defined %} --frip-ref-peaks { sample.frip_ref_peaks } {% endif %} {% if sample.motif is defined %} --motif {% endif %} {% if sample.sob is defined %} --sob {% endif %} - {% if sample.sob is defined %} {% if sample.genome in refgenie and refgenie[sample.genome].tallymer_index is defined %} --search-file { refgenie[sample.genome].tallymer_index.search_file } {% endif %} {% endif %} - {% if sample.sob is defined %} {% if sample.genome in refgenie and refgenie[sample.genome].fasta is defined %} --fasta { refgenie[sample.genome].fasta.fasta } {% endif %} {% endif %} - {% if sample.fasta is defined %} --fasta { sample.fasta } {% elif sample.genome in refgenie and refgenie[sample.genome].fasta is defined %} --fasta { refgenie[sample.genome].fasta.fasta } {% endif %} + {% if sample.sob is defined %} {% if refgenie[sample.genome].tallymer_index is defined %} --search-file { refgenie[sample.genome].tallymer_index.search_file } {% endif %} {% endif %} + {% if sample.sob is defined %} {% if refgenie[sample.genome].fasta is defined %} --fasta { refgenie[sample.genome].fasta.fasta } {% endif %} {% endif %} + {% if sample.fasta is defined %} --fasta { sample.fasta } {% elif refgenie[sample.genome].fasta is defined %} --fasta { refgenie[sample.genome].fasta.fasta } {% endif %} {% if sample.no_scale is defined %} --no-scale {% endif %} {% if sample.prioritize is defined %} --prioritize {% endif %} {% if sample.keep is defined %} --keep {% endif %} {% if sample.no_fifo is defined %} --noFIFO {% endif %} {% if sample.lite is defined %} --lite {% endif %} {% if sample.skipqc is defined %} --skipqc {% endif %} - {% if sample.skip_dedup is defined %} --skip-dedup {% endif %} --pipestat-config {pipestat.config_file} compute: @@ -49,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 diff --git a/tests/refgenie1_validation/RUN_NOTES.md b/tests/refgenie1_validation/RUN_NOTES.md new file mode 100644 index 00000000..9e5563a9 --- /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` + diff --git a/tools/audit_refgenie_surface.sh b/tools/audit_refgenie_surface.sh new file mode 100755 index 00000000..61b7e513 --- /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' \ + .