v3.0.7 test branch (new_human_filter) merge

.github/workflows/ci.yml

@@ -7,6 +7,7 @@ on:
   pull_request:
     branches: [main]
   workflow_dispatch:
+  workflow_call:              # Reused by publish.yml so the released commit is tested
 
 concurrency:
   # Cancel superseded runs on the same ref to save Actions minutes.

.github/workflows/publish.yml

@@ -0,0 +1,97 @@
+name: Publish to PyPI
+
+# Publishing a GitHub Release (tag vX.Y.Z) tests, builds and uploads the package
+# to PyPI via Trusted Publishing - no API token is stored in the repository.
+# A manual run (workflow_dispatch) tests and builds but does not upload.
+on:
+  release:
+    types: [published]
+  workflow_dispatch:
+
+jobs:
+  # ---------------------------------------------------------------------------
+  # Same unit tests and compile check as ci.yml, run against the released commit.
+  # ---------------------------------------------------------------------------
+  test:
+    uses: ./.github/workflows/ci.yml
+
+  # ---------------------------------------------------------------------------
+  # Build the sdist and wheel, check them, and smoke-test the installed wheel.
+  # ---------------------------------------------------------------------------
+  build:
+    name: Build and check distributions
+    needs: test
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: "3.11"
+
+      - name: Check release tag matches pyproject version
+        if: github.event_name == 'release'
+        run: |
+          pkg_version=$(python -c "import tomllib; print(tomllib.load(open('pyproject.toml','rb'))['project']['version'])")
+          tag_version="${GITHUB_REF_NAME#v}"
+          echo "pyproject: $pkg_version  tag: $tag_version"
+          if [ "$pkg_version" != "$tag_version" ]; then
+            echo "::error::Release tag $GITHUB_REF_NAME does not match pyproject.toml version $pkg_version"
+            exit 1
+          fi
+
+      - name: Build sdist and wheel
+        run: |
+          python -m pip install --upgrade pip build twine
+          python -m build
+          twine check --strict dist/*
+
+      - name: Smoke-test the wheel in a clean environment
+        run: |
+          python -m venv /tmp/wheeltest
+          /tmp/wheeltest/bin/pip install dist/*.whl
+          pkg_version=$(python -c "import tomllib; print(tomllib.load(open('pyproject.toml','rb'))['project']['version'])")
+          /tmp/wheeltest/bin/beegees --version | grep -F "$pkg_version"
+          # Packaged workflow, resources and config must be inside the installed package
+          pkg=$(/tmp/wheeltest/bin/python -c "from beegees.utils.configs import get_package_dir; print(get_package_dir())")
+          for f in workflow/Snakefile \
+                   resources/contaminants/human_reference.fasta \
+                   resources/hmm/COI-5P.hmm \
+                   resources/hmm/rbcL.hmm \
+                   config/config.yaml \
+                   config/multiqc_config.yaml \
+                   profiles/slurm/config.yaml; do
+            test -f "$pkg/$f" || { echo "::error::Missing from wheel: $f"; exit 1; }
+          done
+          # `beegees init` must be able to copy the templates out
+          (cd "$(mktemp -d)" && /tmp/wheeltest/bin/beegees init && test -f config/config.yaml && test -d profiles/slurm)
+
+      - uses: actions/upload-artifact@v4
+        with:
+          name: dist
+          path: dist/
+
+  # ---------------------------------------------------------------------------
+  # Upload to PyPI. Only on a published release, never on a manual run.
+  # Requires a Trusted Publisher on pypi.org: owner bge-barcoding, repo BeeGees,
+  # workflow publish.yml, environment pypi.
+  # ---------------------------------------------------------------------------
+  publish:
+    name: Upload to PyPI
+    needs: build
+    if: github.event_name == 'release'
+    runs-on: ubuntu-latest
+    environment:
+      name: pypi
+      url: https://pypi.org/project/beegees/
+    permissions:
+      id-token: write
+    steps:
+      - uses: actions/download-artifact@v4
+        with:
+          name: dist
+          path: dist/
+
+      - name: Publish to PyPI
+        uses: pypa/gh-action-pypi-publish@release/v1

CHANGELOG

@@ -1,5 +1,58 @@
 # Changelog
 
+## v3.0.7 (2026-10-06)
+- `structural_validation.py` now builds the barcode from every significant, 
+  non-overlapping nhmmer envelope instead of only the lowest-E-value 
+  (i.e. best) one. nhmmer reports one envelope per contiguous stretch matching 
+  the profile, so a consensus with a low-coverage hole in the middle comes back 
+  as two envelopes. `parse_nhmmer_result()` kept one and dropped the rest with 
+  no log line at all, truncating the barcode to whichever fragment happened to 
+  score best. Envelopes are now accepted 'greedily' in E-value order, and only 
+  where they overlap nothing already accepted in HMM and sequence coordinates.
+  The span between envelopes is filled with N.
+- An extra envelope is merged only when it spans more model positions than the
+  gap it would open (i.e. is a "net gain"). The gap has no data and is filled 
+  with N's, so a small, distant envelope adds more ambiguity than sequence. 
+  Envelopes inside the accepted span cost nothing, since they convert
+  existing N into bases, and are always merged. The best-E-value envelope 
+  remains the anchor, so no barcode can come out shorter than before the 
+  envelope net-gain fix.
+- Replaced `01_human_cox1_filter.py` with `01_human_mitogenome_filter.py`, to
+  remove human contamination by mapping reads to a human mitogenome reference 
+  and NUMTs with `bwa aln` instead of comparing them positionally against a 
+  hardcoded COX1 sequence.
+- Added `fasta_cleaner.human_reference` to set the mapping reference. `null` uses
+  the packaged human reference.
+- `fasta_cleaner.human_threshold` is removed - mapping has no similarity
+  threshold.
+- Added the `bwa_index_human_ref` rule, which indexes the human reference once per
+  run and is shared by the merge, concat and se filter rules. The index is
+  written to `resources/bwa_index/` in the working directory via `bwa index -p`. 
+  It is removed at the end of the run; set `keep_bwa_index: true` when using a
+  genome-scale/alternative reference.
+- Renamed the `human_cox1_filter` rules and the `rules.human_cox1_filter` resource
+  key to `human_mitogenome_filter`. Output paths, file names and metric names are
+  unchanged.
+- MultiQC report header no longer shows the absolute path of the analysed
+  directory (`show_analysis_paths: false` in `multiqc_config.yaml`), so server
+  directory layouts (e.g. Galaxy job directories) are not exposed to end users.
+- The `file` column of `structural_validation.csv` now holds paths relative to
+  the output directory (e.g. `03_barcode_recovery/concat_mode/consensus/...`)
+  instead of absolute paths.
+- Renamed the `val_barcode_recovery_metrics_merge` rule to
+  `val_barcode_recovery_metrics_combine` to be more descriptive (not only a 
+  product of the 'merge' pre-processing mode).
+- `gene_fetch.type` is now checked at startup with the other required
+  `gene_fetch` keys, so a missing value fails immediately with a clear message
+  instead of a KeyError when the rule is built.
+- Declared `pyyaml` as a package dependency. The `beegees` CLI imports it to
+  read the config, but `dependencies` was empty, so `pip install beegees` gave a
+  CLI that failed at startup with `ModuleNotFoundError: No module named 'yaml'`.
+  Conda installs were unaffected, since Snakemake pulls in pyyaml.
+- Added `.github/workflows/publish.yml`: publishing a GitHub Release runs the CI
+  tests, builds and checks the sdist and wheel, smoke-tests the installed wheel,
+  and uploads to PyPI via Trusted Publishing.
+
 ## v3.0.6 (2026-09-02)
 - The bioconda recipe was missing `bbmap` (`reformat.sh`, used by the
   optional downsampling rules), `samtools`, `wget` (`download_taxdump`) and

README.md

@@ -141,7 +141,7 @@ beegees run --config config/config.yaml --profile local --log-file beegees.log
 3. **Barcode recovery** - protein reference-guided extraction of barcode sequences from preprocessed reads using MitoGeneExtractor, producing initial consensus sequences.
 4. **Consensus preparation** - header standardisation and concatenation of raw consensus sequences into multi-FASTA.
 5. **Consensus cleaning and filtering (`fasta_cleaner`)** - sequential filters applied to MGE read alignments to remove contaminants and outliers before generating cleaned consensus sequences:
-   - Human COI contamination removal, common in museum specimens ([`01_human_cox1_filter.py`](https://github.com/bge-barcoding/BeeGees/blob/main/workflow/scripts/01_human_cox1_filter.py))
+   - Human contamination removal, common in museum specimens, by mapping reads to a human reference with `bwa aln` ([`01_human_mitogenome_filter.py`](https://github.com/bge-barcoding/BeeGees/blob/main/workflow/scripts/01_human_mitogenome_filter.py)). The packaged reference is the human mitogenome plus human NUMTs
    - AT content filtering, targeting suspected fungal/bacterial contamination ([`02_at_content_filter.py`](https://github.com/bge-barcoding/BeeGees/blob/main/workflow/scripts/02_at_content_filter.py))
    - Statistical outlier removal of reads dissimilar to the initial consensus ([`03_statistical_outlier_filter.py`](https://github.com/bge-barcoding/BeeGees/blob/main/workflow/scripts/03_statistical_outlier_filter.py))
    - Optional custom reference-based filtering ([`04_reference_filter.py`](https://github.com/bge-barcoding/BeeGees/blob/main/workflow/scripts/04_reference_filter.py))
@@ -280,10 +280,10 @@ C: Genetic code for Exonerate (https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprin
 t: Consensus threshold (e.g. 0.5 = 50%) (default: 0.5)
 ```
 
-**fasta_cleaner** - applied in order: (01) human COI → (02) AT content → (03) statistical outlier → (04, optional) reference-based → (05) cleaned consensus generation → (06) metrics aggregation
+**fasta_cleaner** - applied in order: (01) human contamination → (02) AT content → (03) statistical outlier → (04, optional) reference-based → (05) cleaned consensus generation → (06) metrics aggregation
 ```
 consensus_threshold: Proportion of bases at each position that must agree to be included in the consensus (e.g. 0.5 = ≥50%)
-human_threshold: Similarity to the human COI reference above which reads are removed (e.g. 0.95)
+human_reference: FASTA of human sequence that reads are mapped against for contamination removal. null uses the packaged human reference (the human mitogenome NC_012920.1 plus human NUMTs from MANUDB), which removes both mitochondrial and NUMT-derived (nuclear mitochondrial insertion) contamination
 at_difference: AT content deviation from the consensus above which reads are removed (e.g. 0.1 = 10%)
 at_mode: Absolute (remove if AT content differs by more than the threshold in either direction), Higher (remove only above threshold), or Lower (remove only below threshold)
 outlier_percentile: Similarity to the consensus below which reads are flagged as statistical outliers and removed (e.g. 90.0)
@@ -606,5 +606,4 @@ Issues, forks and pull requests are welcome. This pipeline was produced by Dan P
 
 # Future developments #
 - Expand supported markers beyond COI-5P and rbcL. Will require marker-specific HMMs, BLAST databases and associated taxonomy files for barcode validation. Next likely marker to be added = MatK.
-- Update 01_human_cox1_filter.py so it does not solely filter aligned reads against human COI, but instead against the whole human mitogenome.
 - Replacement of final BLASTn taxonomic validation and parsing, using the [BOLDistilled](https://www.boldsystems.org/data/boldistilled/) COI database, with a SINTAX approach instead. SINTAX should be faster, does taxonomic traversal, and provides a confidence score for each match.

beegees/config/config.yaml

@@ -58,11 +58,11 @@ barcode_recovery:
 
 
 ## Post-processing of aligned reads for cleaning (using fasta_cleaner)
-# fasta_cleaner undertakes 6 steps - 1:Human cox1 similarity filtering, 2:AT-content filtering compared to consensus, 3:Statistical outlier filtering compared to consensus, 4:Custom reference-based filtering, 5:Consensus Generation, 6:Metrics aggregation. 
+# fasta_cleaner undertakes 6 steps - 1:Human contamination removal by read mapping, 2:AT-content filtering compared to consensus, 3:Statistical outlier filtering compared to consensus, 4:Custom reference-based filtering, 5:Consensus Generation, 6:Metrics aggregation. 
 # See README.md and docstring of each script (beegees/workflow/scripts) for more detailed parameter information
 fasta_cleaner:
   consensus_threshold: 0.5                       # Threshold at which bases at each position must 'agree' to be incldued in the consensus (0-1)
-  human_threshold: 0.95                          # Threshold at which reads are removed due to similarity with human COI (0-1)
+  human_reference: null                          # FASTA of human sequence that reads are mapped against for contamination removal. null = packaged human mitogenome (NC_012920.1)
   at_difference: 0.1                             # Threshold at which reads are removed due to AT content variation (0-1)
   at_mode: "absolute"                            # At content filtering mode (Absolute/Higher/Lower)
   outlier_percentile: 90.0                       # Threshold at which reads are flagged as statistical outliers to the consensus and removed (0-100)
@@ -92,23 +92,31 @@ taxonomic_validation:
 
 
 ## NCBI taxonomy dump directory
-# Used by the download_taxdump rule and multiqc_plots.
-# Defaults to resources/ncbi_taxdump/ inside your working directory.
+# Defaults to resources/ncbi_taxdump/ inside your BeeGees (working) directory.
 # Override this when running BeeGees as a pip-installed package or when a
 # pre-downloaded dump should be reused across runs.
 # taxdump_dir: "/path/to/writable/ncbi_taxdump"
 
 
-## Resource allocation for each rule
+# The BWA index is built once per run by the bwa_index_human_ref rule and written
+# to resources/bwa_index/ in your BeeGees (working) directory. It is deleted at the end of the
+# run because a small reference re-indexes in about a second. Set the following
+# to true when using a genome-scale reference, where `bwa index` takes an hour or
+# more per gigabase and should be kept between runs.
+# keep_bwa_index: false
 #
-# mem_mb is in MB (2048 = 2 GB) and is multiplied by the attempt number on retry.
+# Override where the index is written:
+# bwa_index_dir: "/path/to/writable/bwa_index"
+
+
+## Resource allocation for each rule
+# mem_mb is in MB (1024 = 1 GB) and is multiplied by the attempt number on retry.
 # A rule that dies at 4096 is resubmitted with 8192, etc. Therefore, size the base 
 # value for a typical run and let retries absorb any outliers.
 #
-# PER-SAMPLE rules run once per sample (MitoGeneExtractor runs once per sample per
-# r/s combination). Their cost is fixed no matter how many samples are run, so
-# size them for your deepest single sample - a 1-sample and a 500-sample run need
-# the same values.
+# PER-SAMPLE rules (generally) run once per sample. 
+# Their cost is fixed no matter how many samples are run, so size them for your deepest 
+# single sample - a 1-sample and a 500-sample run need the same values.
 #
 # AGGREGATE rules run once for the whole run, so their memory grows with sample
 # count. The values below carry headroom to ~500 samples; the per-rule comments
@@ -144,9 +152,9 @@ rules:
     threads: 4                                   # rename_headers.py --threads. >4 is likely to have dimishing speed returns
   gzip_merged_clean:                             # Gzip *_merged_clean.fastq files to reclaim disk space
     mem_mb: 2048                                 # Flat memory requirement regardless of sample number
-  human_cox1_filter:                             # 01_human_cox1_filter.py: remove human COX1 contamination
-    mem_mb: 8192                                 # Process pool over per-sample alignments: peak is [threads] x [largest alignment], so it tracks read depth, NOT sample count
-    threads: 4
+  human_mitogenome_filter:                       # 01_human_mitogenome_filter.py: remove human contamination by mapping reads to a human reference
+    mem_mb: 8192                                 # Process pool over per-sample alignments: peak is [threads] x [largest alignment] plus the BWA index, so it tracks read depth and reference size, NOT sample count
+    threads: 4                                   # Alignment files processed in parallel (one bwa thread each)
   at_content_filter:                             # 02_at_content_filter.py: remove reads with aberrant AT content
     mem_mb: 8192                                 # As above: [threads] x [largest alignment]
     threads: 4

beegees/config/multiqc_config.yaml

@@ -4,6 +4,9 @@
 
 title: "BeeGees Pipeline Report"
 
+# Hide absolute report paths
+show_analysis_paths: false
+
 report_comment: >
   This MultiQC report was generated by the
   <a href="https://github.com/bge-barcoding/BeeGees">BeeGees (Barcode Gene Extraction and Evaluation from Genome Skims) pipeline</a>

beegees/profiles/local/config.yaml

@@ -21,6 +21,4 @@ rerun-incomplete: true       # Re-run jobs whose outputs were left incomplete
 #printshellcmds: true        # Echo the shell command for each job
 
 # --- Retries ---
-restart-times: 3             # Retry a failed job up to 3 times. Note: retries re-request
-                             # the same memory, so an OOM kill will simply fail 3x unless the rule scales mem_mb on `attempt` — e.g.
-                             # resources: mem_mb=lambda wc, attempt: 8000 * attempt
+restart-times: 3             # Retry a failed job up to 3 times.